The linear phase [Math Processing Error] FIR algorithm

The following program designs a Type-I length-[Math Processing Error] linear phase [Math Processing Error] filter with fixed transition bands. The code creates even sampling in the bands but it can easily be modified to accomodate for uneven sampling. For other linear phase types the user can modify the definition of C according to (Reference).

% Lp Linear Phase design program. Uses p-homotopy and partial updating.
 
% Author: R. A. Vargas, 4-20-08
 
%%% User parameters
P = 400;                         % Desired p
K = 1.7;                         % Homotopy rate
fp = 0.2; fs = 0.24;             % Normalized passband & stopband
L  = 21;                         % Filter length
 
%%% 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;                        % L2 initial guess
e = C*a - Ad;                    % Initial error
 
%%% Algorithm iteration
c = 1; i = 40; p = 2;
while (c<i),
    c = c+1;
    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
    e = C*a - Ad;                % Error
end
 
%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;