% This is the Normalized Least-Mean-Square Algotithm(NLMS) of two_channel acoustic echo cancellation.
clear all
clc
tic
%-load the  impulse response-------------%
load g_1 Imp                             % one input signal in the receiving room
g1=Imp;
load g_2 Imp                             % another input signal in the receiving room
g2=Imp;
load h_1  Imp                            % one impulse response in the receiving room
h1=Imp;
load h_2  Imp                            % another impulse response in the receiving room 
h2=Imp;
load h_21  Imp                           % one impulse response in the receiving room
h21=Imp;
load h_22  Imp                           % another impulse response in the receiving room 
h22=Imp;

%------read the acoustic signal----------%
y=wavread('english.wav',[1 30000]);
s=y(:,1);                                % speech source signal

%%%%%%%%%----2.white gauss noise signal
% load white_noise wn;
% s=wn(1:30000);

%-----filter the input signal------------%
L1=length(s);                            % length of the input signal
N=length(g1);                            % order of the modeling filters
for i=1:L1-N+1
    s0=s(i+N-1:-1:i);
    s1(i)=s0'*g1';                       % get the signal filtered by impuse response of the transmition room
    s2(i)=s0'*g2';
end

%----preprocess the speech signal--------%
% beta=0.1;                                % factor of preprocessing
% s1=s1+beta*(s1+abs(s1));                 % preprocessing of nonlinear transform
% s2=s2+beta*(s2-abs(s2));

%--generate a background noise signal----%
L=length(s1);                            % length of the useful signal
p1=sum(abs(s1).*abs(s1))/L;
p2=sum(abs(s2).*abs(s2))/L;
a=10*log10((p1+p2)/(10^5));              % background noise(50db)
noise=wgn(1,L,a);
power_n=(p1+p2)/(10^5);
%-------set another parameters needed----%
M=length(h1);
hh1=[h1,h2]';                            % concatenate the two room inpulse response as one vector
hh2=[h21,h22]';
w=zeros(2*M,1);                          % weight vector of the adaptive filter
d=zeros(1,L-M);                          % desired signal
y=zeros(1,L-M);                          % adaptive filter output signal
e=zeros(1,L-M);                          % defference of the y and d signal
MSE=zeros(1,L-M);                        % mean square error 
mis=zeros(1,L-M);                        % misalignment 
stepsize=0.9;                            % adaptive stepsize 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NLMS  Algorithm
for i=1:L-M
    xx1=s1(i+M-1:-1:i);                  % one receiving singal of the receiving room                       
    xx2=s2(i+M-1:-1:i);                  % another receiving singal of the receiving room
    xx=[xx1 xx2]';                       % concatenate the two room inpulse response as one vector
    h=hh1;
%     if i>=20000                          % alter the impulse response in the receiving room to another one
%         h=hh2;
%     end
%     if i==19999
%         a=norm(w);
%     end
    d(i)=h'*xx;%+noise(i);                 % desired signal 
    y=w'*xx;                             % filter signal
    e(i)=d(i)-y;                         % error signal
    q=xx'*xx;                            % norm of input signal
    w=w+stepsize*e(i)*xx/q;              % iterate the coefficient of modeling filters
    mis(i)=10*log10(norm(h-w)/norm(h));  % misalignment
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  draw picture                          
NN=length(e);
block=500;
for i=1:NN-block+1                       % mean square error and smoothed with 300 data
    MSE(i)=e(i:i+block-1)*e(i:i+block-1)'/(d(i:i+block-1)*d(i:i+block-1)');   
    mse(i)=10*log10(MSE(i));
    MSE2(i)=[(d(i:i+block-1)*d(i:i+block-1)')-block*power_n]/[e(i:i+block-1)*e(i:i+block-1)'-block*power_n];   
    mse2(i)=10*log10(MSE2(i));
end                                                    

%-----------plot-------------------------%
figure(1);
plot(mis);                     % plot misalignment curve
ylabel('mis');xlabel('samples');
title('Misalignment of NLMS algorithm');
figure(2);
plot(mse);                     % plot misalignment curve
ylabel('mse');xlabel('samples');   %axis([0 L1 -400 0]);
title('Mean Square Error of NLMS algorithm');
% figure(3)
% plot(mse2);
title('ERLE');
%disp('norm-a=');disp(a);
disp('norm-w = ');disp(norm(w));
disp('norm-h = ');disp(norm(h));

save NLMS_mse mse;
save NLMS_mis mis;

toc