% 

% Matlab codes for the following two papers

 % 

%F. Albu, J. Kadlec, N. Coleman, A. Fagan, "The Gauss-Seidel Fast Affine Projection Algorithm", IEEE Workshop SIPS 2002, pp. 109-114, San Diego, U.S.A, October 2002

% 

%F. Albu, A. Fagan, The Gauss-Seidel pseudo affine projection algorithm and its application for echo cancellation, Conference Record of the Thirty-Seventh Asilomar Conference on Signals, Systems and %Computers, 2003, Vol. 2,  9-12 Nov. 2003, pp. 1303-1306

%

%%% Felix Albu

%%% Valahia University of Targoviste, Romania

%%% Email: felix.albu@valahia.ro

%%% http://falbu.50webs.com

%

if 1 % 1 for executing GSPAP

echo off

clear

randn('seed', 0) ;   rand('seed', 0) ;

L=500;  % length of the adaptive filter

N=10;       % size of the affine projection

mu=1;       % FAP step size

update_parameter=10; % updating parameter

load fisierecho % speech signal+echo file

sizex=fisierecho(:,1); % speech signal

ss=fisierecho(:,2); % noisy echo (signal-to-near-end-noise ratio = 30 dB

maxs=max(sizex);

Delta= maxs*maxs*L/1000;

nrsamp=length(sizex);

eout=zeros(nrsamp,1);   %the fap error over the length of the experiment

% initializations

V_matrix=zeros(N,1);    %newest excitation vector

old_V=zeros(N,1);  %oldest excitation vector

extended_V_matrix=zeros(L+1,1); % the length of this vector has L+1 instead L+N of GSFAP !

W_aux_matrix=zeros(L,1);    % the fap-adaptive filter real coefficients !

p_matrix=[1/Delta; zeros(N-1,1)];

b_vector=[1; zeros(N-1,1)];

rcorr_matrix = Delta * eye ( N, N );

rcorr_buf=[Delta;zeros(N-1,1)]; % it has a N dimension !

inv_vector=[1/Delta;zeros(N-1,1)];

% all the vectors from above are also used by GSFAP

% here is the different initialisation part

% epsilon_buf and neta_buf from GSFAP are not needed

extended_U_matrix=zeros(L+1,1); % this new vector has a L+1 length

den1=Delta;

for iter=1:nrsamp

   if rem(iter, 1000)==0, iter, end;

% shift state vectors

V_matrix(2:N)=V_matrix(1:N-1);

V_matrix(1)=sizex(iter);

old_V(2:N)=old_V(1:N-1);   

extended_V_matrix(2:L+1)=extended_V_matrix(1:L);

extended_V_matrix(1)=V_matrix(1);

old_V(1)=extended_V_matrix(L+1);

% update fap paramaters

rcorr_buf=rcorr_buf+V_matrix(1)*V_matrix(1:N)-old_V(1)*old_V(1:N);

rcorr_matrix(2:N,2:N)=rcorr_matrix(1:N-1,1:N-1);

rcorr_matrix(:,1)=rcorr_buf;

rcorr_matrix(1,:)=rcorr_buf';

inv_vector(2:N)=inv_vector(1:N-1);

inv_vector(1)=1/rcorr_buf(1);

% GS section

if rem(iter,update_parameter)==0

   for i=1:N

      sum1=0;

      for j=1:N

         if j~=i

               sum1=sum1-rcorr_matrix(i,j)*p_matrix(j);

          end

      end

      p_matrix(i)=(b_vector(i)+sum1)*inv_vector(i);  

   end

end

err_buf=ss(iter)-extended_V_matrix(1:L)'*W_aux_matrix;

err_buf1=ss(iter)-extended_V_matrix(1:L)'*W_aux_matrix1;

%[iter err_buf1-err_buf]

% all equations from above are identical in GSFAP

% the different part is following

extended_U_matrix(2:L+1)=extended_U_matrix(1:L);

extended_U_matrix(1)=V_matrix'*p_matrix/p_matrix(1);

den1=den1+extended_V_matrix(1)*extended_U_matrix(1)-extended_V_matrix(L+1)*extended_U_matrix(L+1); % simpler way to compute U'(n)*X(n)+Delta

%den1=den1+extended_U_matrix(1)*extended_U_matrix(1)-extended_U_matrix(L+1)*extended_U_matrix(L+1); % simpler way to compute U'(n)*U(n)+Delta -- a possible alternative

W_aux_matrix=W_aux_matrix+mu*err_buf*extended_U_matrix(1:L)/den1;

eout(iter)=err_buf;     % eout is the FAP output

end

figure, plot(ss), hold on, plot(eout,'r'), legend('noisy echo','residual echo')

end % if

if 0 % 1 for executing GSFAP

%======================================================================================

echo off

clear

randn('seed', 0) ;   rand('seed', 0) ;

L=500;  % length of the adaptive filter

N=10;       % size of the affine projection

mu=1;       % FAP step size

update_parameter=10; % updating parameter

load fisierecho % speech signal+echo file

sizex=fisierecho(:,1); % speech signal

ss=fisierecho(:,2); % noisy echo (signal-to-near-end-noise ratio = 30 dB

maxs=max(sizex);

Delta= maxs*maxs*L/1000;

nrsamp=length(sizex);

% initializations

V_matrix=zeros(N,1);    %newest excitation vector

old_V=zeros(N,1);  %oldest excitation vector

extended_V_matrix=zeros(L+N,1); % this contains N length L excitation vectors

epsilon_buf=zeros(N,1); % the "epsilon" vector in fap

neta_buf=zeros(N,1);    % the E-vector in fap

W_aux_matrix=zeros(L,1);    % the fap-adaptive filter auxiliary coefficients

eout=zeros(nrsamp,1);   %the fap error over the length of the experiment

p_matrix=[1/Delta; zeros(N-1,1)];

b_vector=[1; zeros(N-1,1)];

rcorr_matrix = Delta * eye ( N, N );

rcorr_buf=[Delta;zeros(N-1,1)];

inv_vector=[1/Delta;zeros(N-1,1)];

for iter=1:nrsamp

   if rem(iter, 1000)==0, iter, end;

% shift state vectors

V_matrix(2:N)=V_matrix(1:N-1); 

V_matrix(1)=sizex(iter);   

old_V(2:N)=old_V(1:N-1);   

extended_V_matrix(2:L+N)=extended_V_matrix(1:L+N-1);

extended_V_matrix(1)=V_matrix(1);      

old_V(1)=extended_V_matrix(L+1);

% update fap paramaters

rcorr_buf=rcorr_buf+V_matrix(1)*V_matrix(1:N)-old_V(1)*old_V(1:N);

rcorr_matrix(2:N,2:N)=rcorr_matrix(1:N-1,1:N-1);

rcorr_matrix(:,1)=rcorr_buf;

rcorr_matrix(1,:)=rcorr_buf';

inv_vector(2:N)=inv_vector(1:N-1);inv_vector(1)=1/rcorr_buf(1);

% GS section

if rem(iter,update_parameter)==0

   for i=1:N

      sum1=0;

      for j=1:N

         if j~=i

               sum1=sum1-rcorr_matrix(i,j)*p_matrix(j);

          end

      end

      p_matrix(i)=(b_vector(i)+sum1)*inv_vector(i);  

   end

end

W_aux_matrix=W_aux_matrix+mu*neta_buf(N)*extended_V_matrix(N+1:N+L);

err_buf=ss(iter)-extended_V_matrix(1:L)'*W_aux_matrix;

err_buf=err_buf-mu*rcorr_buf(2:N)'*neta_buf(1:N-1);

epsilon_buf=err_buf*p_matrix;

neta_buf=[0; neta_buf(1:N-1)]+epsilon_buf;

eout(iter)=err_buf;     % eout is the FAP output

end

figure, plot(ss), hold on, plot(eout,'r'), legend('noisy echo','residual echo')

end % if