echocanceler3_4.m

用带leaky因子的GEIGEl算法进行语音状态检测的回声消除程序及其性能测试曲线

clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%产生作为输出缓噪声的伪噪声序列%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% temp=(randn(1,8192)>0);
% PN=zeros(1,16384);
% PN(1:8192)=temp;
% PN(16384:-1:8193)=-temp;
% w_k=zeros(1,16384);
% w_k(1:3715)=1;
% w_k(end-3715:end)=1;
% FFT_value=w_k.*exp(j*PN*pi);
% FFT_temp(1)=0;
% FFT_temp(2:16385)=FFT_value;
% IFFT_temp=ifft(FFT_temp,16385);
% IFFT_value=IFFT_temp(2:16385);
% FINAL_value=fftshift(real(IFFT_value));
% FINAL_value=resample(FINAL_value,2,11);
% PN_value=FINAL_value;
% PN_num=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%产生一个回声路径的冲激响应%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 240; 
t = [0 : N-1]'; 
tau = N./4; 
envelope = exp(-(t./tau));
w_k = randn(N,1).*envelope; 
w_k = [w_k(end-19 : end);w_k(1 : end-20)]; 
w_k = w_k./sqrt(N);
% figure; 
% plot(w_k);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%输入两段语音信号作为远端和近端信号%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% array1=wavread('quiet',100000);%近端语音
% array2=wavread('wang2',100000);%远端语音

% array1=wavread('wang1',300000);
% array2=wavread('wang2',300000);
% array1(1:100000)=wavread('quiet',100000);
% array2(end-69999:end)=wavread('quiet',70000);

% array1=wavread('quiet',300000);
array2=wavread('wang2',300000);

% array1=array1';
array2=array2';
%wavplay(array1,22050);
%wavplay(array2,22050);
M=length(array2); %序列总长度

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%通过远端语音信号与回声路径冲激响应卷积生成回声信号%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=length(w_k);%滤波器阶数
pecho=zeros(1,N);
pecho=conv(array2,w_k); 
pecho=pecho(1:M);%模拟的回声信号
array1= 0.05*mean(abs(pecho))*randn(1,M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%读入N个数据，设置初值%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rin=array2; %远端信号
% Sinn=pecho;
Sinn=pecho+array1;%近端信号——近端说话人的语音信号与回声叠加构成
% wavplay(pecho,8000);
% wavplay(Sinn,8000);

x=Rin(1:N);
w=zeros(1,N);
a=0.0001;
u=0.5;

energ_window=1/32;         %短窗功率估计
near_energ=0;
far_energ(1)=0;
for i=2:N-1
    far_energ(i)=(1-energ_window)*abs(far_energ(i-1))+energ_window*abs(x(i-1));
end;

Dtemp=x(1:N-1)*x(1:N-1)';%滤波器输入信号的能量
sout=zeros(1,M);
ff=0;
dd=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%NLMS算法的自适应滤波%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=N:M; 
    x(N)=Rin(i);
    D=Dtemp+x(N).^2;
    far_energ(N)=(1-energ_window)*abs(far_energ(N-1))+energ_window*abs(x(N-1));
    near_energ=(1-energ_window)*abs(near_energ)+energ_window*abs(Sinn(i-1));
    y=dot(w,x);
    en=Sinn(i)-y;  
    if near_energ>1/2*max(far_energ);
        sout(i)=en;
        dd=dd+1;
%       sout(i)=nlp_new(en); ;%远端模式，启动NLP 
     else;
         w=w+u*en*x/(a+D);    
         sout(i)=en;
%         sout(i)=nlp_new(en); ;%远端模式，启动NLP 
        ff=ff+1;
    end;

    
    %%%%%%%%%更新滤波器输入信号的能量值%%%%%%%%%%
    Dtemp=D-x(1).^2;
    x(1:N-1)=x(2:N);
    far_energ(1:N-1)=far_energ(2:N);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
% figure;
% plot(array1);
% hold on;
% plot(pecho,'r');
% figure;
% plot(Sinn);
figure;
plot(sout);
figure; 
plot(Sinn); hold on; plot(sout,'r');
legend('near end signal','error signal');

wavwrite(Sinn,8000,16,'Sinn');
wavwrite(Rin,8000,16,'Rin');
wavwrite(sout,8000,16,'sout');
%wavplay(sout,22050);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%性能测试%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_Sinn=abs(Sinn).^2;
P_en=abs(sout).^2;
Ave_num=5600;
P_num=fix(M/Ave_num);
ERLE_temp=P_Sinn./P_en;
for j=1:P_num;
    Mark=(j-1)*Ave_num;
    ERLE_win=ERLE_temp(Mark+1:Mark+Ave_num);
    ERLE(j)=sum(ERLE_win)/Ave_num;
end;
ERLE=10*log10(ERLE);
figure;
plot(ERLE);