% Linear Phase Filter Design using the Chebychev Criterion
% min max_{k=1,...,K} | H(omega_k) - H_des(omega_k) |
%
% by Wing-Kin Ma, Dec, 2008

clear all; close all;

n= 200;  %filter size
K= 1000; %no. of freq. sample pts

omega= (0:(pi/(K-1)):pi)';  % set of freq. sample pts.

%--- Set H_des(omega) to have a pass response on [0, omega_p], and 
% a stop response on [ omega_s, pi ]
omega_p= .25*pi;  %passband
omega_s= .28*pi;  %stopband
H_des= zeros(K,1);
ind= find( omega <= omega_s );
H_des(ind)= 1;
% linearly decreasing response for [ omega_p, omega_s ]
ind= find( omega >= omega_p & omega <= omega_s );
ell= length(ind);
H_des(ind(1)+(0:1:(ell-1)))= 1- (1:1:ell)/ell;  

%H_des= H_des+ 10^-2;

%--- CVX setup

A= zeros(K,n+1);
for k=1:K
    A(k,:)= [ 1 2*cos(omega(k)*(1:n)) ];
end;
cvx_begin
    variable h(n+1);
    minimize( norm(A*h- H_des, Inf) );
cvx_end
    

%--- plot the result
H= zeros(K,1);  %freq. response of the designed filter
for k=1:K
    H(k)= [ 1; 2*cos(omega(k)*(1:n)') ]'*h;
end;

figure(1); clf;
plot(omega,H,'linewidth',2); hold on
plot(omega,H_des,'r','linewidth',2);
xlabel('\omega'); ylabel('Frequency Response');
legend('Desired Response','Chebychev FIR Response');

figure(2); clf; hold on;
plot(omega,20*log10(abs(H)),'linewidth',1); hold on;
plot(omega,20*log10(10^(-10)+H_des),'k','linewidth',2);

xlabel('w','fontsize',16); ylabel('H(w) (dB)','fontsize',16);
set(gca,'XLim',[0,omega(end)]);
set(gca,'YLim',[-40,5]);
legend('Chebychev Design','Desired Response');

% hold on;
% plot([omega_s omega_s],[-70,10],'--k','linewidth',1);
% plot([omega(end-1) omega(end-1)],[-70,10],'--k','linewidth',1);
% plot([omega_s omega(end-1)],[-20,-20],'-k','linewidth',2);


