the spectral subtraction method.rar

close all; a1=0;a2=0;a3=0;a4=0;a5=0; nn1=zeros(5,1:26000); N=32;w=zeros(1,32);DD=zeros(1,32);var=0.99997;x=zeros(1,32); e=zeros(4,1:26000);Usf=0.0008; [n1,n2,n3,n4,n5]=dataread5('C:\MATLAB6p5\work\beamform\5ch8k40cm01.dat',80000); % 5ch8k1m4fu3cmcircle.dat',80000);21760);5ch8k1m4fu3cmcircle.dat'定义五个通道含噪声 sp=wavread('C:\MATLAB6p5\work\beamform\speech_dft.wav');%纯语音 for i=1:26000 speech(i)=sp(i*3); end %===初始化 %================================================= nn1(1,1:26000)=n1(1:26000)/50;nn1(2,1:26000)=n2(1:26000)/50; %自适应 nn1(3,1:26000)=n3(1:26000)/50;nn1(4,1:26000)=n4(1:26000)/50; nn1(5,1:26000)=n5(1:26000)/50; xx(1:26000)=nn1(5,1:26000)+speech; for aa=1:4 dd(1:26000)=nn1(aa,1:26000); for k=1:26000 if k<=32 x(1:k)=dd(k:-1:1); else x(1:32)=dd(k:-1:k-32+1); end y=w*x'; e(aa,k)=xx(k)-y; %e(aa,k)自适应输出 var=0.01*xx(k)^2+0.99*var; if k<=32 w(1:k)=w(1:k)+Usf*e(aa,k)*dd(k:-1:1)/var; else w(1:32)=w(1:32)+Usf*e(aa,k)*dd(k:-1:k-32+1)/var; end end end %spectral subtraction hanning=zeros(1,256);soutput=zeros(4,26000);szeros=zeros(1,26000); szeros=zeros(1,26000);Sout0=zeros(1,256); 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 for i=1:8 %======================================= noise1=e(1,i*256-255:i*256).*hanning; Noise1=abs(fft(noise1)); Noise(i,1:256)=Noise1; if i>1 Noise(i,1:256)=Noise(i,1:256)-Noise(1,1:256); end end S0=zeros(4,1:26000); as=4; for k=1:4 as for i=1:200 snoise=e(k,(i-1)*128+1:(i-1)*128+256).*hanning; %重叠目的：抵消窗效应 phase=angle(fft(snoise)); Snoise=abs(fft(snoise)); for m=1:8 for n=1:256 if(Snoise(n)^2-Noise(m,n)^2)<0 Sout0(n)=0; else Sout0(n)=(Snoise(n)^2-Noise(m,n)^2)^0.5; end S0(k,n)=Sout0(n)*(cos(phase(n))+j*sin(phase(n))); end Snoise=Sout0; end sout0=ifft(S0(k,1:256)); szeros(((i-1)*128+1):((i-1)*128+256))=real(sout0); soutput(k,1:26000)=soutput(k,1:26000)+szeros; szeros=zeros(1,26000); end end output=zeros(1,26000); for i=1:4 output=soutput(i,1:26000)+output; end subplot(4,1,1) plot(xx) subplot(4,1,2) plot(e(1,1:26000)) subplot(4,1,3) plot(soutput(1,1:26000)) subplot(4,1,4) plot(output (1:26000)) end end end end