%pmhw.m
%Given d1, d2, extremal frequencies, desired responses, find
%h(k), plot |H(w)| and arg(H(w)) using Parks-McCllelan algorithm.
%d1 = delta 1 = tolerance in the passband, d2 = delta 2 = tolerance in the stopband
%extr = extremal frequencies. Des = desired frequency response
%Example
%>> extr=[0,2*pi/7,3*pi/7,5*pi/7,pi]';Des=[1,1,1,0,0]';
%>> pmhw(0.1,0.2,extr,Des);
function[h,H]=pmhw(d1,d2,extr,Des)
delta=d1;
r=size(extr,1);
q=r-1;
N=2*q+1;
for i=1:r
    if Des(i)>0 W(i)=1;
    else W(i)=d1/d2;
    end
end
%disp(W);
for i=1:r
    for j=1:r
        A(i,j)=cos((j-1)*extr(i));
    end
end
%disp(A);
for i=1:r
    b(i)=Des(i)-(-1)^(i-1)*delta/W(i);
end
%disp(b);
bp=b';
c=A\bp;
disp('c values');
disp(c);
h(1)=c(1);
for i=2:r
    h(i)=c(i)/2;
end
disp('Impulse Response');
disp(h);
%
% Calculating and plotting the frequency response.
%
M = 400;
for n = 1:M+1
   f(n) = (n-1)/M;
   H(n)=c(1);
   q(n) = pi*f(n);   
   
   for k = 2:r
      H(n) = H(n) + c(k)*cos(q(n)*(k-1));
   end

   mag(n) = abs(H(n));
   dB(n) = 20*log10(abs(H(n)));
   if dB(n) < -200 dB(n) = -200; end
   pha(n) = angle(H(n))/pi; 
end
disp('|H(0)|,|H(pi/4)|,|H(pi/2)|,|H(pi)|)');
disp(abs(H(1)));
disp(abs(H(M/4+1)));
disp(abs(H(M/2+1)));
disp(abs(H(3*M/4+1)));
disp(abs(H(M+1)));
figure
plot(f,mag);grid;
title('Magnitude Response')
xlabel('Frequency')
ylabel('Magnitude')
figure
plot(f,dB);grid;
title('Magnitude Response in dB')
xlabel('Frequency')
ylabel('Magnitude in dB')
figure
plot(f,pha);grid;
title('Phase Response')
xlabel('Frequency')
ylabel('Phase in Radians/pi')