%谱减法程序
function [x,enhanced,snr]=pujian(filename,wsize,nlevel)
%if nargin<2
% fprintf('Usage: pujian(filename,noisyfile,outfile) \n\n');
% return;
%end
%读入语音文件
[speech,fs,nbits]=wavread(filename); %读入数据
%参数定义
winsize=wsize; %窗长
n=nlevel; %噪声水平
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
%除去Hamming窗引起的增益
for i=1:size
if hamwin(i)==0
enhanced(i)=0;
else
enhanced(i)=enhanced(i)/hamwin(i);
end
end
%计算增强前后信噪比
SNR1 = 10*log10(sum(speech.^2)/sum((speech'-x).^2)); %加噪语音信噪比
SNR2 = 10*log10(sum(speech.^2)/sum((speech'-enhanced).^2)); %增强语音信噪比
%disp('谱减法信噪比改善量（dB）：')
snr=SNR2-SNR1;
% wavwrite(x,fs,nbits,noisyfile); %输出带噪信号
% wavwrite(enhanced,fs,nbits,outfile); %输出增强语音
% %画波形
% figure('Name','谱减法','NumberTitle','off');
% 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]);
