denoisingnew1.m

小波极大值消噪程序

%主程序
%____模极大值法去噪：
%二进制小波变换（离散平稳小波变换）并找到模极大序列,对极大值序列处理后重建信号;
%参数c, O_d4, O_d3可以调节。

clear;
%snr=5;
%init=2055615866;
%[xref,xn]=wnoise(1,10,snr,init);
%load  origindata.mat;
load accel1.txt;
xref=accel1(1:4096,1);   %数据长度应为2的幂次


x = xref; 
b=x;



signal=x;


points=4096;        level=5;    sr=360;   num_inter=6;   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.8; 
%此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%修改添加处
D5_wpeak=wpeak(level,:);
M=max(D5_wpeak);
Thr=C*M/level; %阈值计算，可参考论文："3mm波段脉冲雷达系统研究和小波去噪分析"。
D5_wpeak=D5_wpeak.*(abs(D5_wpeak)>Thr);

%模极大值的处理方式：
%在尺度j上极大值点位置，构造一个搜索区域，
%在尺度j-1中，极大值点落在该区域的点保留，其他的置0；
D4_wpeak=wpeak(level-1,:);
D5_p=(D5_wpeak~=0);
O_d5=10;%该参数确定在上一级搜索极大值的范围，可以调整。
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=10;%该参数确定在上一级搜索极大值的范围，可以调整。
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=10;%该参数确定在上一级搜索极大值的范围，可以调整。
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(1,:);
D2_p=(D2_wpeak~=0);
D1_wpeak=D1_wpeak.*D2_p;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%修改添加处
wpeak=[D1_wpeak' D2_wpeak' D3_wpeak' D4_wpeak' D5_wpeak']; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wpeak=wpeak';


%____重构信号：
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); % 计算重建信号   
     
xcrr=xref-pswa'; % 重建误差    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%修改处－进行转秩处理
figure,
subplot(411)
plot(xref(1:points),'r');
%axis([1 points -2 8]);
subplot(412)
plot(x(1:points),'r');
%axis([1 points -2 8]);
subplot(413)
plot(pswa(1:points)); 
%axis([1 points -2 8]);
subplot(414)
plot(xcrr(1:points)); 
%axis([1 points -2 8]);
all=pswa';
save denoisedaccel1.txt all -ascii
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0=1:length(b);
figure,
subplot(211)
n=length(b);
f1=5000;
for i=1:n
    %t1(i)=i/31.2;
    t2(i)=i/100;
end
plot(t2,b)
subplot(212)
b=hilbert(b);
[tfr,t,f,wt]=tfrscalo(b);
%[pxxx,f1]=psd(b,512,100)
%plot(f1,pxxx)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0=1:length(all);
figure,
subplot(211)
%plot(t0,all)
n=length(all);
f1=5000;
a=all(1:n);
for i=1:n
    %t1(i)=i/31;
    t1(i)=i/100;
end
plot(t1,a)
subplot(212)
a=hilbert(a);
[tfr,t,f,wt]=tfrscalo(a);
%[pxxx1,f2]=psd(a,512,100)
%plot(f2,pxxx1)