function [e,h] = nlms(h,x,y,mu_tilde,start_time,stop_time)
%NLMS [e,h] = nlms(h,x,y,mu_tilde,start_time,stop_time) adjusts the filter h 
% such that the output of h when driven by x is close (in the least-squares 
% sense) to y.  The resulting error signal and final value of h are returned.
% The normalized lms has a time-varying step-size inversely proportional 
% (proportionality constant, mu_tilde) to the input power.  The first four 
% arguments are required, the last two must both be specified or not.
%   
%  If start_time is specified, the process begins adaptation at the 
%  start_time-th sample else adaptation begins immediately.  Delaying 
%  adaptation can be used to allow startup transients to settle before 
%  adaptation begins.
%
%  If stop_time is specified, the process stops adaptation at stop_time-th 
%  sample else adaptation continues through the last sample.  stop_time must 
%  be > N.  For adaptation to continue through the last sample, set stop_time 
%  = len(x)+1.  Stopping adaptation can be used to study the error without 
%  continued adaptation.    
%
%  Note: This function will work with complex-valued h,x, and y.

%  Reference: Haykin, Adaptive Filter Theory, p. 356

%------------------
% Check valid input
%------------------
if ((nargin ~= 4) & (nargin ~= 6))
   error('Must have 4 or 6 input arguments.')
end;

[N,tmp] = size(h);
if (tmp ~= 1)
   error('h must be a column vector.');
end;

signal_len = length(x);
signal_len_y = length(y);
if (signal_len ~= signal_len_y)
   error('x and y must have same lengths.')
end;

if (mu_tilde <= 0 | mu_tilde >= 2)
   error('must have 0 < mu_tilde < 2');
end;
if (mu_tilde > 2)
   disp('Warning: mu_tilde exceeds mean-square convergence stability bound.');
   disp('         See IEEE ASSP Feb. 1985 p.225.')
end;

if (nargin == 4)  % no options are specified, use defaults
   start_time = 1;
   stop_time = signal_len+1;
else  % verify valid options
   if ((start_time < 1) | (start_time >= stop_time))
     error('Must have 1 <= start_time < stop_time.');
   end;
   if ((stop_time <= start_time) | (stop_time > signal_len+1))
     error('Must have start_time < stop_time <= signal_len+1.');
   end;
   if (stop_time <= N)
      error('Must have stop_time > N.');
   end;
end;

%-----------
% Initialize
%-----------
a = 0.01;  % positive constant prevents 0 division on calculation of mu
signal_power = 0;

%-------
% Malloc
%-------
e = zeros(signal_len,1);

%----------------
% Adaptive Update
%----------------
% Startup Period (fill adaptive FIR filter) and Possible Adapt Period
for k = 1:N-1
   e(k) = y(k) - h' * [x(k:-1:1).' zeros(1,N-k)].';
   signal_power = x(k) * conj(x(k)) + signal_power;
   if (start_time <= k)
      h = h + mu_tilde / (a + signal_power) * ...
         [x(k:-1:1).' zeros(1,N-k)].' * e(k)';
      start_time = N;
   end;
end;

% Possible Pre-Adapt Period
if (start_time > N)
   e(N) = y(N) - h' * x(N:-1:1);
   signal_power = x(N) * conj(x(N)) + signal_power;
   if (start_time > N+1)
      for k = N+1:start_time-1
         e(k) = y(k) - h' * x(k:-1:k-N+1);
         signal_power = x(k) * conj(x(k)) + signal_power - x(k-N) * ...
            conj(x(k-N));
      end;
   end;
end;

% Adaptation Period
if (start_time == N)
   e(N) = y(N) - h' * x(N:-1:1);
   signal_power = x(N) * conj(x(N)) + signal_power;
   h = h + mu_tilde / (a + signal_power) * x(N:-1:1) * e(N)';
   start_time = N+1;
end;
for k = start_time:stop_time-1
   e(k) = y(k) - h' * x(k:-1:k-N+1);
   signal_power = x(k) * conj(x(k)) + signal_power - x(k-N) * conj(x(k-N));
   h = h + mu_tilde / (a + signal_power) * x(k:-1:k-N+1) * e(k)';
end;

% Possible Post-Adapt Period
if (stop_time <= signal_len)
   for k = stop_time:signal_len
      e(k) = y(k) - h' * x(k:-1:k-N+1);
   end;
end;

