%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the Block Least-Mean-Square Algotithm(BLMS) of two_channel acoustic echo cancellation.
% the impulse response of the receiving room and the transimission room are   order.
% the order of the modeling filter is 28 for each channel.
% when length of input=40000,the algorithm needs  1.5470s to calculate.
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  initializing parameter variable
%load source_white_noise           % white noise  input signal 
%source=source(1:42000);
load sound                       % speech input signal
source=sound(1:42000);
load h1                           % one impulse response in the receiving room
load h2                           % another impulse response in the receiving room 
load g1                           % one impulse response in the transimission room
load g2                           % another impulse response in the transmission room 
load gg1                          % one new impulse response in the transmission room
load gg2                          % another new impulse response in the transmission room 
length_source=length(source);     % length of source signal
N=length(h1);                     % order of the modeling filters 
h_filter=[h1 h2];                 % concatenate the two room inpulse response as one vector             
for i=N:length_source 
    x=source(i:-1:i-N+1); 
    x1(i)=x*g1';                      % signal filtered by g1
    x2(i)=x*g2';                      % signal filtered by g2
    if i>=26624                       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         x1(i)=x*gg1';                % the impulse response in the transmission room changed when iterated 35000 times
         x2(i)=x*gg2';                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
end
h1=h1';h2=h2';
zeros_m=zeros(N,1);
h1=[h1;zeros_m];h2=[h2;zeros_m];  
W1=zeros(1,2*N)';   W2=zeros(1,2*N)';% initializing  the modeling filters  
stepsize=0.9;                        % adaptive stepsize 
x_old1=x1(1:N);    x_old2=x2(1:N);
x1=x1(N+1:length_source);x2=x2(N+1:length_source);
length_source=length(x1);             % length of source signal
%P1=0.00000001*ones(2*N,1);P2=0.00000001*ones(2*N,1);
P=0.00000001*ones(2*N,1);
gama=0.01;
H1=fft(h1);H2=fft(h2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BLMS  Algorithm
tic
for i=1:length_source/N
    x_new1=x1(N*i-N+1:N*i);      x_new2=x2(N*i-N+1:N*i); 
    xx1=[x_old1 x_new1]';        xx2=[x_old2 x_new2]';
    x_old1=x_new1;               x_old2=x_new2;
    U1=fft(xx1);                 U2=fft(xx2);
    Q_d1=U1.*H1;                 Q_d2=U2.*H2;
    Q_y1=U1.*W1;                 Q_y2=U2.*W2;
    y1=ifft(Q_y1);               y2=ifft(Q_y2);
    d1=ifft(Q_d1);               d2=ifft(Q_d2);
    d1=d1(N+1:2*N);              d2=d2(N+1:2*N);
    y1=y1(N+1:2*N);              y2=y2(N+1:2*N);
    d(N*i-N+1:N*i,1)=d1+d2;
    e(N*i-N+1:N*i,1)=d1+d2-y1-y2;
    temp_e=[zeros_m;e(N*i-N+1:N*i)];
    E=fft(temp_e);
    P=gama*P+(1-gama)*(U1.*conj(U1)+U2.*conj(U2));
    D=1./P;
    R1=D.*conj(U1).*E;
    Q1=ifft(R1);
    Q1=Q1(1:N);
    temp_Q1=[Q1;zeros_m];
    p1=fft(temp_Q1);
    R2=D.*conj(U2).*E;
    Q2=ifft(R2);
    Q2=Q2(1:N);
    temp_Q2=[Q2;zeros_m];
    p2=fft(temp_Q2);
    W1=W1+stepsize*p1; 
    W2=W2+stepsize*p2; 
    w1=ifft(W1);                 w2=ifft(W2);
    w1=w1(1:N);                  w2=w2(1:N);
    w=[w1;w2]';
    mis(i)=norm(h_filter-w)/norm(h_filter);     % misalignment
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  draw picture
e=e.';
d=d.';
NN=length(e);
for i=1:NN-599                           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    mse(i)=e(i:i+599)*e(i:i+599)'/(d(i:i+599)*d(i:i+599)');    % mean square error and smoothed with 100 data
end                                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
figure(1)    
plot(w1,w2)                       % plot misalignment curve
ylabel('w')
figure(2)
plot(e)                       
ylabel('e')