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

Jump to: Source code Text output Plots Library index

    % "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