%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    %
%            Fast LMS Algorithm                      %
%                                                    %
% Written By: Sundar Sankaran and A. A. (Louis) Beex %
%             DSP Research Laboratory                %
%             Dept. of Electrical and Comp. Engg     %
%             Virginia Tech                          %
%             Blacksburg VA 24061-0111               %
%                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

NoOfData = 8000 ;	% Set no of data points used for training
M = 32 ;		% Set the adaptive filter order

Mu = 0.01 ;		% Set the step-size constant
Gamma = 0.9 ;		% Forgetting factor
Delta = 0.01 ;		% R_est initialized to Delta*I

u = randn(NoOfData, 1) ;% Input assumed to be white
h = rand(M, 1) ;	% System picked randomly
d = filter(h, 1, u) ;	% Generate output (desired signal)

% Initialize fast-lms

W = zeros(2*M,1) ;
p = Delta*ones(2*M,1) ;

y = zeros(M, 1) ;
e = zeros(M, 1) ;

for k = 2 : floor(length(u)/M) - 1 ;

	U = fft(u((k-1)*M:(k+1)*M-1)) ;
	Y = ifft(U.*W) ;
	y(k*M:(k+1)*M-1) = real(Y(M+1:2*M)) ;
	
	e(k*M:(k+1)*M-1) = d(k*M:(k+1)*M-1) - y(k*M:(k+1)*M-1) ;
	E = fft([zeros(M,1); e(k*M:(k+1)*M-1)]) ;

	p = Gamma*p+(1-Gamma)*abs(U).^2 ;
	p_inv = 1./p ;

	PHI = ifft(p .* conj(U) .* E) ;
	phi = real(PHI(1 : M )) ;
	W = W + Mu / (2*M) * fft([phi; zeros(M,1)]) ;

end ;

% Plot results

figure ;
plot(20*log10(abs(e))) ;
title('Learning Curve') ;
xlabel('Iteration Number') ;
ylabel('Output Estimation Error in dB') ;