% 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('chinese_2.wav',[1 20000]);
s=y(:,1);                                % speech source signal
% K=40000;
% k=1:K;
% s=0.1*sin(2*pi*k/K)';

%-----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^4));              % background noise(40db)
noise=wgn(1,L,a);

%-------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.5;                            % 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>=(L-M)/2                        % alter the impulse response in the receiving room to another one
        h=hh2;
    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)=norm(h-w)/norm(h);            % misalignment
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  draw picture                          
NN=length(e);
block=500;
E=[zeros(1,fix(block/2)),e,zeros(1,fix(block/2))];
D=[zeros(1,fix(block/2)),d,zeros(1,fix(block/2))];
for i=1:NN                               % 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)');   
end                                                   

%-----------plot-------------------------%
figure(1);
plot(10*log10(mis));                     % plot misalignment curve
ylabel('mis');xlabel('samples');
title('Misalignment of NLMS algorithm');
figure(2);
plot(10*log10(mse));                     % plot misalignment curve
ylabel('mse');xlabel('samples');
title('Mean Square Error of NLMS algorithm');

disp('norm-w = ');disp(norm(w));
disp('norm-h = ');disp(norm(h));

toc