目 录
1. 小波消噪原理
2. 小波阈值消噪步骤
3. 参数选择
(1)小波基的选择
(2)分解尺度的选择
(3) 阈值的选择
(4)阈值函数的选择
4. 语音消噪中的实例运用
5. 小波消噪matlab代码分析
(1) 函数介绍
(2) 参数选择
(3)SNR择优(重点!)
(4) 择优消噪
(5) 代码实现
1. 小波消噪原理
小波变换是一种信号的时间——尺度(时间——频率)分析方法,它具有多分辨分析的特点,而且在时频两域都具有表征信号局部特征的能力,是一种窗口大小固定不变但其形状可改变,时间窗和频率窗都可以改变的时频局部化分析方法。即在低频部分具有较低的时间分辨率和较高的频率分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率。
傅里叶是将信号分解成一系列不同频率的正余弦函数的叠加,同样小波变换是将信号分解为一系列的小波函数的叠加(或者说不同尺度、时间的小波函数拟合),而这些小波函数都是一个母小波经过平移和尺度伸缩得来的。
小波变换就是把某一被称为基本小波的函数作位移后,再在不同尺度下,与待分析信号作内积,即:
上式中,α>0,成为尺度因子,作用是对基本小波
函数作伸缩,反应位移,其值可正可负,两者都是连续的变量,故又称为连续小波变换。在不同尺度下小波的持续时间随值的加大而增宽,幅度则
与反比减少,但波形不变。(具体详细原理可以参考该视频)
2. 小波阈值消噪步骤
小波消噪的基本步骤为(图1):
(1)分解:选定一种层数为N的小波对信号进行小波分解;
(2)阈值处理:分解后选取一阈值,用阈值函数对各层系数进行量化;
(3)重构:用处理后的系数重构信号。
小波消噪步骤
:X→ca3,cd3,cd2,cd1;小波重构:ca3,cd3,cd2,cd1→X。其中ca为低频信息、近似分量,cd为高频、细节分量。
小波分解示意图
3. 参数选择
小波阈值去噪的基本问题包括四个方面:小波基的选择,分解尺度的选择,阈值的选择,阈值函数的选择。下面我们依次介绍各方面的作用及选择依据。
(1)小波基的选择
小波基经典小波函数主要有haar小波、dbN小波、coifN小波、symlet小波、meyer小波等15中小波。通常我们希望所选取的小波满足以下条件:正交性、高消失矩、紧支性、对称性或反对称性。但事实上具有上述性质的小波是不可能存在的,因为小波是对称或反对称的只有Haar小波,并且高消失矩与紧支性是一对矛盾,所以在应用的时候一般选取具有紧支性的小波以及根据信号的特征来选取较为合适的小波。(本部分参考博文小波滤波小波基、阈值、阈值函数的选择)
由于小波基函数在处理信号时各有特点,且没有任何一种小波基函数可以对所有类型信号都取得最优的去噪效果。一般来讲,db小波系和sym小波系在语音去噪中是经常会被用到的两族小波基。
(2)分解尺度的选择
在小波分解中,分解尺度j的选择也是非常重要的一步。j取得越大,则噪声和信号表现的不同特性越明显,越有利于二者的分离。但另一方面,分解尺度越大,重构到的信号失真也会越大,在一定程度上又会影响最终去噪的效果。因此在应用时要格外注意处理好两者之间的矛盾,选择一个合适的分解尺度。
在语音信号去噪中,由于人的日常发声频率一般在85—1100HZ之间,因此只需保证分解后的最低层频率,在0—250HZ左右即可。例如(重点!!):对于一段采样频率为8000HZ的音频,其所包含的最大语音频率应该为4000HZ,对其进行尺度为4的小波分解后,最底层频率系数为0—250HZ。这时对各层系数进行阈值处理,就已经可以保证信号的噪声得到充分分离了。
(3) 阈值的选择
直接影响去噪效果的一个重要因素就是阈值的选取,不同的阈值选取将有不同的去噪效果。目前常见的阈值选择方法有:固定阈值估计、极值阈值估计、无偏似然估计以及启发式估计等。
一般来讲,极值阈值估计和无偏似然估计方法比较保守,当噪声在信号的高频段分布较少时,这两种阈值估计方法去噪效果较好,可以将微弱的信号提取出来。而固定阈值估计法和启发式阈值估计法去噪比较彻底,在去噪时显得更为有效,但是也容易把有用的高频信号误认为噪声而去除掉。
(4)阈值函数的选择
确定了高斯白噪声在小波系数(域)的阈值门限之后,就需要有阈值函数对这个含有噪声系数的小波系数进行过滤,去除高斯噪声系数,常用的阈值函数有软阈值和硬阈值方法:
- 硬阈值法:当小波系数的绝对值小于给定阈值时,令其为0;大于阈值时,则令其保持不变。
- 软阈值法:当小波系数的绝对值小于给定的阈值时,令其为0;大于阈值时,令其都减去阈值。
硬阈值函数在均方误差意义上优于软阈值法,但是信号会产生附加震荡,产生跳跃点,不具有原始信号的平滑性。
软阈值估计得到的小波系数整体连续性较好,从而使估计信号不会产生附加震荡,但是优于会压缩信号,会产生一定的偏差, 直接影响到重构的信号与真实信号的逼近程度。
4. 语音消噪中的实例运用
综上所述,针对人声特点给出参数选择见下表:
参数选择 | |
小波基 | db小波系 |
分解尺度 | 4 |
阈值 | 启发式阈值选择 |
阈值函数 | 通过SNR择优 |
5. 小波消噪matlab代码分析
(1) 函数介绍
XD = wden(x,tptr,sorh,scal,n,’wname’)
返回信号 x 的去噪版本 XD。该函数使用指定的正交或双正交小波 wname 对 X 进行 N 级小波分解,以获得小波系数。阈值选择规则TPTR应用于小波分解。SORH 和 SCAL 定义了规则的应用方式。
(2) 参数选择
x为需要进行小波分解的对象;tptr为阈值选择标准;sorh为阈值函数选择;scal规定了阈值处理随噪声水平的变化;n为小波分解尺度;wname为使用小波基类型。
其中scal中有三种阈值处理随噪声水平的变化,当scal=one时,阈值处理不随噪声水平变化;当scal=sln时,阈值处理根据第一层小波分解的噪声水平估计进行调整;当scal=min时,阈值处理根据每一层小波分解的噪声水平估计进行调整。
根据第四节的表,对应到代码实现部分我们给出参数设置如下表:
参数设置 | |
tptr阈值选择标准 | heursure 启发式估计 |
sorh 阈值函数 | 通过SNR择优 |
scal阈值处理随噪声水平 | 通过SNR择优 |
n 小波分解尺度 | 4 |
wname 小波基类型 | db小波基() |
(3)SNR择优(重点!)
因为具体应用时无法判断哪种阈值函数更优、哪种阈值处理方法更好,所以需要通过计算信噪比,择优选出最优方法进行消噪。通过计算各种方案消噪后的信噪比择优使用,例如下图
择优消噪
本部分选择了添加高斯随机噪声和单频噪声进行测试。
高斯随机噪声
通过调用matlab中的max函数选取SNR最大的方案进行小波消噪,给出消噪的时频效果图如下:
由上图可见,通过选取最优方案的小波消噪效果较为理想。同时通过SNR比较,该种方案较其他方案不论是信噪比还是均方误差,效果均为突出。
② 单频噪声
由图可见,小波消噪法也能较为明显的滤去加入信号中的单频噪声。但效果相较于针对性的单频陷波器去噪而言,还是差一点。
(5) 代码实现
下面这部分是将未加噪的音频导入到Matlab工作区。
%% 实现小波去噪
%% 载入噪声信号数据
clc;clear all;
[y,fs]=audioread('M3.wav'); %样本数据 和 采样率
y=y(:,1); %取一个声道(如果是双声道)
result={}; %用于存放结果
下面这部分是对数据进行处理:
%% 原始数据准备
N=length(y); %长度
t=(0:N-1)/fs; %对应时间
Y=fft(y); %对信号进行傅氏变换
f=fs/N*(0:round(N/2)-1); %显示实际频点的一半
下面是加噪部分,可选择高斯或单频
%% 加噪
f1=3000;k=0:N-1;
noise=0.1*sin(2*pi*f1*k/fs); %加噪声为3000HZ的正弦信号
yn=y+noise'; %对原始的语音信号进行加噪处理
% noise=0.05*randn(N,1); %使用randn函数产生与音频信号等长度的高斯随机噪声信号
% yn=y+noise; %对原始的语音信号进行加噪处理
Yn=fft(yn,N);
下面是最重要的一部分,各种方式的选择及SNR计算和比较:
%% 硬阈值处理
yh1=wden(yn,'heursure','h','one',4,'db4');%硬阈值去噪处理后的信号序列
yh2=wden(yn,'heursure','h','sln',4,'db4');%硬阈值去噪处理后的信号序列;
yh3=wden(yn,'heursure','h','mln',4,'db4');%硬阈值去噪处理后的信号序列
%% 软阈值处理
ys1=wden(yn,'heursure','s','one',4,'db4');%软阈值去噪处理后的信号序列
ys2=wden(yn,'heursure','s','sln',4,'db4');%软阈值去噪处理后的信号序列
ys3=wden(yn,'heursure','s','mln',4,'db4');%软阈值去噪处理后的信号序列
result(1,:)=[{yh1},{yh2},{yh3},{ys1},{ys2},{ys3}]; %元胞第一行放消噪后信号
%% 计算信噪比SNR
Psig=sum(y.*y);
Pnoih1=sum((y-yh1).*(y-yh1)); %硬阈值 one
Pnoih2=sum((y-yh2).*(y-yh2)); %硬阈值 sln
Pnoih3=sum((y-yh3).*(y-yh3)); %硬阈值 mln
Pnois1=sum((y-ys1).*(y-ys1)); %软阈值 one
Pnois2=sum((y-ys2).*(y-ys2)); %软阈值 sln
Pnois3=sum((y-ys3).*(y-ys3)); %软阈值 mln
SNRh1=10*log10(Psig/Pnoih1);
SNRh2=10*log10(Psig/Pnoih2);
SNRh3=10*log10(Psig/Pnoih3);
SNRs1=10*log10(Psig/Pnois1);
SNRs2=10*log10(Psig/Pnois2);
SNRs3=10*log10(Psig/Pnois3);
result(2,:)=[{SNRh1},{SNRh2},{SNRh3},{SNRs1},{SNRs2},{SNRs3}]; %元胞第二行放各snr
%% 计算均方根误差RMSE
RMSEh1=sqrt(Pnoih1);
RMSEh2=sqrt(Pnoih2);
RMSEh3=sqrt(Pnoih3);
RMSEs1=sqrt(Pnois1);
RMSEs2=sqrt(Pnois2);
RMSEs3=sqrt(Pnois3);
result(3,:)=[{RMSEh1},{RMSEh2},{RMSEh3},{RMSEs1},{RMSEs2},{RMSEs3}];
%% 输出结果
infor={'-----六种阈值设定方式的降噪处理结果-------';...
['硬阈值 one处理的SNR=',num2str(SNRh1),',RMSE=',num2str(RMSEh1)];...
['硬阈值 sln处理的SNR=',num2str(SNRh2),',RMSE=',num2str(RMSEh2)];...
['硬阈值 mln处理的SNR=',num2str(SNRh3),',RMSE=',num2str(RMSEh3)];...
['软阈值 one处理的SNR=',num2str(SNRs1),',RMSE=',num2str(RMSEs1)];...
['软阈值 sln处理的SNR=',num2str(SNRs2),',RMSE=',num2str(RMSEs2)];...
['软阈值 mln处理的SNR=',num2str(SNRs3),',RMSE=',num2str(RMSEs3)]};
disp('-----六种阈值设定方式的降噪处理结果-------');
disp(['硬阈值 one处理的SNR=',num2str(SNRh1),',RMSE=',num2str(RMSEh1)]);
disp(['硬阈值 sln处理的SNR=',num2str(SNRh2),',RMSE=',num2str(RMSEh2)]);
disp(['硬阈值 mln处理的SNR=',num2str(SNRh3),',RMSE=',num2str(RMSEh3)]);
disp(['软阈值 one处理的SNR=',num2str(SNRs1),',RMSE=',num2str(RMSEs1)]);
disp(['软阈值 sln处理的SNR=',num2str(SNRs2),',RMSE=',num2str(RMSEs2)]);
disp(['软阈值 mln处理的SNR=',num2str(SNRs3),',RMSE=',num2str(RMSEs3)]);
最后是绘图部分:
%% 消噪前后频域绘图
[M,U]=max([result{2,:}]) %获取最大SNR的值和列索引
yf=result{1,U};
Yf=fft(yf);
subplot(311);
plot(f,abs(Y(1:round(N/2))));
xlabel('Frequency/Hz');ylabel('Amplitude');
title('加噪前信号频谱');grid;
subplot(312);
plot(f,abs(Yn(1:round(N/2))));
xlabel('Frequency/Hz');ylabel('Amplitude');
title('加噪后信号频谱');grid;
subplot(313);
plot(f,abs(Yf(1:round(N/2))));
xlabel('Frequency/Hz');ylabel('Amplitude');
title('择优消噪后信号频谱');grid;
%% 消噪前后时域绘图
figure(5);
subplot(311);
plot(t,y,'g'); %绘制时域波形
xlabel('Time/s');ylabel('Amplitude'); %添加坐标轴
title('加噪前时域波形');grid;
subplot(312);
plot(t,yn,'g'); %绘制时域波形
xlabel('Time/s');ylabel('Amplitude'); %添加坐标轴
title('加噪后时域波形');grid;
subplot(313);
plot(t,yf,'g'); %绘制时域波形
xlabel('Time/s');ylabel('Amplitude'); %添加坐标轴
title('择优消噪后时域波形');grid;