The adaptive linear phase [Math Processing Error] FIR algorithm

The following code expands on the program from "The linear phase [Math Processing Error] FIR algorithm". If at a given iteration the error increases, the idea is to take a step back in[Math Processing Error] and then take a smaller step. This is achieved by reducing the homotopy step parameter K in the program.

% Lp Linear Phase design program. Uses p-homotopy and partial updating.
% This code uses the adaptive algorithm.
 
% Author: R. A. Vargas, 4-20-08
 
%%% User parameters
P = 200;                         % Desired p
K = 1.7;                         % Homotopy rate
fp = 0.2; fs = 0.248;            % Normalized passband & stopband
L  = 21;                         % Filter length
dk = 0.9;                        % K update factor
 
%%% Initialization
NF = 2^ceil(log2(10*L));         % Number of frequency samples
Np = round(NF*fp/(.5 - (fs - fp))); % No. of passband samples
dp = fp/(Np - 1);
Ns = NF - Np;                    % No. of stopband samples
ds = (0.5 - fs)/(Ns - 1);
fd = [(0:Np - 2)*dp,fp,(0:Ns - 2)*ds + fs,0.5]'; % Frequencies
Ad = [ones(Np,1); zeros(Ns,1)];  % Desired response
M = floor((L+1)/2)-1;
C = cos(2*pi*fd*(0:M));          % Fourier matrix
a = C\Ad; a0 = a;                % L2 initial guess
e = C*a - Ad; en = e;            % Initial error
 
%%% Algorithm iteration
c = 1; maxiter = 40; p = 2;
while (c<maxiter),
   p = min(P,K*p);               % p-homotopy update
   w = abs(e).^((p-2)/2);        % Lp weight update
   W = diag(w/sum(w));           % Normalized weights
   x = (W*C)\(W*Ad);             % WLS solution
   a = (x + (p-2)*a) ./(p-1);    % Partial update
   en = C*a - Ad;                % Error
   if (norm(en,P) <= norm(e,P)) | p>=P
      c = c+1;
      e = en;
      a0 = a;
   else
      a = a0;
      p = p/K;
      K = dk*K;                  % Update homotopy step
   end
end
 
%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;