songjianhui 发表于 2009-4-23 08:13:56

为什么模极大值法无法检测奇异点

原始图形如下,添加了噪声,使用的是db3小波,做了5层分解。
可是使用模极大值法检测不到第二个波的波头,不知道为什么?
请达人指教。
源程序在如下
plot(cable(:,1),cable(:,2))
tcable=cable(:,1);
ycable=cable(:,2);
length(ycable)

signal_length=33856;
t(1:signal_length,1)=0;
Ncable=length(tcable);
for i=1:Ncable
t(i,1)=tcable(i,1);
end

y(1:signal_length,1)=0;
for i=1:Ncable
y(i,1)=ycable(i,1);
end
noise=0.01*rand(1,length(y));
y_noise=y+noise';

signal=y_noise;
%plot(t,signal)
points=length(signal);      level=5;    sr=360;   num_inter=6;   wf='db3';
%所处理数据的长度    分解的级数   抽样率    迭代次数      小波名称
offset=0;


%____进行二进制小波变换(离散平稳小波变换),并给出各级波形:
=wfilters(wf);%Lo_D:分解低通滤波器
                                 %Hi_D:分解高通滤波器
                                 %Lo_R:重建低通滤波器
                                 %Hi_R:重建高通滤波器
                                 
= swt(signal,level,Lo_D,Hi_D);%swa小波概貌
                                        %swd小波系数

figure;%figure1
subplot(level,1,1); plot(real(signal)); grid on;axis tight;
for i=1:level
    subplot(level+1,2,2*(i)+1);
    plot(swa(i,:)); axis tight;grid on;xlabel('time');
    ylabel(strcat('a   ',num2str(i)));
    subplot(level+1,2,2*(i)+2);
    plot(swd(i,:)); axis tight;grid on;
ylabel(strcat('d   ',num2str(i)));
end
%以上内容是对信号进行离散小波变换。
%尺度为5



%____求小波变换的模极大值及其位置,并按级给出小波变换模极大的波形:
% swa:小波概貌;swd:小波细节;
% ddw:局部极大位置; wpeak:小波变换的局部极大序列。
ddw=zeros(size(swd));%建立与swd维数相同的矩阵,矩阵值都为零
pddw=ddw;
nddw=ddw;
posw=swd.*(swd>0);%.*数组乘法,将swd中小于零的数置零
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
negw=swd.*(swd<0);%将swd中大于零的数置零
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
ddw=pddw|nddw;%|:或
ddw(:,1)=1;
ddw(:,points)=1;
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;%第一列的值加上1e-10,不知道出于什么目的,目前看来不会对结果产生影响
wpeak(:,points)=wpeak(:,points)+1e-10;%同理,将最后一列的值加上1e-10,不知道出于什么目的,目前看来同样不会对结果产生影响
%以上内容为:寻找每级尺度上小波变换系数对应的模极大值点

wpeak_threshold=wpeak(level,:)/max(abs(wpeak(level,:)));
threshold=0.1;%阈值

for i=1:points;
   if (abs(wpeak_threshold(i))<threshold)
      wpeak(level,i)=0;
    end
end;
D5_wpeak=wpeak(level,:);

%对最大尺度2.^J上的极大值点设定阈值T0,将低于阈值T0的模极大值点去掉,


%____进行模极大值的处理:
%%C=0.8;
%此参数需要调节,为了在最大尺度上设定合适阈值,以确定最大尺度上该保留的模极大值点。
%%D5_wpeak=wpeak(level,:);
%%M=max(abs(D5_wpeak));
%%Thr=C*M/level; %阈值计算,可参考论文:"3mm波段脉冲雷达系统研究和小波去噪分析"。
%%for i=1:points
%%    if(abs(D5_wpeak(i))<Thr)
%%      wpeak(level,i)=0;
%%    end
%%end
%%D5_wpeak=wpeak(level,:);

figure;%figure2
subplot(level+1,1,1); plot(real(signal)); grid on;axis tight;
for i=1:level
    subplot(level+1,1,i+1);
    plot(wpeak(i,:)); axis tight;grid on;
ylabel(strcat('j=   ',num2str(i)));
end


%步骤4
%模极大值的处理方式:
%在尺度j上极大值点位置,构造一个搜索区域,
%在尺度j-1中,极大值点落在该区域的点保留,其他的置0;
D4_wpeak=wpeak(level-1,:);
sousu=6;
D5_p=(D5_wpeak~=0);
O_d5=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
for P_d5=O_d5:(length(D5_wpeak)-O_d5);
    if D5_p(P_d5)==1;
      for i=1:O_d5-1;
      D5_p(P_d5-i)=1;
      end ;
    end;   
end;
D4_wpeak=D4_wpeak.*D5_p;

D3_wpeak=wpeak(level-2,:);
D4_p=(D4_wpeak~=0);
O_d4=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
for P_d4=O_d4:(length(D4_wpeak)-O_d4);
    if D4_p(P_d4)==1;
      for i=1:O_d4-1;
      D4_p(P_d4-i)=1;
      end ;
    end;   
end;
D3_wpeak=D3_wpeak.*D4_p;

D2_wpeak=wpeak(level-3,:);
D3_p=(D3_wpeak~=0);
O_d3=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
for P_d3=O_d3:(length(D3_wpeak)-O_d3);
    if D3_p(P_d3)==1;
      for i=1:O_d3-1;
      D3_p(P_d3-i)=1;
      end ;
    end;   
end;
D2_wpeak=D2_wpeak.*D3_p;

%第一层单独处理,在第二层极大值点位置上,保留第一层相应极大值点:
D1_wpeak=wpeak(level-4,:);
D2_p=(D2_wpeak~=0);
D1_wpeak=D1_wpeak.*D2_p;

wpeak=;
wpeak=wpeak';

%____重构信号:
pswa=swa(level,:); % pswa: 为待重建的信号
wframe=(wpeak~=0); %迭代初始化
w0=zeros(1,points);
=swt(w0,level,Lo_D,Hi_D);
w2=d;% w2为待重建小波
    for j=1:num_inter
       w2=Py_Pgama(d,wpeak,wframe,1,sr);   % 先进行Py投影和 Pgama投影
       w0=iswt(pswa,w2,Lo_R,Hi_R);         % 再进行Pv投影
       =swt(w0,level,Lo_D,Hi_D);      % Pv
    end
   pswa=iswt(swa(level,:),w2,Lo_R,Hi_R); % 计算重建信号   
   
xcrr=y'-pswa; % 重建误差
figure;%figure3
subplot(411)
plot(y(1:points),'r');
axis tight;
subplot(412)
plot(signal(1:points),'r');
axis tight;
subplot(413)
plot(pswa(1:points));
axis tight;
subplot(414)
plot(xcrr(1:points));
axis tight;
figure%figure4
subplot(611)
plot(signal);
subplot(612)
plot(D5_wpeak);
subplot(613)
plot(D4_wpeak);
subplot(614)
plot(D3_wpeak);
subplot(615)
plot(D2_wpeak);
subplot(616)
plot(D1_wpeak);
=max(abs(D2_wpeak))
=max(abs(D1_wpeak))

xiekai1225 发表于 2009-5-12 14:43:18

先用MATLAB自带的小波工具箱选择不同小波函数试试效果,你这个分解没有产生奇异点啊

阳光zzl 发表于 2010-3-31 20:39:44

可以找到突变点吗?我想了解了解,能不能说一下啊

上帝的孩子 发表于 2010-4-21 10:15:30

回复 3# 阳光zzl


    我也来请教~·

prosavo 发表于 2010-7-21 11:47:20

我也对模极大值感兴趣

primwolf 发表于 2010-8-8 19:00:47

对小波的模极大值,也需要学习的~有些不明白~
页: [1]
查看完整版本: 为什么模极大值法无法检测奇异点

招聘斑竹