小波主要程序去躁

%% 基于小波变换的阈值去噪
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');