function [h,rs,be] = firelp(n,wo,del);
% [h,rs,be] = firelp(n,wo,del);
% Lowpass linear phase FIR filter design with specified ripple sizes
% parameters
%  h   : filter of length 2*n+1
%  rs  : reference set upon convergence
%  be  : band edges
%  wo  : half-magnitude frequency
%  del : [ripple size in passband, stopband] 
% subprograms needed
%  localmax, frefine, firfbe
% % EXAMPLE
%  n  =  17;
%  wo = 0.21*pi;
%  del = [0.03 0.01];
%  [h,rs,be] = firelp(n,wo,del);

% ------------------ initialize some constants ------------------------------
L = 2^ceil(log2(10*n));				% grid size
SN = 1e-12;                                     % Small Number (stopping criterion)
w  = [0:L]*pi/L;                                % frequency axis
n1 = max([1,round(n*wo/pi)]);                   % number of reference frequencies in pass band
n2 = n-n1;                                      % number of reference frequencies in stop band
rs = [wo*[0:n1-1]/n1, wo+(pi-wo)*[1:n2]/n2]';
rs = sort([rs;wo]);                              % initial reference set
Zn = zeros(2*(L+1-n)-3,1);
ideal = [1 0];
up = ideal + del;				% upper bounds
lo = ideal - del;				% lower bounds
Ho = 0.5;					% prescribed value of A(w) at w0
PF = 1;                                 	% PF : flag : Plot Figures

% begin looping
Err = 1;
while Err > SN

   % --------------- calculate interpolation values -------------------------
   Yp = up(1)*(1-(-1).^(1:n1))/2 + lo(1)*((-1).^(1:n1)+1)/2;
   Yp = Yp(n1:-1:1);
   Ys = lo(2)*(1-(-1).^(1:n2))/2 + up(2)*((-1).^(1:n2)+1)/2;
   Y  = [Yp, Ho, Ys];
   
   % --------------- calculate cosine coefficient ---------------------------
   a = cos(rs*[0:n])\Y';
   H = real(fft([a(1);a(2:n+1)/2;Zn;a(n+1:-1:2)/2])); H = H(1:L+1);

   % Plot Figures if PF == 1
   if PF
	figure(1),
	plot(w/pi,H), 
	hold on, 
	plot(rs/pi,Y,'o'), 
	hold off, 
	axis([0 1 -.2 1.2])
   end

   % --------------- determine new reference set-----------------------------
   ri = sort([localmax(H); localmax(-H)]);   
   lr = length(ri);
   while lr > n %  new reference set too big
      if (ri(2)-1)*pi/L > wo
         ri(1) = [];
      elseif (ri(lr-1)-1)*pi/L < wo
         ri(lr) = [];
      else
         x1 = abs(H(ri(1))-H(ri(2)))/del(1);
         x2 = abs(H(ri(lr))-H(ri(lr-1)))/del(2);
         if x1 < x2
            ri(1) = [];
         else
            ri(lr) = [];
         end
      end
      lr = lr - 1;
   end
   rs = (ri-1)*pi/L;
   rs = frefine(a,rs);
   rsp = rs(rs<wo);                % passband reference set frequencies
   rss = rs(rs>wo);                % stopband reference set frequencies
   rs = [rsp; wo; rss];            % new reference set
   n1 = length(rsp);               % number of passband reference set frequencies
   n2 = length(rss);               % number of stopband reference set frequencies

   % --------------- calculate frequency response over reference set --------
   Hp = cos(rsp*[0:n]) * a;
   Hs = cos(rss*[0:n]) * a;

   % --------------- calculate maximum constraint violation -----------------
   Err = max([max(Hp)-up(1), lo(1)-min(Hp), max(Hs)-up(2), lo(2)-min(Hs)]);
   fprintf(1,'    Err = %20.15f\n',Err);
end

% ------------------ calculate band edges -----------------------------------
be = firfbe(a,[wo,wo],[lo(1),up(2)]);

% --------------- calculate interpolation values -------------------------
Yp = up(1)*(1-(-1).^(1:n1))/2 + lo(1)*((-1).^(1:n1)+1)/2;
Yp = Yp(n1:-1:1);
Ys = lo(2)*(1-(-1).^(1:n2))/2 + up(2)*((-1).^(1:n2)+1)/2;
Y  = [Yp, Ho, Ys];

figure(1), plot(w/pi,H), hold on, plot(rs/pi,Y,'o'), axis([0 1 -.2 1.2]), hold off
xlabel('w/pi'), ylabel('H'), title('Frequency Response Amplitude')

h = [a(n+1:-1:2)/2; a(1); a(2:n+1)/2];           % h : filter coefficients


