function [err,y,h,k]=gall(x,d,M,beta,alpha,epsi,k0,h0)
%
% [err,y,h,k]=gall(x,d,M,beta,alpha,epsi,k,h)
%
% Gradient adaptive Laguerre Lattice Filter. Implements an adaptive
% filter with the option to include an adaptable pole which can be
% useful to model system with long impulse response. Step size is
% optional. The filter operates in two modes: 1) Learning mode
% and 2)Output (Freeze mode).
%
% ***Mode 1:The arguments for the Learning mode are:
%
% x -measurement data
% d -desired response or reference signal (same size as x)
% M -filter order
% beta -Forgetting factor (0<=beta<=1). beta =1 remembers all the data.
% alpha -Scalar value of the magnitude of the pole in the filter (0<= alpha <=1).
% Alpha of 0 is equivalent to the standard classical GAL filter.
%epsi -small positive constant that initializes the algorithm (epsi <<1)
%k0 -optional 1xM+1 initial reflection coefficients vector
%h0 -optional 1xM+2 initial FIR ladder coefficient vector
%
% Function returns:
%
% err -learning error curve for the algorithm
% y -prediction curve (dhat)
% h -FIR coefficients of the ladder
% k -Lattice coefficients
%
%
% ***Mode 2:The arguments for the Ouput (Freeze) mode are:
%
% x -measurement data
% d -desired response or reference signal (same size as x)
% M -empty []
% beta -[]
% alpha -Scalar value of the magnitude of the pole in the filter (0<= alpha <=1).
% Alpha of 0 is equivalent to the standard classical GAL filter.
% k -Reflection coefficients of the trained lattice.
% h -Coefficients of the second stage FIR filter.
% Function returns:
%
% err -Error between filter's output and reference signal
% y -Filter's output
% h -[]
% k -[]
%
% Reference: Fezjo & Lev-Ari (1997) IEEE Trans. SP
%
% Copyright (C) 2010 Ikaro Silva
%
% This library is free software; you can redistribute it and/or modify it under
% the terms of the GNU Library General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option) any
% later version.
%
% This library is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE. See the GNU Library General Public License for more
% details.
%
% You should have received a copy of the GNU Library General Public License along
% with this library; if not, write to the Free Software Foundation, Inc., 59
% Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% You may contact the author by e-mail (ikaro@ieee.org).
%
% %%******Example: Robuts Echo cancellation
% %Algorithm performs standard echo cancelation at t<0 and than
% %a sudden cross-talk inteference is added at t=0.
%
% N=2^12;
% h=[1.0000 -1.8976 1.1870 -0.0642 -0.2718 0.0858];
% h2=[1.0000 -2.0134 1.4571 -0.2602 -0.2377 0.1015];
% fs=20;
% epsi=10^-2;
% w=10;
% t1=-50;
% t=linspace(t1,t1+N/fs,N);
% [trash,N0]=min(abs(t));
% S=zeros(1,N);
%
% for i=1:20
% %%%SIGNAL GENERATION
% %No cross-talk Phase
% s=sawtooth(t,0.5);
% n=randn(1,N);
% d=s+2*n;
% r=filter(1,h,n*2);
%
% %Add Cross-Talk Component to reference signal
% xtalk=r;
% q=10;
% talk=filter(ones(q,1)./q,1,s(N0:end));
% xtalk(N0:end)=xtalk(N0:end) + talk*std(xtalk)./std(talk);
% dtalk=d;
%
% %%%TECHNIQUE%%%%
% %Cross Talk Resistance - using Two GAL Filters
% x1=xtalk(1:N0-1);
% d1=dtalk(1:N0-1);
% %Initial filter- Estimates Noise TF and Freeze
% [err,n_est3,w3,k_xtr]=gall(x1,d1,w,1,0,10^-3,[],[]);
% [s_est3,n_est3,trash,trash]=gall(xtalk,dtalk,[],[],0,[],k_xtr,w3);
% s_est3=s_est3(1:end-1,end);
%
%
% %Estimate XTalk Interference
% [n_est,xtalk_est,w3,k_xtr]=gall(s_est3(N0:end-1,end),n_est3(N0:end),w,0.9975,0.001,10^-3,[],[]);
% s_est=s_est3';
%
% %s_est4(N0:end)=[zeros(1,w-1) dtalk(N0:end-w+1)]' - n_est4(1:end,end);
% s_est(N0:end)=dtalk(N0:end)' - n_est(1:end,end);
% s_est(N0:end)=[s_est(N0+w:end) zeros(1,w)];
% S=s_est./i + (i-1).*S./i;
% end
%
% %Plot
% figure
% plot(t,S,'r');hold on ;plot(t,s)
% legend('Gall XTR', 'Desired Response')
N=length(x);
if(~isempty(M))
%Multistage lattice predictor
f=zeros(1,M+1);
b=zeros(1,M+1);
bt=zeros(1,M+1);
Q=zeros(1,M+1)+epsi;
delta=zeros(1,M+1);
y=zeros(N,1);
D=zeros(1,M+2);
err=zeros(N+1,M+2);
b_old=b;
bt_old=bt;
if(isempty(k0))
k=zeros(1,M+1);
else
k=k0;
end
if(isempty(h0))
h=zeros(1,M+2);
else
h=h0;
end
for n=1:N
f(1)=alpha*b_old(1) + sqrt(1-alpha^2)*x(n);
b(1)=f(1);
Q(1) = beta.*Q(1) + b(1).^2;
err(n,1)=d(n);
bt_old=bt;
b_old=b;
for m=2:M+1
%Lattice Section
bt(m-1)= b_old(m-1) + alpha.*(bt_old(m-1) - b(m-1));
Q(m-1)=beta.*Q(m-1) + (f(m-1).^2 + bt(m-1).^2)./2;
delta(m)=beta.*delta(m) + f(m-1).*conj(bt(m-1));
k(m) = delta(m)./Q(m-1);
f(m)= f(m-1) - k(m)*bt(m-1);
b(m)= bt(m-1) - conj(k(m))*f(m-1);
end
Q(end)=beta.*Q(end) + (f(end).^2 + bt(end).^2)./2;
for m=1:M+1
%Ladder Section
D(m+1)= beta.*D(m+1) + err(n,m).*conj(b(m));
h(m+1)= D(m+1)./Q(m);
err(n,m+1)=err(n,m) - h(m+1).*conj(b(m));
end
y(n)= h(2:end)*conj(b)';
end
else
%Multistage lattice predictor
k=k0;
h=h0;
M=length(k)-1;
f=zeros(1,M+1);
b=zeros(1,M+1);
bt=zeros(1,M+1);
y=zeros(N,1);
err=zeros(N+1,M+2);
b_old=b;
bt_old=bt;
for n=1:N
f(1)=alpha*b_old(1) + sqrt(1-alpha^2)*x(n);
b(1)=f(1);
err(n,1)=d(n);
bt_old=bt;
b_old=b;
for m=2:M+1
%Lattice Section
bt(m-1)= b_old(m-1) + alpha.*(bt_old(m-1) - b(m-1));
f(m)= f(m-1) - k(m)*bt(m-1);
b(m)= bt(m-1) - conj(k(m))*f(m-1);
end
for m=1:M+1
%Ladder Section
err(n,m+1)=err(n,m) - h(m+1).*conj(b(m));
end
y(n)= h(2:end)*conj(b)';
end
h=[];
k=[];
end
%err=err(:,end);