标签：系数 阈值 3sigma 噪声 小波 源码 语音 函数
一、简介
通常情况下， 我们在从设备上采集到的信号都是具有一定的噪声的，大多数情况下，可认为这种噪声为高斯白噪声。被噪声污染的信号=干净的信号+噪声。
为什么要使用阈值：由于信号在空间上（或者时间域）是有一定连续性的，因此在小波域，有效信号所产生的小波系数其模值往往较大；而高斯白噪声在空间上（或者时间域）是没有连续性的，因此噪声经过小波变换，在小波阈仍然表现为很强的随机性，通常仍认为是高斯白噪的。那么就得到这样一个结论：在小波域，有效信号对应的系数很大，而噪声对应的系数很小。 刚刚已经说了，噪声在小波域对应的系数仍满足高斯白噪分布。如果在小波域，噪声的小波系数对应的方差为 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
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。