空域相關(guān)濾波SSNF算法程序

1、% 空域相關(guān)濾波 % % % %參考文獻：%Wavelet transform domain filters:a spatially selective noise filtration technique%不能選正交小波， 文獻中提出用二次樣條小波 ( quadratic-spline ).此處用雙正交小波： bior 1.5%優(yōu)點：%1.對于信噪比高的信號濾波效果好；%2.對于邊沿的保護強過閾值濾波， 不會產(chǎn)生閾值濾波情況下的過于平滑與Gibbs 現(xiàn)象。%缺點：%1.由于對邊沿信號沒做任何處理，所以邊沿可能會有脈沖噪聲保留下來；%2.計算相關(guān)系數(shù)中， 如果計算出來的小波系數(shù)點位置偏差大，

2、則相關(guān)系數(shù)計算受影響；%3.需要迭代運算， 迭代的噪聲能量閾值選取很重要， 這里以開始段無信號處估計噪聲；%4.需要迭代運算，所以運算量比閾值法大；%5.受分解層次影響，在大尺度上小波系數(shù)點位置偏差更大，相關(guān)系數(shù)計算不準確。%需要具體調(diào)整的地方：%1.分解的尺度；%2.選定用什么信號作為噪聲的估計；%3.設(shè)定停止迭代的噪聲能量閾值參數(shù)cc 。% 空域相關(guān)濾波 % % close all;clc;clear;snr=5;init=2055615866;xref,x=wnoise(1,10,snr,init);signal=x;points=1024; level=5; wf='bior

3、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);% %進行二進制小波變換(離散平穩(wěn)小波變換) ，并給出各級波形： Lo_D,Hi_D,Lo_R,H

4、i_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(s

5、wd(i,:); axis tight;grid on;% ylabel(strcat('d ',num2str(i);% end% %小波系數(shù)的處理：Swd_n=swd;swd_org=swd;mask_n=zeros(size(Swd_n); %先把系數(shù)處理矩陣設(shè)置為全0。for j=1:(level-1)%在 1:(level-1) 分解層次上對高頻系數(shù)處理，最后一層無法求相關(guān)系數(shù)， 所以不作處理。Noise_d1=swd_org(j,:);Noise_d1=Noise_d1(1:80);Noise_var=var(Noise_d1); %以信號的前 80 個只含有噪聲的

6、點估計噪聲在各層的方差。Pw_var=var(swd_org(j,:);Corr=swd_org(j,:).*swd_org(j+1,:); %定義相關(guān)系數(shù)為相鄰兩層的乘積。cc=1.7;%用以設(shè)定停止迭代的噪聲能量閾值，需要根據(jù)情況調(diào)節(jié)。%while Pw_var>cc*Noise_varPw=sum(abs(swd(j,:).A2); %計算小波能量Pcorr=sum(abs(Corr)42); %計算相關(guān)系數(shù)能量Corr_ new=Corr.*(Pw/Pcorr)A0.5); % 歸一化corr_mod=abs(Corr_new);w_mod=abs(swd(j,:);swd_n=

7、swd(j,:).*(corr_mod>w_mod);swd_n1=(swd_n=0);mask_n(j,:)=mask_n(j,:)+swd_n1; %將選出的點賦給系數(shù)處理矩陣相應位置。 swd_n0=ones(size(swd_n1);swd_n0=swd_n0-swd_n1;swd(j,:)=swd(j,:).*swd_n0; %將高頻系數(shù)選出大值后的地方置0。Pw_var=var(swd(j,:);Corr_new=Corr_new.*swd_n0; %將相關(guān)系數(shù)選出大值后的地方置0。Corr=Corr_new;endendmask_max=ones(1,length(mask

8、_n);mask_n=mask_n(1:(level-1),:);mask_max; %最后一層系數(shù)處理矩陣全置 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; % 空域法處理后的高頻系數(shù)。subplo

9、t(level+1,1,1); plot(real(signal); grid on;axis tight;title(' 空域法處理后的高頻系數(shù) ');for i=1:level subplot(level+1,1,i+1); plot(Swd_reg(i,:); axis tight;grid on;ylabel(strcat('j= ',num2str(i);endfigure; % 高頻系數(shù)處理前后的比較。for i=1:level subplot(level,2,2*(i)-1); plot(swd_org(i,:); axis tight;grid

10、on; ylabel(strcat('d ',num2str(i); subplot(level,2,2*(i);plot(Swd_reg(i,:); axis tight;grid on;ylabel(strcat('j= ',num2str(i);endfigure; % 信號濾波前后比較。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); axis tight;grid on