谱相减MATLAB代码以及信噪比计算

实验二 语音信号的频域处理

一、 实验目的、要求

（1）掌握语音信号频域分析方法 （2）了解语音信号频域的特点

（3）了解谱减法作为频域语音增强的原理与编程实现 （3）了解谱减法的缺点，并分析产生该缺点的原因 二、实验原理

语音虽然是一个时变、非平稳的随机过程。但在短时间内可近似看作是平稳的。因此如果能从带噪语音的短时谱中估计出“纯净”语音的短时谱，即可达到语音增强的目的。由于噪声也是随机过程，因此这种估计只能建立在统计模型基础上。利用人耳感知对语音频谱分量的相位不敏感的特性，这类语音增强算法主要针对短时谱的幅度估计。

短时话幅度估计概述 设一帧加窗后的带噪语音为

()()()

01y n s n d n n N =+≤≤- (2.1)

其中()s n 为纯净语音，()d n 假设为平稳加性高斯噪声。

将()y n 在一组基{()}k n φ上展开，使展对系数为各不相关的随机变量。设()y n 的相关函数为(,)y R n m ，由K －L 展开得知{()}k n φ满足

1

()()(,)()N k y k m K n R n m m λφφ-==∑ (2.2)

则()y n 的展开式为

1

1

0()()()()N k k K N k k n y n Y n

Y y n n φφ-=-=?

=????=??

∑∑ (2.3) 如果()y n 的相关长度小于帧长N ，则()k n φ的近似函数为

2()k nk n j N π???

=

??

? (2.4)

可见()y n 的展开过程实际上相当于离散博里叶交换，其展开系数（为傅里叶变换系数。由()()()y n s n d n =+，则有：k k k Y S N =+。

其中[]||exp k k k Y Y j θ=、[]||exp k k k S S j α=、k N 分别为()y n 、()s n 及()d n 的傅里叶交换系数。由于假设噪声是高斯分布的，其傅里叶系数k N 相当于多个高

斯样本的加权和，故可认为仍然为高斯分布。其均值为0，方差可通过无语音时对噪声的分析而获得。

语音增强的任务就是利用已知的噪声功率谱信息，从k Y 中估计出k S 。由于人耳对相位不敏感，故只需估计出||k S ，然后借用带噪语音的相位，进行傅里叶逆变换就可得到增强的语音。基于短时谱幅度估计的方法的原理图如图2.1所示。

图2.1 短时谱估计原理图

在短时谱幅(STSA)估计基础上，人们提出了许多语音增强算法。基本思想是利用输入带噪语音短时谱幅来估计清洁语音短时谱幅，结合带噪语音相位信息，得到增强信号。

运用短时傅里叶变换(STFT)和重叠相加是短时谱估计技术中最常用的方法。输入信号)(n y 的短时谱幅|)(|k Y 与一个修正因子相乘，得到增强信号谱幅度|)(?|k S

（通常修正因子与噪声信号)(n d 谱幅估计密切相关）；或将含噪语音谱幅减去噪声谱幅估计，得到增强信号谱幅，这些方法统称为相减类型算法(subtractive －type algorithms)。噪声谱幅估计可以通过有声、无声检测获取。一

般认为噪声与信号不相关，增强信号谱幅|)(?|k S

是清洁语音信号)(n s 的谱幅估计。通常假设人耳对含噪语音相位产生的畸变不敏感，所以进行IDFT 恢复成时域信号时，含噪语音的相位一般不作处理。

三、使用仪器、材料

微机（带声卡）、耳机，话筒。

四、 实验步骤

(1) 分析含噪语音信号的频谱(幅度谱和相位谱) (2) 获取噪声信号频谱(幅度谱) (3) 对含噪语音信号进行分帧并进行加窗处理 (4) 将含噪语音信号谱和噪声谱作为输入,进行谱减法 (5) 回复增强语音信号帧. (6) 对比输入信号与增强信号波形,分析算法对其产生的影响. 五、实验过程原始记录（数据，图表，计算）

假设)(n y 为含噪语音离散时间序列，由清洁语音信号)(n s 和非相关加性噪声信号)(n d 组成。)(n y 可表示为

)(n y ＝)(n s ＋)(n d (2.5)

将输入信号按帧处理，前后帧之间重叠（一般为50%），对每帧含噪语音进行加窗处理，然后进行FFT 变换，变换到频率域。含噪语音的能量谱可以表示为

2|)(|k Y ≈2|)(|k S ＋2|)(|k D (2.6) 其中)(k Y 为

)(k Y ＝∑-=-1

02)(N n N

kn j e

n y π＝|)(|k Y )(k i e ? (2.7)

其中)(k ?为含噪语音)(k Y 的相位。

由于没法直接得到含噪语音中噪声能量谱2|)(|k D ，一般将无声阶段的数帧噪声信号进行能量谱加权平均得到噪声能量谱估计2|)(|k D 。假设噪声与语音信号不相关，语音能量谱估计为

2|)(?|k S

＝2|)(|k Y －2|)(?|k D (2.8) 其中清洁语音能量谱估计2|)(?|k S

由含噪语音能量谱减去噪声能量谱估计得到。由于噪声能量谱估计与含噪语音中噪声能量谱之间存在差异，式(3.4)可能出现负值，为了避免能量谱出现负值，将这些负值设为零，这一处理称为半波整流

(half-wave rectification)。通过半波整流，清洁语音能量谱估计2|)(?|k S

可表示为 0

|)(?|0|)(?|0

|)(?||)(?|2222

<>???=k S k S

k S k S

当当 (2.9) 结合含噪语音相位信息，通过逆离散傅里叶变换(IDFT)得到时域清洁语音信

号的估计信号)(?n s

。 )(?n s

＝IDFT()(|)(?|k j e k S ?) (2.10) 其中增强语音信号频谱)(?k S

也可以通过时变滤波器)(k G 重建 )(?k S ＝)(k G )(k Y (2.11) 式中时变滤波器（或称为增益函数）)(k G 可表示为

否则

如果βαβα+<

????

?

?

??

???-=1

|)(||)(?|)|)(||)(?|()|)(||)(?|1()(2

25.0225.02

2k Y k D k Y k D k Y k D

k G (2.12)

六：实验结果，及分析

(a) 原始含噪语音

(b) 谱相减增强语音

(c) 原始含噪语音

(d) 谱相减增强语音

图2.2 谱减法结果分析

实验结果：图2.2中（a），（c）为不同信噪比输入含噪语音信号，图2.2中（b），（d）为对

应增强语音信号。当信噪比较低时，降噪效果较差。并且引入音乐噪声。

算法缺陷分析：

（1）不可避免的引入音乐噪声。要有效地滤除含噪语音中的噪声，需要准确地估计含噪语音中噪声的频谱。噪声谱估计越准确，增强信号谱中音乐噪声越小。然而，由于不能直接得到噪声谱，在绝大多数谱相减算法中，通过加权平均无声阶段噪声谱得到噪声谱估计，这种噪声谱估计与含噪语音中的瞬时噪声谱存在差异，噪声平稳性越差，差异越大，由于这种差异的存在，谱相减不可避免地引入音乐噪声。

(2) 半波、全波整流

式(3.4)中负能量值的产生是由于噪声谱估计发生了错误。这些负值用半波整流（被设定为0）或全波整流（被设定为绝对值），这样处理并没有纠正这种错误，可能导致时域信号更进一步的失真。

(3) 用含噪语音的相位作为增强语音的相位

在生成增强时域信号时，含噪语音的相位并没有作任何修改。这是基于这样一个事实，相位的失真对语音质量下降造成的影响不大。当信噪比比较高（>5dB）时，相位失真确实对语音质量的影响不大，然而，当信噪比较低时（<0dB）这种由于相位失真造成的语音质量下降是可以感觉得到的。

七、实验参考程序代码

谱相减原程序代码

hanning=zeros(1,256);

speech=zeros(1,32000);

dd=wavread('D:\语音信号处理实验2\noise.wav');%噪声；

sp=wavread('D:\语音信号处理实验2\speech.wav');%纯语音；

hanning=zeros(1,256);

noise=dd(10000:45000);

speech=sp(25000:60000);

e2=noise+speech;

szeros=zeros(1,32000);

soutput1=zeros(1,32000);

j=sqrt(-1);

a=3,b=0.01;

Snoise=zeros(1,256);

phase=zeros(1,256);

for n=1:256

hanning(n)=(1/2)*(1-cos((2*pi)*(n-1)/255));

end

% to get noise spectral

Noise=zeros(1,256);Noise1=zeros(1,256);

for i=1:3

Noise1(1:256)=abs(fft(e2(1+i*256:256*(i+1)).*hanning'));

Noise=Noise+Noise1/3;

% spectral subtraction algrithm

for i=1:245 %%第i帧

snoise=e2((i-1)*128+1:(i-1)*128+256).*hanning';

phase=angle(fft(snoise)); % to get noisy signal spectral phase

Snoise=abs(fft(snoise)); % to get noisy signal spectral amplitude

for n=1:256

if(Snoise(n)^2-Noise(n)^2)<0 %the power of enhanced speech must be positive

Sout0(n)=0.1*Snoise(n);

else

Sout0(n)=(Snoise(n)^2-(Noise(n)*1.5)^2)^0.5; %power type spectral subtraction

end

S0(n)=Sout0(n)*(cos(phase(n))+j*sin(phase(n))); % to get enhanced speech spectral amplitude

end

sout0=ifft(S0);

szeros(((i-1)*128+1):((i-1)*128+256))=real(sout0);

soutput1=soutput1+szeros;

szeros=zeros(1,32000);

end

%the following is for SNR calculation

sp_energe=zeros(1,256);sn_energe=zeros(1,256);SN=zeros(1,256);

in_SNR1=zeros(1,125);out_SNR1=zeros(1,125);snoise=zeros(1,256);

for i=1:240

snoise=speech((i-1)*128+1:(i-1)*128+256).*hanning';%第i帧清洁语音存放到snoise

SN=noise((i-1)*128+1:(i-1)*128+256).*hanning'; %第i帧清洁语音存放到SN

%in_SNR1(i)

sp_energe(i)=snoise(1:256)'*snoise(1:256);%第i帧清洁语音snoise的能量存放到sp_energe(i) sn_energe(i)=SN(1:256)'*SN(1:256);%第i帧皂声SN的能量存放到sn_energe(i)

in_SNR1(i)=10*log10(sp_energe(i)/sn_energe(i));%第i帧含噪语音信号输入信噪比in_SNR1(i)

SN=soutput1((i-1)*128+1:(i-1)*128+256).*hanning; %第i帧输出信号存放SN

sn_energe1(i)=SN(1:256)*SN(1:256)';%第i帧输出信号SN能量

out_SNR1(i)=10*log10(sp_energe(i)/abs(sn_energe1(i)-sp_energe(i)));%第i帧增强信号输出信噪比end

plot(e2);

hold on

plot(soutput1,'r');

plot(speech,'g');

plot(in_SNR1) hold on

plot(out_SNR1,'r')

figure

plot(hanning);

%sound(speech); %sound(e2); sound(soutput1);