spectrasubtractitera.m

语音处理中谱相减MATLAB程序

% 可选择 是否作为 函数 使用， 或 作为 单纯 迭代谱减法
%
%
function [sout,snd]=spectrasubtractiterative(y,num)

%y=rpure/max(rpure); % 做 函数使用是，用Y 参数，这里 单纯减噪声
%%% y=(c+rpure)/max(c+rpure);% 语音加  噪声

figure;
fs=4000;
[Pxx,w] = periodogram(y,[],'twosided',512,fs);
psdplot(Pxx,w,'Hz','','Sample PSD Plot');



sn=y(1:num);%   初始化 第一次参数, 函数参数
sout=y(1:num);

%  迭代次数  NM
for nm=1:2
    k=nm;
sn=sout;
szeros=zeros(1,num);sout=zeros(1,num);
j=sqrt(-1);
Snoise=zeros(1,256);
phase=zeros(1,256);
for n=1:256
   hanning(n)=(1/2)*(1-cos((2*pi)*(n-1)/255));% hanning wondow defination (see P231)
end
%以下检测无声帧，因为每次采样，前面都留出1秒左右（这里约4000个点）作为噪音的输入值，

noise1=sn(10:265).*hanning;
Noise1=abs(fft(noise1));
noise2=sn(300:555).*hanning;
Noise2=abs(fft(noise2));
noise3=sn(600:855).*hanning;
Noise3=abs(fft(noise3));

Noise=1.5*(Noise1+Noise2+Noise3)/3;  %噪音取均值，注意这里 的 过减因子1.5
%%%有时 过减因子 取得太小， 会造成信号较大畸变， 计算得的信噪比并不增加
%%%%

% i 长度239

  for i=1:fix(num/128-1)
      
      m=i; %一般开始一段无语音，不处理，最后一段的hanning 直接取 1
      
%    if m==1%这段要根据信号的长度count,重新编程
 %     snoise=sn((m-1)*128+1:(m-1)*128+256).*hanning; 
  %    snoise(1:128)=sn(1:128);
  %  else
  
    if m==fix(num/128-1)    %以下3条语句是对最后一个hanning窗的处理，即前128点重叠相加
    snoise=sn((m-1)*128+1:(m-1)*128+256).*hanning; % m=0，(1～256）x hanning ; m=1,(128～384）x hanning...,各段间有128个点的重叠
     snoise(129:256)=sn(num-127:num);   %后128点保持不变（不乘hanning窗，因为没有另128点与之累加）
    else
      snoise=sn((i-1)*128+1:(i-1)*128+256).*hanning; 
   end
    
   phase=angle(fft(snoise));
   Snoise=abs(fft(snoise));   
   
     for n=1:256     %公式参照语音信号处理P396  
        if(Snoise(n)^2-Noise(n)^2)<0
        Sout0(n)=0; % 因为频谱不可能是负数。
        else
        Sout0(n)=(Snoise(n)^2-Noise(n)^2)^0.5;
        end
        S0(n)=Sout0(n)*(cos(phase(n))+j*sin(phase(n)));
     end
    
   sout0=ifft(S0);
   szeros=zeros(1,num); %szeros先初始化
   szeros(((i-1)*128+1):((i-1)*128+256))=real(sout0); %在每段上是时域信号的实部
   sout=sout+szeros;  % 每一帧内时域信号的累加 ？？？？？？？？？？？？？？？？
   
  end
%snr=snratio(sout,y,240);
%figure;plot(sn);hold on;plot(sout,'m');

end
snd=sout/max(sout);
figure;plot(y);
hold on;
plot(sout,'m');