% "Filter design" lecture notes (EE364) by S. Boyd
% (figures are generated)
%
% Designs a linear phase FIR lowpass filter such that it:
% - minimizes maximum stopband attenuation
% - has a constraint on the maximum passband ripple
%
% This is a convex problem (when sampled it can be represented as an LP).
%
%   minimize   max |H(w)|                     for w in the stopband
%       s.t.   1/delta <= |H(w)| <= delta     for w in the passband
%
% where H is the frequency response function and variable is
% h (the filter impulse response). delta is allowed passband ripple.
%
% Written for CVX by Almir Mutapcic 02/02/06

%********************************************************************
% user's filter specifications
%********************************************************************
% filter order is 2n+1 (symmetric around the half-point)
n = 10;

wpass = 0.12*pi;        % passband cutoff freq (in radians)
wstop = 0.24*pi;        % stopband start freq (in radians)
max_pass_ripple = 1;    % (delta) max allowed passband ripple in dB
                        % ideal passband gain is 0 dB

%********************************************************************
% create optimization parameters
%********************************************************************
N = 30*n;                              % freq samples (rule-of-thumb)
w = linspace(0,pi,N);
A = [ones(N,1) 2*cos(kron(w',[1:n]))]; % matrix of cosines

% passband 0 <= w <= w_pass
ind = find((0 <= w) & (w <= wpass));    % passband
Lp  = 10^(-max_pass_ripple/20)*ones(length(ind),1);
Up  = 10^(max_pass_ripple/20)*ones(length(ind),1);
Ap  = A(ind,:);

% transition band is not constrained (w_pass <= w <= w_stop)

% stopband (w_stop <= w)
ind = find((wstop <= w) & (w <= pi));   % stopband
As  = A(ind,:);

%********************************************************************
% optimization
%********************************************************************
% formulate and solve the linear-phase lowpass filter design
cvx_begin
  variable h(n+1,1);

  minimize( max( abs( As*h ) ) )
  subject to
    % passband bounds
    Lp <= Ap*h;
    Ap*h <= Up;
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
else
  fprintf(1,'The minimum attenuation in the stopband is %3.2f dB.\n\n',...
          20*log10(cvx_optval));
  % construct the full impulse response
  h = [flipud(h(2:end)); h];
end

%********************************************************************
% plots
%********************************************************************
figure(1)
% FIR impulse response
plot([0:2*n],h','o',[0:2*n],h','b:')
xlabel('t'), ylabel('h(t)')

figure(2)
% frequency response
H = exp(-j*kron(w',[0:2*n]))*h;
% magnitude
subplot(2,1,1)
plot(w,20*log10(abs(H)),...
     [0 wpass],[max_pass_ripple max_pass_ripple],'r--',...
     [0 wpass],[-max_pass_ripple -max_pass_ripple],'r--');
axis([0,pi,-50,10])
xlabel('w'), ylabel('mag H(w) in dB')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w'), ylabel('phase H(w)')
 
Calling sedumi: 756 variables, 240 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 240, order n = 757, dim = 757, blocks = 1
nnz(A) = 6948 + 0, nnz(ADA) = 5844, nnz(L) = 3042
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.72E+02 0.000
  1 :  -2.64E-01 1.54E+02 0.000 0.5670 0.9000 0.9000   8.59  1  1  7.7E+01
  2 :  -1.41E-01 6.07E+01 0.000 0.3934 0.9000 0.9000   5.26  1  1  9.9E+00
  3 :  -4.93E-02 3.86E+01 0.000 0.6360 0.9000 0.9000   3.88  1  1  2.9E+00
  4 :  -3.46E-02 2.87E+01 0.000 0.7427 0.9000 0.9000   2.55  1  1  1.7E+00
  5 :  -2.32E-02 1.56E+01 0.000 0.5430 0.9000 0.9000   1.97  1  1  7.5E-01
  6 :  -2.01E-02 7.84E+00 0.000 0.5038 0.9000 0.9000   1.37  1  1  3.5E-01
  7 :  -1.79E-02 2.57E+00 0.000 0.3285 0.9000 0.9000   1.17  1  1  1.1E-01
  8 :  -1.75E-02 8.33E-01 0.000 0.3235 0.9000 0.9000   1.04  1  1  3.5E-02
  9 :  -1.75E-02 1.83E-01 0.000 0.2192 0.9080 0.9000   1.01  1  1  8.1E-03
 10 :  -1.75E-02 3.45E-02 0.000 0.1892 0.9055 0.9000   1.00  1  1  1.6E-03
 11 :  -1.75E-02 4.21E-03 0.000 0.1219 0.9088 0.9000   1.00  1  1  2.0E-04
 12 :  -1.75E-02 2.19E-05 0.000 0.0052 0.9990 0.9990   1.00  1  1  
iter seconds digits       c*x               b*y
 12      0.1  15.4 -1.7476196636e-02 -1.7476196636e-02
|Ax-b| =   3.2e-16, [Ay-c]_+ =   5.0E-17, |x|=  6.1e-01, |y|=  3.2e-01

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    5.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1.122018e+00,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.33849.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.0174762
Problem is Solved
The minimum attenuation in the stopband is -35.15 dB.