%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This function is an implementation of the 
%   generalized NLMS algorithm (SPL-05) and the sign algorithm 
%   with adaptive learning rate for ICASSP-06 
%   submission
%
%    AUTHORS  :  Santhanam Balasubramaniam
%    DATE     :  10/05/05
%
%    USAGE [filco,out,MSE,eta] = gnlms(x,d,alpha,rho,ntap);
%
%    ifreq1,ifreq2 : estimated IF tracks
%    eta           : learning rate 
%    MSE           : Instantaneous Mean Square Error (dB)
%    ntap          : Number of taps in the filter
%    alpha         : Step size parameter (0 - 2)
%    rho           : Offset parameter learning rate (typ: 0-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [filco,out,MSE,eta] = gnlms(x,d,alpha,rho,ntap)
format long
x = x(:); d = d(:);
N = length(x); X = convmtx(x,ntap);
% Initialization
W(1,:) = zeros(1,ntap); 
[M,N] = size(X); beta(1) = 0.01; 
out = zeros(1,N); MSE(1) = 1;
for k = 2:1:M-ntap+1
    out(k) = W(k-1,:)*X(k,:).';
    e(k) = d(k) - out(k);
    cfactor(k) = e(k)*e(k-1)*X(k,:)*X(k-1,:).';
    num = rho*alpha*cfactor(k); 
    den = X(k-1,:)*X(k-1,:).' + beta(k-1);
    beta(k) = beta(k-1) - num/den; 
    clear num den cfactor out
    eta(k) = alpha / ((X(k,:) * X(k,:).') + beta(k));
    W(k,:) = W(k-1,:) + eta(k)*conj(X(k,:))*e(k);
    MSE(k) = e(k)^2; 
end
filco = W;
clear beta e  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%