【语音去噪】基于matlab小波硬阈值语音降噪【含Matlab源码 532期】

一、简介

通常情况下， 我们在从设备上采集到的信号都是具有一定的噪声的，大多数情况下，可认为这种噪声为高斯白噪声。被噪声污染的信号=干净的信号+噪声。

为什么要使用阈值：由于信号在空间上（或者时间域）是有一定连续性的，因此在小波域，有效信号所产生的小波系数其模值往往较大；而高斯白噪声在空间上（或者时间域）是没有连续性的，因此噪声经过小波变换，在小波阈仍然表现为很强的随机性，通常仍认为是高斯白噪的。那么就得到这样一个结论：在小波域，有效信号对应的系数很大，而噪声对应的系数很小。 刚刚已经说了，噪声在小波域对应的系数仍满足高斯白噪分布。如果在小波域，噪声的小波系数对应的方差为 sigma，那么根据高斯分布的特性，绝大部分（99.99%）噪声系数都位于[-3sigma，3sigma]区间内（切比雪夫不等式， 3sigma准则）。因此，只要将区间[-3sigma，3sigma]内的系数置零（这就是常用的硬阈值函数的作用），就能最大程度抑制噪声的，同时只是稍微损伤有效信号。将经过阈值处理后的小波系数重构，就可以得到去噪后的信号。 常用的软阈值函数，是为了解决硬阈值函数 “一刀切” 导致的影响（模小于3sigma的小波系数全部切除，大于3sigma全部保留，势必会在小波域产生突变，导致去噪后结果产生局部的抖动，类似于傅立叶变换中频域的阶跃会在时域产生拖尾）。软阈值函数将模小于 3sigma 的小波系数全部置零，而将模大于 3sigma 的做一个比较特殊的处理，大于 3sigma 的小波系数统一减去 3sigma，小于 -3sigma 的小波系数统一加 3sigma。经过软阈值函数的作用，小波系数在小波域就比较光滑了，因此用软阈值去噪得到的图象看起来很平滑，类似于冬天通过窗户看外面一样，像有层雾罩在图像上似的。

比较硬阈值函数去噪和软阈值函数去噪：硬阈值函数去噪所得到的峰值信噪比（PSNR）较高，但是有局部抖动的现象；软阈值函数去噪所得到的 PSNR 不如硬阈值函数去噪，但是结果看起来很平滑，原因就是软阈值函数对小波系数进行了较大的 “社会主义改造”，小波系数改变很大。因此各种各样的阈值函数就出现了，其目的我认为就是要使大的系数保留，小的系数被剔出，而且在小波域系数过渡要平滑。
如何估计小波域噪声方差 sigma 的估计，这个很简单：把信号做小波变换，在每一个子带利用 robust estimator 估计就可以（可能高频带和低频带的方差不同）。 robust estimator 就是将子带内的小波系数模按大小排列，然后取最中间那个，然后把最中间这个除以0.6745 就得到噪声在某个子带内的方差 sigma。利用这个 sigma，然后选种阈值函数，就可以去去噪了，在 matlab 有实现 api 可使用。

在小波分析中经常用到近似和细节，近似表示信号的高尺度，即低频信息；细节表示信号的低尺度，即高频信息。对含有噪声的信号，噪声分量的主要能量集中在小波解的细节分量中。
在以上过程中，小波基和分解层数的选择，阈值的选取规则，和阈值函数的设计，都是影响最终去噪效果的关键因素。
确定了高斯白噪声在小波系数（域）的阈值门限之后，就需要有个阈值函数对这个含有噪声系数的小波系数进行过滤，去除高斯噪声系数，常用的阈值函数有软阈值和硬阈值方法，很多文献论文中也有在阈值函数进行一些大量的改进和优化。

二、源代码

clear all; clc; close all;

[xx, fs] = wavread('C5_4_y.wav');           % 读入数据文件
xx=xx-mean(xx);                         % 消除直流分量
x=xx/max(abs(xx));                      % 幅值归一化
N=length(x);
%-------------------------加入指定强度的噪声---------------------------------
SNR=5;
s=awgn(x,SNR,'measured','db');               % 叠加噪声
wname='db7';

jN=6;  %分解的层数
snrs=20*log10(norm(x)/norm(s-x));
signal=Wavelet_Hard(s,jN,wname);
signal=signal/max(abs(signal));
snr1=SNR_Calc(x,s);            % 计算初始信噪比
snr2=SNR_Calc(x,signal);            % 计算降噪后的信噪比
snr=snr2-snr1;
fprintf('snr1=%5.4f   snr2=%5.4f   snr=%5.4f\n',snr1,snr2,snr);
% 作图
time=(0:N-1)/fs;                        % 设置时间
subplot 311; plot(time,x,'k'); grid; axis tight;
title('纯语音波形'); ylabel('幅值')
subplot 312; plot(time,s,'k'); grid; axis tight;
%小波硬阈值函数
function signal=Wavelet_Hard(s,jN,wname)
[c,l]=wavedec(s,jN,wname);
%高频分量的索引
first = cumsum(l)+1;
first1=first;
first = first(end-2:-1:1);
ld   = l(end-1:-1:2);
last = first+ld-1;
%--------------------------------------------------------------------------
%硬阈值
cxdhard=c;
for j=1:jN                                      %j是分解尺度
    flk = first(j):last(j);                      %flk是di在c中的索引
    thr(j)=sqrt(2*log((j+1)/j))*median(c(flk))/0.6745;
    for k=0:(length(flk)-1)                 %k是位移尺度
        djk2=c(first(j)+k);
        absdjk=abs(djk2);
        thr1=thr(j);
        %阈值处理
       if absdjk<thr1
           djk2=0;
       end 

三、运行结果

四、备注

版本：2014a

标签：系数,阈值,3sigma,噪声,小波,源码,语音,函数
来源： https://blog.51cto.com/u_15287606/2975944