hemitter.m

用于信号的处理

%____模极大值法去噪：
%二进制小波变换（离散平稳小波变换）并找到模极大序列,对极大值序列处理后重建信号;
%参数c, O_d4, O_d3可以调节。
close all;
clc;
clear;
%snr=2;
%init=2055615866;
%[xref,x]=wnoise(2,10,snr,init);
%signal=x;
s1=load ('D:\Program Files\MATLAB\R2006b\work\伸腕.txt');
signal=s1(1:512);
points=512;        level=5;    sr=360;   num_inter=20;   wf='db3';
%所处理数据的长度    分解的级数   抽样率    迭代次数        小波名称
offset=0;


%____进行二进制小波变换（离散平稳小波变换），并给出各级波形：
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);
[swa,swd] = swt(signal,level,Lo_D,Hi_D);
figure;
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


%____求小波变换的模极大值及其位置,并按级给出小波变换模极大的波形:
% swa:小波概貌;  swd:小波细节;
% ddw:局部极大位置; wpeak:小波变换的局部极大序列。
ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
posw=swd.*(swd>0);
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);
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;
wpeak(:,points)=wpeak(:,points)+1e-10;
figure;
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


%____进行模极大值的处理：
%____进行模极大值的处理：
C=0.48; 
%此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。
D5_wpeak=wpeak(level,:);
M=max(D5_wpeak);
Thr=C*M/level; %阈值计算。
absd5=abs(D5_wpeak);
thri=absd5>Thr;
D5_wpeak=D5_wpeak.*thri;

D4_wpeak=wpeak(level-1,:);
D5_wp=find(D5_wpeak);
%D4_p=(D4_wpeak~=0);
D5_p(1:points)=0;
O_d5=3;%该参数确定在上一级搜索极大值的范围，可以调整。
%for P_d4=O_d4:(length(D3_wpeak)-O_d4);
%    if D4_p(P_d4)==1;
for P_d5=3:length(D5_wp)-2;
    D5_p(D5_wp(P_d5))=1; 
        for i=-1*O_d5+1:O_d5-1;
        D5_p(D5_wp(P_d5)-i)=1;
        end ;
     
end;
D4_wpeak=D4_wpeak.*D5_p;
%模极大值的处理方式：
%在尺度j上极大值点位置，构造一个搜索区域，
%在尺度j-1中，极大值点落在该区域的点保留，其他的置0；
D3_wpeak=wpeak(level-2,:);
D4_wp=find(D4_wpeak);
%D4_p=(D4_wpeak~=0);
D4_p(1:points)=0;
O_d4=3;%该参数确定在上一级搜索极大值的范围，可以调整。
%for P_d4=O_d4:(length(D3_wpeak)-O_d4);
%    if D4_p(P_d4)==1;
for P_d4=3:length(D4_wp)-2;
    D4_p(D4_wp(P_d4))=1; 
        for i=-1*O_d4+1:O_d4-1;
        D4_p(D4_wp(P_d4)-i)=1;
        end ;
     
end;
D3_wpeak=D3_wpeak.*D4_p;

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

%第一层单独处理，在第二层极大值点位置上，保留第一层相应极大值点：
%D1_wpeak=wpeak(1,:);
%D2_p=(D2_wpeak~=0);
%D1_wpeak=D1_wpeak.*D2_p;
D1_wpeak=wpeak(level-4,:);
D3_p=(D3_wpeak~=0);
D2_wp=find(D2_wpeak);
O_d2=2;%该参数确定在上一级搜索极大值的范围，可以调整。
D2_p(1:points)=0;
%for P_d3=O_d3:(length(D3_wpeak)-O_d3);
 %   if D3_p(P_d3)==1; 
 %       for i=-1*O_d3+1:O_d3-1;
  %%      D3_p(P_d3-i)=1;
  %      end ;
  %  end;     
%end;
for P_d2=2:length(D2_wp)-2;
    D3_p(D2_wp(P_d2))=0; 
        for i=-1*O_d2+1:O_d2-1;
        D2_p(D2_wp(P_d2)-i)=1;
        end ;
     
end;
D1_wpeak=D1_wpeak.*D2_p;
wpeak=[D1_wpeak' D2_wpeak' D3_wpeak' D4_wpeak' D5_wpeak'];
wpeak=wpeak';


%Hermite插值重构信号
%在较低尺度，保证精确性，采用s=3的5阶插值公式；
pswa = wpeak;
for level_i=1:level;
  pswai=pswa(level_i,:); 
   mmp=find(pswai);
   for l=1:(length(mmp)-1);
       for pswaix=(mmp(l)+1):(mmp(l+1)-1);
 %           pswai(pswaix)=Hermitei3(pswaix,mmp(l),mmp(l+1),mmp(l+2),pswai(mmp(l)),pswai(mmp(l+1)),pswai(mmp(l+2)));
  pswai(pswaix)=Hermitei2(pswaix,mmp(l),mmp(l+1),pswai(mmp(l)),pswai(mmp(l+1)));
        end 
   end
   pswa(level_i,:)=pswai;
end
%在较高尺度，起平滑作用，采用s=2的3阶插值公式；    
%for level_i=3:level;
%   pswai=pswa(level_i,:); 
%   mmp=find(pswai);
%  for l=1:(length(mmp)-1);
%           pswai(pswaix)=Hermitei2(pswaix,mmp(l),mmp(l+1),pswai(mmp(l)),pswai(mmp(l+1)));
%        end
%   end
%   pswa(level_i,:)=pswai;
%end
%小波逆变换
denoise=iswt(swa(level,:),pswa,Lo_R,Hi_R);
%显示，比较
xd = wden(signal,'sqtwolog','s','sln',level,'sym8');%阈值去噪

pswa=swa(level,:); % pswa: 为待重建的信号
wframe=(wpeak~=0); %迭代初始化
w0=zeros(1,points);
[a,d]=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投影
       [a,d]=swt(w0,level,Lo_D,Hi_D);      % Pv
    end
     pswa=iswt(swa(level,:),w2,Lo_R,Hi_R); % 计算重建信号   
figure,
subplot(411)
plot(signal(1:points),'k');
axis([1 points -1.5 1.5]);grid on;
title('伸腕采集信号');xlabel('t/ms');ylabel('s/mv');
subplot(412)
plot(denoise(1:points),'k');
axis([1 points -1.5 1.5]);grid on;
title('小波模极大值去噪Hermite重构结果');xlabel('t/ms');ylabel('s/mv');
subplot(413)
plot(pswa(1:points),'k');
axis([1 points -1.5 1.5]);grid on;
title('小波模极大值去噪交替投影结果');xlabel('t/ms');ylabel('s/mv');
subplot(414)
plot(xd(1:points),'k');
axis([1 points -1.5 1.5]);grid on;
title('小波阈值去噪结果');xlabel('t/ms');ylabel('s/mv');