传统语音增强——基本的维纳滤波语音降噪算法

一、维纳滤波的基本原理
基本维纳滤波就是用来解决从噪声中提取信号问题的一种过滤(或滤波)方法。它基于平稳随机过程模型,且假设退化模型为线性空间不变系统的。实际上这种线性滤波问题,可以看成是一种估计问题或一种线性估计问题。基本的维纳滤波是根据全部过去的和当前的观察数据来估计信号的当前值,它的解是以均方误差最小条件下所得到的系统的传递函数H(z)或单位样本响应h(n)的形式给出的,因此更常称这种系统为最佳线性过滤器或滤波器。设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位样本响应h(n)或传递函数H(z)的表达式,其实质是解维纳-霍夫(Wiener-Hopf)方程。

设带噪语音信号为x(n)=s(n)+v(n)
其中,x(n)表示带噪信号;v(n)表示噪声。则经过维纳滤波器h(n)的输出响应y(n)为

理论上,x(n)通过线性系统h(n)后得到的y(n)应尽量接近于s(n),因此y(n)为s(n)
的估计值,可用 \hat{s}(n) 表示,即

从上式可知,卷积形式可以理解为从当前和过去的观察值x(n),x(n-1),x(n-2)…x(n-m),…来估计信号的当前值\hat{s}(n)。因此,用h(n)进行滤波实际上是一种统计估计问题。

\hat{s}(n)按最小均方误差准则使\hat{s}(n)和s(n)的均方误差ξ=E[e^2(n)]=E[{s(n)-\hat{s}(n)}^2]
达到最小。使ξ最小的充要条件是ξ对于h(n)的偏导数为零,即

上式整理可得

这就是正交性原理或投影原理。可得

已知,s(n)和d(n)是联合宽平稳的,因此令x(n)的自相关函数Rx(m-l)=E{x(n-m)x(n-l)},s(n)与x(n)的互相关函数Rsx(m)=E{s(n)x(n-m)},则上式可变为

上式称为维纳滤波器的标准方程或维纳-霍夫(Wiener-Hopf)方程。如果已知Rsx(m)和Rx(m-l),那么解此方程即可求得维纳滤波器的冲激响应。

当l从0到N-1取有限个整数值时,设滤波器冲激响应序列的长度为N,冲激响应矢量为

滤波器输入数据矢量为

则滤波器的输出为

 这样,上式所示的维纳-霍夫方程可写成Q = Rh
其中,Q=E[x(n)s(n)]是s(n)与x(n)的互相关函数,它是一个N维列矢量;R是x(n)的自相关函数,是N阶方阵R=E[x(n)x^T(n)],则最优的维纳滤波器的冲激响应为hopt=R^(-1)Q

如果进行傅里叶变换可得

式中,Px(k)为x(n)的功率谱密度;Psx(k)为x(n)与s(n)的互功率谱密度。
由于v(n)与s(n)互不相关,即Rsv(k)=0,则可得

此时,上可变为

该式为维纳滤波器的谱估计器。此时,\hat{s}(n)的频谱估计值为

此外,H(k)还可以写为

式中,λs(k)和λv(k)分别为第k个频点上的信号和噪声的功率谱。
传统的维纳滤波法需要估计出纯净语音信号的功率谱,一般用类似谱减法的方法得到,
即用带噪语音功率谱减去估计到的噪声功率谱,这种方法会存在残留噪声大的问题。

二、维纳滤波语音增强实验

基本维纳滤波函数Weina_Norm

名称:Weina_Norm
功能:基本维纳滤波算法。
调用格式:
enhanced = Weina_Norm(x, wind, inc, NIS,alpha,beta)
说明:输入参数x是输入的含噪语音信号;wlen为窗函数或窗长;inc是帧移;NIS是前导无话段帧数;alpha和beta是谱减法的参数。enhanced是降噪后的信号。

函数程序如下:

% 维纳滤波 enhanced=Weina_Norm(x,wind,inc,NIS,alpha,beta)
% MS估计噪声功率谱,需要估计纯净信号功率Ps
% x:输入语音信号
% framesize:帧长
% overlap:帧重叠长度
% NIS:无声帧帧数
% alpha,beta:抑制参数
% ---------------------------------------------------------------------------------------------------------------
 %%
function enhanced=Weina_Norm(x,wind,inc,NIS,alpha,beta)
    nwin=length(wind);           % 取窗长
    if (nwin == 1)              % 判断窗长是否为1,若为1,即表示没有设窗函数
       framesize= wind;               % 是,帧长=win
       wnd=hamming(framesize);                      % 设置窗函数
    else
       framesize = nwin;              % 否,帧长=窗长
       wnd=wind;
    end
    
    y=enframe(x,wnd,inc)';             % 分帧
    framenum=size(y,2);                           % 求帧数
    y_fft = fft(y);                         % FFT
    y_a = abs(y_fft);                       % 求取幅值
    y_phase=angle(y_fft);                   % 求取相位角
    y_a2=y_a.^2;                            % 求能量
    noise=mean(y_a2(:,1:NIS),2);               % 计算噪声段平均能量
    signal=zeros(framesize,1);
    for i=1:framenum
         frame=y(:,i);                                 %取一帧数据
         y_fft=fft(frame);                 %对信号帧y_ham进行短时傅立叶变换,得到频域信号y_fft
         y_fft2=abs(y_fft).^2;     %计算频域信号y_fft每帧的功率谱y_w
     
         %带噪语音谱减去噪声谱
         for k=1:framesize
                if   abs( y_fft2(k) )  >=alpha*noise(k)%(k,i)
                      signal(k)=y_fft2(k)-alpha*noise(k);%(k,i);
                      if signal(k)<0
                          signal(k)=0;
                      end
                else
                      signal(k)=beta*noise(k);%*0.01;
                end
                
         end
         %计算H(W)
         Hw=( signal./(signal+1*noise) ).^1 ;
         %维纳滤波器输出
         yw(:,i)=Hw.*y_fft;
         yt(:,i)=ifft(yw(:,i));
    end
  %采用相位,反而信噪比低
    enhanced=filpframe(yt',wnd,inc);

信噪比计算函数SNR_Calc

名称:SNR_Calc
功能:计算信噪比。
调用格式:
snr=SNR_Calc(x,xn)

说明:输入信号x是输入的纯净语音信号;xn是输入的含噪信号。输出参数snr是计算的信噪比。

函数程序如下:

function snr=SNR_Calc(I,In)
% 计算带噪语音信号的信噪比
% I 是纯语音信号
% In 是带噪的语音信号
% 信噪比计算公式是
% snr=10*log10(Esignal/Enoise)
I=I(:)';                             % 把数据转为一列
In=In(:)';
Ps=sum((I-mean(I)).^2);              % 信号的能量
Pn=sum((I-In).^2);                      % 噪声的能量
snr=10*log10(Ps/Pn);                 % 信号的能量与噪声的能量之比,再求分贝值

 案例、用基本的维纳滤波算法给语音减噪

程序如下:

clear all; clc; close all;

[xx, fs] = audioread('C5_3_y.wav');           % 读入数据文件
xx=xx-mean(xx);                         % 消除直流分量
x=xx/max(abs(xx));                      % 幅值归一化
IS=0.25;                                % 设置前导无话段长度
wlen=200;                               % 设置帧长为25ms
inc=80;                                 % 设置帧移为10ms
SNR=5;                                  % 设置信噪比SNR
NIS=fix((IS*fs-wlen)/inc +1);           % 求前导无话段帧数
alpha=2;
beta=0.01;
signal=awgn(x,SNR,'measured','db');               % 叠加噪声
output=Weina_Norm(x,wlen,inc,NIS,alpha,beta) ;
output=output/max(abs(output));
len=min(length(output),length(x));
x=x(1:len);
signal=signal(1:len);
output=output(1:len);

snr1=SNR_Calc(x,signal);            % 计算初始信噪比
snr2=SNR_Calc(x,output);            % 计算降噪后的信噪比
snr=snr2-snr1;
fprintf('snr1=%5.4f   snr2=%5.4f   snr=%5.4f\n',snr1,snr2,snr);

% 作图
time=(0:len-1)/fs;                        % 设置时间
subplot 311; plot(time,x,'k'); grid; axis tight;
title('纯语音波形'); ylabel('幅值')
subplot 312; plot(time,signal,'k'); grid; axis tight;
title(['带噪语音 信噪比=' num2str(SNR) 'dB']); ylabel('幅值')
subplot 313; plot(time,output,'k');grid;%hold on;
title('滤波后波形'); ylabel('幅值'); xlabel('时间/s');

 运行结果如下:

实验使用到的语音数据下载链接如下:

传统语音增强——最小方均(LMS)自适应滤波算法-数据集文档类资源-CSDN下载

参考文献:语音信号处理实验教程;梁瑞宇、赵力、魏昕(编著) 

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。

相关推荐


学习编程是顺着互联网的发展潮流,是一件好事。新手如何学习编程?其实不难,不过在学习编程之前你得先了解你的目的是什么?这个很重要,因为目的决定你的发展方向、决定你的发展速度。
IT行业是什么工作做什么?IT行业的工作有:产品策划类、页面设计类、前端与移动、开发与测试、营销推广类、数据运营类、运营维护类、游戏相关类等,根据不同的分类下面有细分了不同的岗位。
女生学Java好就业吗?女生适合学Java编程吗?目前有不少女生学习Java开发,但要结合自身的情况,先了解自己适不适合去学习Java,不要盲目的选择不适合自己的Java培训班进行学习。只要肯下功夫钻研,多看、多想、多练
Can’t connect to local MySQL server through socket \'/var/lib/mysql/mysql.sock问题 1.进入mysql路径
oracle基本命令 一、登录操作 1.管理员登录 # 管理员登录 sqlplus / as sysdba 2.普通用户登录
一、背景 因为项目中需要通北京网络,所以需要连vpn,但是服务器有时候会断掉,所以写个shell脚本每五分钟去判断是否连接,于是就有下面的shell脚本。
BETWEEN 操作符选取介于两个值之间的数据范围内的值。这些值可以是数值、文本或者日期。
假如你已经使用过苹果开发者中心上架app,你肯定知道在苹果开发者中心的web界面,无法直接提交ipa文件,而是需要使用第三方工具,将ipa文件上传到构建版本,开...
下面的 SQL 语句指定了两个别名,一个是 name 列的别名,一个是 country 列的别名。**提示:**如果列名称包含空格,要求使用双引号或方括号:
在使用H5混合开发的app打包后,需要将ipa文件上传到appstore进行发布,就需要去苹果开发者中心进行发布。​
+----+--------------+---------------------------+-------+---------+
数组的声明并不是声明一个个单独的变量,比如 number0、number1、...、number99,而是声明一个数组变量,比如 numbers,然后使用 nu...
第一步:到appuploader官网下载辅助工具和iCloud驱动,使用前面创建的AppID登录。
如需删除表中的列,请使用下面的语法(请注意,某些数据库系统不允许这种在数据库表中删除列的方式):
前不久在制作win11pe,制作了一版,1.26GB,太大了,不满意,想再裁剪下,发现这次dism mount正常,commit或discard巨慢,以前都很快...
赛门铁克各个版本概览:https://knowledge.broadcom.com/external/article?legacyId=tech163829
实测Python 3.6.6用pip 21.3.1,再高就报错了,Python 3.10.7用pip 22.3.1是可以的
Broadcom Corporation (博通公司,股票代号AVGO)是全球领先的有线和无线通信半导体公司。其产品实现向家庭、 办公室和移动环境以及在这些环境...
发现个问题,server2016上安装了c4d这些版本,低版本的正常显示窗格,但红色圈出的高版本c4d打开后不显示窗格,
TAT:https://cloud.tencent.com/document/product/1340