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

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

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);

[swa,swd] = swt(signal,level,Lo_D,Hi_D);%swd是细节系数，swa是近似系数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_var>cc*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_mod>w_mod);%(corr_mod>w_mod)返回0或者1

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;s qtwolog;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('滤波误差信号');