小波阈值去躁

%% 基于小波变换的阈值去噪

clc;

clear;

close all;

%% 产生仿真信号

Fs=100; %数据采样率Hz

t=(1:1/Fs:4096*1/Fs)'; %对数据进行采样,将t转置为1列

N = length(t); %数据的采样数目

f1 =0.8; %信号的频率

f2=0.05;

x=2*sin(2*pi*f1*t+cos(2*pi*f2*t)); %产生原始信号

nt=0.9*randn(N,1); %高斯白噪声生成

y=x+nt; %含噪信号

%% 用db4小波对含噪信号进行5层分解并提取系数

[c,l]=wavedec(y,5,'db4'); 

%取第5层低频近似系数

ca5=appcoef(c,l,'db4',5);

%取各层高频细节系数

cd5=detcoef(c,l,5);

cd4=detcoef(c,l,4);

cd3=detcoef(c,l,3);

cd2=detcoef(c,l,2);

cd1=detcoef(c,l,1);

thr=thselect(y,'sqtwolog'); % 阈值获取

%% 进行软阈值处理

ysoft5=wthresh(cd5,'s',thr);

ysoft4=wthresh(cd4,'s',thr);

ysoft3=wthresh(cd3,'s',thr);

ysoft2=wthresh(cd2,'s',thr);

ysoft1=wthresh(cd1,'s',thr);

c1=[ca5;ysoft5;ysoft4;ysoft3;ysoft2;ysoft1];

Y=waverec(c1,l,'db4');

%% 画图

% 时域波形对比图

figure;

subplot(3,1,1);plot(t,x);xlabel('t/s');ylabel('幅值');title('模拟信号x(t)');

subplot(3,1,2);plot(t,y);xlabel('t/s');ylabel('幅值');title('含噪信号y(t)');

subplot(3,1,3);plot(t,Y);xlabel('t/s');ylabel('幅值');title('去噪后信号x(t)');

% 频谱对比图

figure;

fs=5; % 采样频率为5Hz

[f,A] = PinPu(x,fs);

subplot(3,1,1);plot(f,A);xlabel('频率/Hz');ylabel('幅值/V');

[f,A] = PinPu(y,fs);

subplot(3,1,2);plot(f,A);xlabel('频率/Hz');ylabel('幅值/V');

[f,A] = PinPu(Y,fs);

subplot(3,1,3);plot(f,A);xlabel('频率/Hz');ylabel('幅值/V');

%% 降噪指标

% 降噪前信噪比

p1=sum(abs(y).^2)/N;

p2=sum(abs(nt).^2)/N;

SNR(1)=10*log10(p1/p2);

% 降噪后信噪比

p3=sum(abs(Y).^2)/N;

p4=sum(abs(Y-x).^2)/N;

SNR(2)=10*log10(p3/p4);

% 均方根误差

RMSE=sqrt(mean((Y-x).^2));

% 相关系数

CR=corr(Y,x,'type','Pearson');