模极大值.rar_transform denoise_小波去噪 正弦_模 小波_正弦信号去噪_高斯去噪

%____模极大值法去噪： %二进制小波变换（离散平稳小波变换）并找到模极大序列,对极大值序列处理后重建信号; %参数c, O_d4, O_d3可以调节。 close all; clc; clear; snr=2; %init=2055615866; %[xref,x]=wnoise(2,10,snr,init); f=100;fs=2000 x=sin((0:1999)/fs*(2*pi*f)); x_noise=awgn(x,-5); signal=x_noise; points=2000; level=4; sr=360; num_inter=6; wf='sym1'; %所处理数据的长度 分解的级数 抽样率 迭代次数 小波名称 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);%posw是swd大于零的部分组成的矩阵 pdw=((posw(:,1:points-1)-posw(:,2:points))<0);%pdw是posw的第i列与第i+1列向同行的数据相比较 %比较结果小则为一，相反为零 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; %wpeak1=wpeak;%测试程序用 wpeak(:,1)=wpeak(:,1)+1e-10; wpeak(:,points)=wpeak(:,points)+1e-10; %wpeak2=wpeak; %wpeak12=wpeak1-wpeak2; 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; %此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。 D4_wpeak=wpeak(level,:); M=max(D4_wpeak); Thr=C*M/level; %阈值计算，可参考论文："3mm波段脉冲雷达系统研究和小波去噪分析"。 D4_wpeak=D4_wpeak.*(abs(D4_wpeak)>Thr);%D4_wpeak是D4_wpeak阈值处理后的模极大值点 %模极大值的处理方式： %在尺度j上极大值点位置，构造一个搜索区域， %在尺度j-1中，极大值点落在该区域的点保留，其他的置0； D3_wpeak=wpeak(level-1,:); D4_p=(D4_wpeak~=0); O_d4=3;%该参数确定在上一级搜索极大值的范围，可以调整。 for P_d4=O_d4:(length(D4_wpeak)-O_d4); if D4_p(P_d4)==1;%第四层的模极大值点位置P_d4 for i=1:O_d4-1; D4_p(P_d4-i)=1;%将第四层的模极大值点前面的第1，2个点变为零，确定搜索范围 end ; end; end; D3_wpeak=D3_wpeak.*D4_p;%将第三层在模极大值点之外的点置0 D2_wpeak=wpeak(level-2,:); D3_p=(D3_wpeak~=0); O_d3=3;%该参数确定在上一级搜索极大值的范围，可以调整。 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']; 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); % 计算重建信号 % %____分别计算重建小波以及原信号的信噪比 werr=w2-swd; % % 原信号的小波变换（swd)和重建后的小波变换（w2）的比较 figure, for m=1:level wsnr(m)=20*log10(norm(swd(m,:))/norm(werr(m,:))) subplot(level+1,1,m); plot(swd(m,:)),hold on, plot(w2(m,:),'r');grid on;ylabel(strcat('j=',num2str(m))),axis tight; end err=pswa(1:points)-signal(1:points); snr=20*log10(norm(signal)/norm(err))