【数字信号去噪】基于matlab小波模极大值数字信号去躁【含Matlab源码 2602期】.zip

N=1024;n1=7; n=0:N-1; Fs=500; t=n/Fs; x=1*sin(2*pi*10*t); y=1.4*rand(1,1024); for i=1:length(y) if y(1,i)>0.96 y(i)=4; else if y(1,i)>0.92 y(1,i)=2; else y(1,i)=0; end end end a=0; for i=1:length(y) if y(1,i)>0 a=a+1; else a=a; end end a signal=x+y; xref=x; points=1024; level=4; 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)); title('含噪声信号'),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; %此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。 D4_wpeak=wpeak(level,:); M=max(D4_wpeak); Thr=C*M/level; %阈值计算，可参考论文："3mm波段脉冲雷达系统研究和小波去噪分析"。 D4_wpeak=D4_wpeak.*(abs(D4_wpeak)>Thr); %模极大值的处理方式： %在尺度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; 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-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投影 level:1层 sr=360采样率 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=signal-pswa; % 重建误差 figure, subplot(411) plot(xref(1:points),'r'); title('无噪声信号') axis([1 points -2 8]); subplot(412) plot(signal(1:points),'r'); title('含噪声信号') axis([1 points -2 8]); subplot(413) plot(pswa(1:points)); title('重建信号') axis([1 points -2 8]); subplot(414) plot(xcrr(1:points)); title('误差') axis([1 points -2 8]); % %____分别计算重建小波以及原信号的信噪比 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 xh=wden(signal,'heursure','h','one',5,'sym7'); xs=wden(signal,'heursure','s','one',5,'sym7'); figure subplot(511) plot(x(1:points),'r'); title('不含噪声信号') axis([1 points -5 5]); subplot(512) plot(signal(1:points),'r'); title('含噪声信号') axis([1 points -5 5]); subplot(513) plot(pswa(1:points),'r'); title('极值去噪') axis([1 points -5 5]); subplot(514) plot(xh(1:points)); title('硬阈值去噪') axis([1 points -5 5]); subplot(515) plot(xs(1:points)); title('软阈值去噪') axis([1 points -5 5]); snr0=20*log10(norm(signal)/norm(y)) snr=20*log10(norm(pswa)/norm(pswa-x)) snr1=20*log10(norm(xh)/norm(xh-x)) snr2=20*log10(norm(xs)/norm(xs-x))