% 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')