空域相关滤波SSNF算法程序.doc

%空域相关滤波% %参考文献： %Wavelet transform domain filters:a spatially selective noise filtration technique %不能选正交小波，文献中提出用二次样条小波（quadratic-spline）.此处用双正交小波：bior 1.5 %优点：% 1.对于信噪比高的信号滤波效果好； % 2.对于边沿的保护强过阈值滤波，不会产生阈值滤波情况下的过于平滑与Gibbs现象。 %缺点：% 1.由于对边沿信号没做任何处理，所以边沿可能会有脉冲噪声保留下来； % 2.计算相关系数中，如果计算出来的小波系数点位置偏差大，则相关系数计算受影响； % 3.需要迭代运算，迭代的噪声能量阈值选取很重要，这里以开始段无信号处估计噪声； % 4.需要迭代运算，所以运算量比阈值法大； % 5.受分解层次影响，在大尺度上小波系数点位置偏差更大，相关系数计算不准确。 %需要具体调整的地方：% 1.分解的尺度； % 2.选定用什么信号作为噪声的估计； % 3.设定停止迭代的噪声能量阈值参数cc。 %空域相关滤波% close all; clc; clear; snr=5; init=2055615866; xref,x=wnoise(1,10,snr,init); signal=x; points=1024; level=5; wf=bior 1.5; %sym8,bior 1.5 % % %阈值消噪： S_thr=wden(signal,rigrsure,s,mln,5,sym8); % % rigrsure;heursure;sqtwolog; minimaxi % %one;sln;mln % subplot(211);plot(signal); % title(阈值滤波); % subplot(212);plot(S_thr); % %进行二进制小波变换（离散平稳小波变换），并给出各级波形： 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 % %小波系数的处理： Swd_n=swd; swd_org=swd; mask_n=zeros(size(Swd_n); %先把系数处理矩阵设置为全0。 for j=1:(level-1) %在1:(level-1)分解层次上对高频系数处理，最后一层无法求相关系数，所以不作处理。 Noise_d1=swd_org(j,:); Noise_d1=Noise_d1(1:80); Noise_var=var(Noise_d1); %以信号的前80个只含有噪声的点估计噪声在各层的方差。 Pw_var=var(swd_org(j,:); Corr=swd_org(j,:).*swd_org(j+1,:); %定义相关系数为相邻两层的乘积。 cc=1.7; %用以设定停止迭代的噪声能量阈值，需要根据情况调节。% while Pw_varcc*Noise_var Pw=sum(abs(swd(j,:).2); %计算小波能量 Pcorr=sum(abs(Corr).2); %计算相关系数能量 Corr_new=Corr.*(Pw/Pcorr)0.5); %归一化 corr_mod=abs(Corr_new); w_mod=abs(swd(j,:); swd_n=swd(j,:).*(corr_modw_mod); swd_n1=(swd_n=0); mask_n(j,:)=mask_n(j,:)+swd_n1; %将选出的点赋给系数处理矩阵相应位置。 swd_n0=ones(size(swd_n1); swd_n0=swd_n0-swd_n1; swd(j,:)=swd(j,:).*swd_n0; %将高频系数选出大值后的地方置0。 Pw_var=var(swd(j,:); Corr_new=Corr_new.*swd_n0; %将相关系数选出大值后的地方置0。 Corr=Corr_new; end end mask_max=ones(1,length(mask_n); mask_n=mask_n(1:(level-1),:);mask_max; %最后一层系数处理矩阵全置1。 Swd_reg=swd_org.*mask_n; signal_n=iswt(swa,Swd_reg,wf); % S_mix=wden(signal_n,sqtwolog,s,sln,5,sym8); %rigrsure;heursure;sqtwolog;minimaxi xcrr=signal_n-xref; %求滤波误差信号。 % %画图： figure; %空域法处理后的高频系数。 subplot(level+1,1,1); plot(real(signal); grid on;axis tight; title(空域法处理后的高频系数); for i=1:level subplot(level+1,1,i+1); plot(Swd_reg(i,:); axis tight;grid on; ylabel(strcat(j= ,num2str(i); end figure; %高频系数处理前后的比较。 for i=1:level subplot(level,2,2*(i)-1); plot(swd_org(i,:); axis tight;grid on; ylabel(strcat(d ,num2str(i); subplot(level,2,2*(i); plot(Swd_reg(i,:); axis tight;grid on; ylabel(strcat(j= ,num2str(i); end figure; %信号滤波前后比较。 subplot(3,1,1); plot(signal); axis tight;grid on; axis(0 1000 -5 22); title(原始信号); subplot(3,1,2); plot(signal_n); axis tight;grid on; axis(0 1000 -5 22); title(空域法滤波后信号); subplot(3,1,3); plot(xcrr); axis tight;grid on; axis(0 1000 -5 22); title(滤波误差信号); figure; %空域法滤波与阈值滤波的比较。 subplot(3,1,1); plot(signal