www.pudn.com > MMSE-denoising.rar > MMSE.m
% 本程序为经典MMSE方法，引自Y. Ephraim and D.Malah “Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator ,”
%IEEE trans.Acoust.Speech Signal Processing, Vol.ASSP-32,no.6, pp. 1109-1121,1984
%1984年，Y. Ephraim and D.Malah提出的MMSE方法是基于语音短时谱估计的增强方法中的代表算法，一直被称为well kown，可见其重要程度。
%从1984年至今，他们一直在语音增强领域作着不懈的努力，有着突出的贡献，成为语音增强领域的著名人物
%Y. Ephraim在IEEE上及相关刊物发表有关语音增强文章七十余篇。1992年，受IEEE约稿，Ephraim发表了一篇长达四十多页，参考文献300多篇的文章，全面的总结语音增强发展的过程
%二人的研究成果和个人简介均可可从网上得到
%******************初始化部分****************************************%
clear;
d=wavread('E:\学位论文写作\音源\原始\white__05dB.wav'); %读入带噪语音,要求为8Kz采样
d=d';
le=length(d); %求输入带噪语音序列长度
s=fix(le/64); %帧间重叠3/4,所以算出所需循环次数
hn=hanning(256)'; %得到STFT分析窗,汉宁窗,帧长256点
hm(1:256)=hanning(256)'; %得到STFT合成窗,汉宁窗,帧长256点
xd(1:256)=0; %噪声谱估计赋初值
mn(1:256)=0;
x0(1:s*64)=0; %增强后语音赋初值
snpr1=0.98; %snpr1代表前一帧的先验信噪比,这是初始化设为0.98，也可以设为别的值，经过更新之后，会达到较真实的值
%****************估计噪声*****************************************%
%假设带噪语音开始部分为纯噪声
for n=1:1:20 %用前5帧带噪语音进行噪声估计，周期图法
d1(1:256)=d((n-1)*64+1:(n-1)*64+256).*hn;
mn=abs(fft(d1,256)); %本算法是基于短时谱的语音增强，所以不考虑噪声相位
mn=mn.*mn;
xd=xd+mn;
end
xd=xd./20; %对几帧的噪声谱进行了平均，求出噪声的幅度谱
%****************主体增强部分****************************************%
for n=1:1:s-3
d1(1:256)=d((n-1)*64+1:(n-1)*64+256).*hn; %时域加窗
Xd(1:256)=fft(d1,256); %FFT
Man=abs(Xd); %带噪语音幅度
En=angle(Xd); %带噪语音相位
snpo=Man.^2./xd; %计算后验信噪比，利用"D-D"法
snpr=(0.98.*snpr1+0.02.*max(snpo-1,0))*1.1;%计算先验信噪比，实验结果表明乘以1.1增强效果较好
v=min(snpr.*snpo./(1+snpr),100); %v有时出现奇异值，所以将它控制在100之内
Gm=0.5.*sqrt(v.*pi)./snpo.*exp(-v./2).*((1+v).*besseli(0,v./2)+v.*besseli(1,v./2));%计算有音时的增益函数
gam=4*exp(v)./(1+snpr); %计算似然比,q设为0.2
G=gam./(1+gam).*Gm; %计算考虑有／无音概率的增益函数;增益函数还可考虑在无声时赋一个较小的值
M=Man.*G;
Mn=M.*exp(i.*En);
snpr1=M.^2./xd; %先验信噪比更新
x(1:256)=ifft(Mn(1:256));
x1=real(x).*hm; %进行IFFT后，有时会出现较小的复数值，所以加real函数.同时加合成窗
x0((n-1)*64+1:(n-1)*64+256)=x0((n-1)*64+1:(n-1)*64+256)+x1;
end
%****************归一化部分****************************************%
g0=hn.*hm; %由于在STFT分析合成时加了分析窗和合成窗,所以进行归一化,抵消加窗造成的影响
gl(1:64)=g0(1:64)+g0(65:128)+g0(129:192)+g0(193:256);
gain=repmat(gl,1,s);
x0=x0./gain;
wavwrite(real(x0),8000,8,'E:\学位论文写作\音源\结果\MMSE\m109.wav');%输出增强语音