基于小波空域相關法去噪MATLAB源程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上close all;clc;clear;snr=5;init=;xref,x=wnoise(1,10,snr,init);signal=x;points=1024;  level=5;  wf='bior 1.5' %sym8,bior 1.5Lo_D,Hi_D,Lo_R,Hi_R=wfilters(wf);swa,swd = swt(signal,level,Lo_D,Hi_D);%swd是細節(jié)系數(shù)，swa是近似系數(shù)Swd_n=swd; swd_org=swd;mask_n=zeros(size(Swd_n);

2、 %先把系數(shù)處理矩陣設置為全0。for j=1:(level-1)     %在1:(level-1)分解層次上對高頻系數(shù)處理，最后一層無法求相關系數(shù)，所以不作處理。    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,

3、:).*swd_org(j+1,:); %定義相關系數(shù)為相鄰兩層的乘積。          cc=1.7; %_用以設定停止迭代的 噪聲能量閾值，需要根據(jù)情況調(diào)節(jié)。_%    while Pw_var>cc*Noise_var    Pw=sum(abs(swd(j,:).2); %計算小波能量    Pcorr=sum(abs(Corr).2); %計算相關系數(shù)能量    Corr_new=Corr.*(Pw/Pcorr)0.5); %歸一化

4、60;        corr_mod=abs(Corr_new);    w_mod=abs(swd(j,:);    swd_n=swd(j,:).*(corr_mod>w_mod);%(corr_mod>w_mod)返回0或者1    swd_n1=(swd_n=0);    mask_n(j,:)=mask_n(j,:)+swd_n1;  %將選出的點賦給系數(shù)處理矩陣相應位置。    swd_n0=ones

5、(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; %將相關系數(shù)選出大值后的地方置0。    Corr=Corr_new;    endend      mask_max=ones

6、(1,length(mask_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;minimaxixcrr=signal_n-xref; %求濾波誤差信號。%畫圖：figure; %空域法處理后的高頻系

7、數(shù)。subplot(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    subpl

8、ot(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);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 1