pujianfa.zip_ pujianfa.rar_减谱_噪声估计_改进谱减

[s1,fs,nbits]=wavread('weizhi.wav'); s=(s1(:,1)+s1(:,2))/2; figure(2) subplot(321);plot(s);title('原始语音'); Wp=[0.025*pi,0.85*pi]; Ws=[0.010*pi,0.95*pi]; Ap=10; As=60; p=0.0001;%p为白噪声信号的方差 alpha=0.1;%预加重系数 beta=0.95;%去加重系数 M=128; %产生白噪声 n_rand=randn(1,length(s)); n_rand_mean=mean(n_rand); n_rand_var=var(n_rand); n=n_rand-n_rand_mean; n=(n*sqrt(p/n_rand_var))'; x=s+n; subplot(322);plot(x);title('带噪语音'); SNR1=10*log10(var(s)/var(x-s))%计算信噪比 aa1=length(x); [N,Wn]=buttord(Wp/pi,Ws/pi,Ap,As);%产生带通滤波器 [b,a]=butter(N,Wn,'bandpass'); [H,w]=freqz(b,a,length(x)); mag=abs(H); db=20*log10(mag/max(mag)); pha=angle(H); grd=grpdelay(b,a,w);%群延迟 %带噪信号通过带通滤波器 x_bap=filter(b,a,x); subplot(323);plot(x_bap);title('经过预处理后的信号'); %进行预加重 y(1)=x_bap(1); for j=2:length(x_bap) y(j)=x_bap(j)-alpha*x_bap(j-1); end subplot(324);plot(y);title('预加重的信号');axis([0,6*10^4,-0.5,0.5]); %分帧并加窗处理 zhenshu=floor(length(y)/(M/2))-1; ham=hamming(M); for j=1:zhenshu for k=1:M j1=(M/2)*(j-1); y_frame(j,k)=y(j1+k); y_ham(j,k)=y_frame(j,k).*ham(k); end end [a,b]=size(y_ham); %对信号进行傅立叶变换并求幅度谱和相位谱 y_fft=fft(y_ham',128); [c,d]=size(y_fft); y_fft_angle=angle(y_fft); y_fft_mag=abs(y_fft); y_w=sum(y_fft_mag); [e,f]=size(y_w); %噪声估计 %根据前20帧计算平均功率 ave=0; for j=1:20 ave=ave+y_w(j); end ave=ave/20; EMAX=max(y_w(1:20)); EMIN=min(y_w(1:20)); T=min(0.03*(EMAX-EMIN)+ave,4*ave); for j=1:length(y_w) if y_w(j)<T noise_w(j)=y_w(j); else noise_w(j)=T; end end %中值平滑 aa=length(noise_w);%3686 for j=2:(aa-1) array=[noise_w(j-1),noise_w(j),noise_w(j+1)]; array=sort(array); noise_w_mid(j)=array(2); end noise_w_mid(1)=ave; noise_w_mid(length(noise_w))=ave; cc=length(noise_w_mid);%3686 %线性平滑 noise_w_line(1)=ave; noise_w_line(2)=ave; bb=length(noise_w_mid);%3686 for j=3:(bb-2) noise_w_line(j)=1/9*noise_w_mid(j-2)+2/9*noise_w_mid(j-1)+3/9*noise_w_mid(j)+2/9*noise_w_mid(j+1)+3/9*noise_w_mid(j+2); end noise_w_line(length(noise_w_mid)-1)=ave; noise_w_line(length(noise_w_mid))=ave; dd=length(noise_w_line);%3686 %noise_delta为残差噪声，noise_w_line为新噪声，noise_w为原噪声 noise_delta=noise_w-noise_w_line; %用中值平滑平滑残差噪声 ee=length(noise_delta); for j=2:(ee-1) array=[noise_delta(j-1),noise_delta(j),noise_delta(j+1)]; array=sort(array); noise_delta_mid(j)=array(2); end noise_delta_mid(1)=0; noise_delta_mid(length(noise_delta))=0; %用线性平滑更新noise_delta_mid noise_delta_line(1)=0; noise_delta_line(2)=0; ff=length(noise_delta_mid);%3686 for j=3:(ff-2) noise_delta_line(j)=1/9*noise_delta_mid(j-2)+2/9*noise_delta_mid(j-1)+3/9*noise_delta_mid(j)+2/9*noise_delta_mid(j+1)+3/9*noise_delta_mid(j+2); end noise_delta_line(length(noise_delta_mid)-1)=0; noise_delta_line(length(noise_delta_mid))=0; %噪声功率谱noise_w_line和残差噪声功率谱noise_delta_line之和，得到noise_w_com noise_w_com=noise_w_line+noise_delta_line; gg=length(noise_w_com); %谱减 y_w_dec=y_w-noise_w_com; for j=1:length(y_w_dec) if y_w_dec(j)<0 y_w_dec(j)=0; else y_w_dec(j)=y_w_dec(j); end end hh=length(y_w_dec); %利用相位谱 for j=1:M for k=1:floor(length(y)/(M/2))-1 y_percent(j,k)=abs(y_fft(j,k))/y_w(k); end end ii=length(y_percent);%到这还运行的挺快 for k=1:length(y_w_dec) for j=1:M y_fft_new(j,k)=y_w_dec(k)*y_percent(j,k)*exp(i*y_fft_angle(j,k)); end end [aaa,bbb]=size(y_fft_new); %做反傅立叶变换 y_new_frame=real(ifft(y_fft_new)); jj=length(y_new_frame) %根据信号帧y_new_frame恢复原语音信号y_new h=1; for j=1:M/2 y_new(h)=y_new_frame(j,1); h=h+1; end%到这时速度有所下降但还行 for j=2:floor(length(y)/(M/2))-1 for k=1:M/2 y_new(h)=y_new_frame(k,j)+y_new_frame(k+M/2,j-1); h=h+1; end end kk=length(y_new);%上面这个循环运行的就有点慢了 %做去加重 for j=1:length(y_new)-1 x_new(j)=y_new(j)+beta*y_new(j+1); end x_new(length(y_new))=y_new(length(y_new)); x_new=x_new(1:length(x)); kkk=length(x_new) subplot(325);plot(x_new);title('增强语音');axis([0,6*10^4,-0.5,0.5]); %输出和播放 wavwrite(s,'yuanshiyuyin'); wavplay(s,fs); wavwrite(x,'daizaoyuyin'); wavwrite(x_new,'zengqiangyuyin'); wavplay(x,fs); wavplay(x_new,fs);