% "Filter design" lecture notes (EE364) by S. Boyd
% "FIR filter design via spectral factorization and convex optimization"
% by S.-P. Wu, S. Boyd, and L. Vandenberghe
% (a figure is generated)
%
% Designs a log-Chebychev filter magnitude design given as:
%
%   minimize   max| log|H(w)| - log D(w) |   for w in [0,pi]
%
% where variables are impulse response coefficients h, and data
% is the desired frequency response magnitude D(w).
%
% We can express and solve the log-Chebychev problem above as
%
%   minimize   max( R(w)/D(w)^2, D(w)^2/R(w) )
%       s.t.   R(w) = |H(w)|^2   for w in [0,pi]
%
% where we now use the auto-correlation coeffients r as variables.
%
% As an example we consider the 1/sqrt(w) spectrum shaping filter
% (the so-called pink-noise filter) where D(w) = 1/sqrt(w).
% Here we use a logarithmically sampled freq range w = [0.01*pi,pi].
%
% Written for CVX by Almir Mutapcic 02/02/06

% parameters
n = 40;      % filter order
m = 15*n;    % frequency discretization (rule-of-thumb)

% log-space frequency specification
wa = 0.01*pi; wb = pi;
wl = logspace(log10(wa),log10(wb),m)';

% desired frequency response (pink-noise filter)
D = 1./sqrt(wl);

% matrix of cosines to compute the power spectrum
Al = [ones(m,1) 2*cos(kron(wl,[1:n-1]))];

% solve the problem using cvx
cvx_begin
  variable r(n,1)   % auto-correlation coefficients
  variable R(m,1)   % power spectrum

  % log-chebychev minimax design
  minimize( max( max( [R./(D.^2)  (D.^2).*inv_pos(R)]' ) ) )
  subject to
     % power spectrum constraint
     R == Al*r;
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

% spectral factorization
h = spectral_fact(r);

% figures
figure(1)
H = exp(-j*kron(wl,[0:n-1]))*h;
loglog(wl,abs(H),wl,D,'r--')
set(gca,'XLim',[wa pi])
xlabel('freq w')
ylabel('mag H(w) and D(w)')
legend('optimized','desired')
 
Calling sedumi: 3641 variables, 2400 equality constraints
------------------------------------------------------------
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
Split 41 free variables
eqs m = 2400, order n = 3083, dim = 4283, blocks = 601
nnz(A) = 5400 + 49200, nnz(ADA) = 9600, nnz(L) = 6000
Handling 82 + 0 dense columns.
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.04E-01 0.000
  1 :   3.49E+01 8.25E-02 0.000 0.2042 0.9000 0.9000  -0.98  1  1  6.4E+01
  2 :   2.06E+02 7.52E-03 0.000 0.0911 0.9900 0.9900  -0.82  1  1  2.8E+01
  3 :   9.16E+01 3.21E-03 0.000 0.4272 0.9000 0.9000   1.69  1  1  7.7E+00
  4 :   3.04E+01 9.91E-04 0.000 0.3088 0.9000 0.9000   3.48  1  1  8.0E-01
  5 :   8.04E+00 4.32E-04 0.000 0.4361 0.9000 0.9000   4.95  1  1  1.0E-01
  6 :   3.44E+00 2.23E-04 0.000 0.5151 0.9000 0.9000   2.91  1  1  3.2E-02
  7 :   2.17E+00 1.01E-04 0.000 0.4548 0.9000 0.9000   1.53  1  1  1.3E-02
  8 :   1.68E+00 5.00E-05 0.000 0.4938 0.9000 0.9000   1.08  1  1  6.4E-03
  9 :   1.45E+00 3.11E-05 0.000 0.6218 0.9000 0.9000   0.92  1  1  4.1E-03
 10 :   1.27E+00 1.52E-05 0.000 0.4900 0.9000 0.9000   0.91  1  1  2.1E-03
 11 :   1.23E+00 8.33E-06 0.000 0.5466 0.9000 0.9000   1.02  1  1  1.1E-03
 12 :   1.20E+00 2.57E-06 0.000 0.3089 0.9016 0.9000   1.00  1  1  3.9E-04
 13 :   1.19E+00 7.43E-07 0.000 0.2887 0.9044 0.9000   1.00  1  1  1.2E-04
 14 :   1.19E+00 2.24E-07 0.000 0.3009 0.9012 0.9000   1.00  1  1  3.8E-05
 15 :   1.19E+00 2.61E-08 0.000 0.1168 0.9000 0.0000   1.00  1  1  1.3E-05
 16 :   1.19E+00 4.21E-09 0.000 0.1611 0.9292 0.9000   1.00  5  5  3.5E-06
 17 :   1.19E+00 2.24E-10 0.176 0.0532 0.9900 0.8463   1.00  6  6  3.3E-07
 18 :   1.19E+00 3.95E-12 0.000 0.0176 0.9900 0.9900   1.00 14 14  6.2E-09

iter seconds digits       c*x               b*y
 18      8.5   Inf  1.1873313477e+00  1.1873324283e+00
|Ax-b| =   2.5e-08, [Ay-c]_+ =   2.1E-09, |x|=  2.7e+02, |y|=  1.1e+00

Detailed timing (sec)
   Pre          IPM          Post
4.000E-02    8.540E+00    3.000E-02    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=7, |skip| = 2, ||L.L|| = 3.72784.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.18733
Problem is Solved