基本谱减法

clc
clear all
%-----------------------读入语音文件-----------------------------
[speech,fs,nbits]=wavread('yuyinshiyan2.wav');%读入数据
%----------------------参数定义---------------------------------
speech = speech(:,1);%双通道数据改为单通道数据
winsize=256;%窗长
n=0.04;%噪声水平
size=length(speech);%语音长度
numofwin=floor(size/winsize);%帧数
ham=hamming(winsize)';%产生汉明窗
hamwin=zeros(1,size);%定义汉明窗的长度
enhanced=zeros(1,size);%定义语音增强的长度
x=speech'+ n*randn(1,size);%产生带噪信号
noisy=n*randn(1,winsize);%噪声估计
N=fft(noisy);%对噪声进行傅里叶变换
nmag=abs(N);%噪声功率谱
%-------------------------------分帧----------------------------
for q=1:2*numofwin-1%对带噪语音帧间重叠一半取值
frame=x(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2);
hamwin(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)=...
hamwin(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)+ham;%加窗
y=fft(frame.*ham);%对带噪语音进行傅里叶变换
mag=abs(y);%带噪语音功率谱
phase=angle(y);%带噪语音相位
%-----------------------------幅度谱减---------------
for i=1:winsize
if mag(i)-nmag(i)>0
clean(i)=mag(i)-nmag(i);
else
clean(i)=0;
end
end
%-------------------在频域中重新合成语音------------------
spectral=clean.*exp(j*phase);
enhanced(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)=...
enhanced(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)+real(ifft(spectral));
end
%---------------------除去汉明窗引起的增益-----------------------
for i=1:size
if hamwin(i)==0
enhanced(i)=0;
else
enhanced(i)=enhanced(i)/hamwin(i);
end
end
%------------------------------计算增强前后信噪比---------------------------
SNR1=10*log10(var(speech')/var(noisy));%加噪语音的信噪比
SNR2=10*log10(var(speech')/var(enhanced-speech'));%增强语音信噪比
wavwrite(x,fs,nbits,'noisy.wav');%输出带噪信号
wavwrite(enhanced,fs,nbits,'enhanced.wav');%输出增强语音信号
%-------------------------画波形-----------------------------------
figure(1);
subplot(3,1,1);plot(speech');
title('原始语音波形');
xlabel('样点数');ylabel('幅度');
axis([0 2.5*10^4 -0.3 0.3]);
subplot(3,1,2);plot(x);
title('加噪语音波形');
xlabel('样点数');ylabel('幅度');
axis([0 2.5*10^4 -0.3 0.3]);
subplot(3,1,3);plot(enhanced);
title('语音增强波形');
xlabel('样点数');ylabel('幅度');axis([0 2.5*10^4 -0.3 0.3]);