Design a 1/f spectrum shaping (pink-noise) filter

% "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: 3000 variables, 641 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 = 641, order n = 1803, dim = 3602, blocks = 602
nnz(A) = 27149 + 0, nnz(ADA) = 51481, nnz(L) = 26061
Handling 1 + 1 dense columns.
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            5.19E+01 0.000
  1 :  -5.87E+00 2.04E+01 0.000 0.3939 0.9000 0.9000  -0.97  1  1  5.9E+02
  2 :  -3.67E+01 3.72E+00 0.000 0.1819 0.9000 0.9000  -0.86  1  1  4.1E+02
  3 :  -5.23E+01 1.81E-01 0.000 0.0488 0.9900 0.9900   0.23  1  1  2.5E+01
  4 :  -9.93E+00 8.50E-02 0.000 0.4687 0.9000 0.9000   8.65  1  1  2.3E+00
  5 :  -2.83E+00 3.12E-02 0.000 0.3665 0.9000 0.9000   3.70  1  1  3.4E-01
  6 :  -1.65E+00 1.14E-02 0.000 0.3645 0.9000 0.9000   1.34  1  1  1.1E-01
  7 :  -1.31E+00 4.78E-03 0.000 0.4205 0.9000 0.9000   0.96  1  1  5.0E-02
  8 :  -1.22E+00 2.30E-03 0.000 0.4809 0.6154 0.9000   0.95  1  1  2.4E-02
  9 :  -1.21E+00 1.28E-03 0.000 0.5593 0.9000 0.6822   1.00  1  1  1.4E-02
 10 :  -1.19E+00 3.98E-04 0.000 0.3102 0.9000 0.9000   1.00  1  1  4.2E-03
 11 :  -1.19E+00 1.31E-04 0.000 0.3291 0.9000 0.9000   1.00  1  1  1.4E-03
 12 :  -1.19E+00 4.30E-05 0.000 0.3281 0.9000 0.9000   1.00  1  1  4.6E-04
 13 :  -1.19E+00 1.17E-05 0.000 0.2713 0.9024 0.9000   1.00  1  1  1.2E-04
 14 :  -1.19E+00 2.32E-06 0.000 0.1987 0.9084 0.9000   1.00  1  1  1.9E-05
 15 :  -1.19E+00 5.26E-07 0.000 0.2270 0.9088 0.9000   1.00  1  1  3.6E-06
 16 :  -1.19E+00 4.17E-08 0.465 0.0792 0.9900 0.9901   1.00  1  1  3.1E-07
 17 :  -1.19E+00 1.02E-08 0.000 0.2451 0.9000 0.9097   1.00  1  1  7.9E-08
 18 :  -1.19E+00 1.41E-10 0.000 0.0138 0.9901 0.9900   1.00  2  2  1.0E-09

iter seconds digits       c*x               b*y
 18      1.4   9.6 -1.1873310834e+00 -1.1873310837e+00
|Ax-b| =   5.0e-10, [Ay-c]_+ =   1.2E-11, |x|=  1.2e+01, |y|=  2.6e+02

Detailed timing (sec)
   Pre          IPM          Post
5.000E-02    1.390E+00    1.000E-02    
Max-norms: ||b||=3.183099e-01, ||c|| = 2,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 6.28317.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.18733
Problem is Solved