小波阈值去噪

在对信号进行小波分解时需要选择合适的小波基，由于没有任何一种小波基可以对不同类型的信号达到最优的分解效果，因此，如何选择小波基成为小波分解的一个重点。针对现实中的信号，小波基的选择一般要考虑以下几个因素：支撑长度、对称性、消失矩、正则性、相似性。针对一维信号，例如语音信号，通常选择dB小波和sym小波。

3.2分解层数的选择

在对信号进行小波分解时，分解的层数取得越大，则噪声和信号表现的不同特性越明显，越有利于二者的分离，但是分解的层数越大，经过重构的信号失真也会越大，在一定程度上会对信号去噪的效果产生较差的影响。因此，如何选择分解层数以解决信噪分离效果和重构信号失真之间的矛盾呢？

小波分解的频段范围与采样频率有关。若进行N层分解，则各个频段范围为：

假设原始信号X的采样频率为1000Hz，则信号的最大频率为500，对该信号做3层小波分解，则各个频段范围如下图所示。

3.3阈值的选择

在小波域，有效信号对应的系数很大，而噪声对应的系数很小并且仍满足高斯分布。因此可以通过设定阈值将信号在小波域某段区间内的系数置零，就能最大程度的抑制噪声，同时只是稍微损伤有效信号。

阈值选择规则基于模型 $\large y=f(t)+e$ ， $\large e$ 是高斯白噪声。

目前常见的阈值选择方法有：无偏风险估计阈值、极大极小阈值、固定阈值、启发式阈值。

MATLAB自带的自适应阈值选择函数，调用格式如下：

thr=thselect(X,str);

根据字符串str定义的阈值选择方法求信号X的自适应阈值thr。

str= ‘rigrsure‘：无偏风险估计阈值；

str= ‘minimaxi’：极大极小阈值；

str= ‘sqtwolog’：固定阈值；

str= ‘heursure’：启发式阈值。

一般来讲，极大极小阈值和无偏风险估计阈值比较保守，当噪声在信号的高频段分布较少时，这两种阈值估计方法效果较好，可以将微弱的信号提取出来。而固定阈值和启发式阈值去噪比较彻底，在去噪时显得更为有效，但是也容易把有用的信号误认为噪声去掉。

3.4阈值函数的选择

在确定了高斯白噪声在小波域的阈值门限之后，就需要有个阈值函数对含有噪声系数的小波系数进行过滤，去除高斯噪声系数，常用的阈值函数有硬阈值函数和软阈值函数。

3.4.1硬阈值函数

当小波系数的绝对值大于给定阈值时，小波系数不变；小于阈值时，小波系数置零。

3.4.2软阈值函数

当小波系数的绝对值大于给定阈值时，令小波系数减去阈值；小于阈值时，小波系数置零。

MATLAB自带的阈值处理函数，调用格式如下：

Y=wthresh(X,SORH,thr);

返回输入向量或矩阵X经过软阈值(SORH=’s’)或硬阈值(SORH=’h’)处理后的信号Y，thr是阈值。

4.MATLAB代码

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