# Filter design by convex iteration

### From Wikimization

(Difference between revisions)

(→demonstration in Matlab) |
(→demonstration in Matlab) |
||

Line 55: | Line 55: | ||

==demonstration in Matlab== | ==demonstration in Matlab== | ||

<pre> | <pre> | ||

- | % CVX by Almir Mutapcic in 2006. | + | % CVX code by Almir Mutapcic in 2006. |

% Adapted in 2010 for impulse response peak-minimization by convex iteration. | % Adapted in 2010 for impulse response peak-minimization by convex iteration. | ||

% | % |

## Revision as of 02:30, 17 October 2010

where denotes impulse response.

For a low pass filter, frequency domain specifications are:

To minimize peak magnitude of , the problem becomes

But this problem statement is nonconvex.

So instead,

is a positive semidefinite rank 1 matrix.

Summing along each of subdiagonals gives entries of the autocorrelation function of .

In particular, the main diagonal of holds squared entries of .

Minimizing is therefore equivalent to minimizing .

Define and define as a zero matrix having vector along the superdiagonal when is positive or along the subdiagonal when is negative.

By spectral factorization, , an equivalent problem is:

Excepting the rank constraint, this problem is convex.

## demonstration in Matlab

% CVX code by Almir Mutapcic in 2006. % Adapted in 2010 for impulse response peak-minimization by convex iteration. % % S.-P. Wu, S. Boyd, and L. Vandenberghe, % "FIR Filter Design via Spectral Factorization and Convex Optimization": % % minimize max |H(w)| for w in stopband % s.t. 1/delta <= |H(w)| <= delta for w in passband % % Change variables via spectral factorization method and get: % % minimize max R(w) for w in stopband % s.t. (1/delta)^2 <= R(w) <= delta^2 for w in passband % R(w) >= 0 for all w % % where R(w) is squared magnitude of frequency response % (and Fourier transform of autocorrelation coefficients r). % Variables are coeffients r and G := hh' where h is impulse response. % delta is allowed passband ripple. % delta2 is specified stopband attenuation. clear all, clc, close all, fclose('all'); rand('twister',sum(100*clock)); randn('state',sum(100*clock)); %********************************************************************* % filter specs (for a low-pass filter) %********************************************************************* % number of FIR coefficients (including zeroth) N = 32; wpass = 0.12*pi; % end of the passband wstop = 0.20*pi; % start of the stopband delta = 1; % maximum passband ripple in dB (+/- around 0 dB) delta2 = -20; % stopband attenuation desired in dB %********************************************************************* % optimization parameters %********************************************************************* % rule-of-thumb discretization (from Cheney's Approximation Theory) m = 15*N; w = linspace(0,pi,m)'; % omega % A is the matrix used to compute the power spectrum % A(w,:) = [1 2*cos(w) 2*cos(2*w) ... 2*cos(N*w)] A = [ones(m,1) 2*cos(kron(w,[1:N-1]))]; % passband 0 <= w <= w_pass ind = find((0 <= w) & (w <= wpass)); % passband Lp = 10^(-delta/20)*ones(length(ind),1); Up = 10^(+delta/20)*ones(length(ind),1); Ap = A(ind,:); % stopband (w_stop <= w) ind = find((wstop <= w) & (w <= pi)); % stopband Sp = 10^(delta2/20)*ones(length(ind),1); As = A(ind,:); % make I matrices B = zeros(N,N^2); for i=0:N-1 C = zeros(N,N); C = spdiags(ones(N,1),i,C); B(i+1,:) = vect(C)'; end %initial direction vector W = eye(N); % peak impulse response peak = 0.125; %******************************************************************** % optimization %******************************************************************** % formulate and solve the magnitude design problem iteration = 1; while 1 tic fprintf('\nMinimizing impulse response peak: iteration %d\n', iteration); cvx_quiet(true); % cvx_solver('sdpt3'); cvx_begin variable r(N,1); variable G(N,N) symmetric; minimize(trace(W'*G)); % passband constraints Ap*r >= Lp.^2; Ap*r <= Up.^2; %stopband constraint As*r <= Sp.^2; % nonnegative-real constraint A*r >= 0; % relate r to h r == B*vect(G); G == semidefinite(N); %minimize peak of h diag(G) <= peak^2; cvx_end toc %computer new direction vector [v,d,q] = svd(G); f = diag(d); fprintf('first few eigenvalues of G:\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n', f(1:7)); W = v(:,2:N)*v(:,2:N)'; % W = v(:,2:N)*diag(rand(N-1,1))*v(:,2:N)'; rankG = sum(diag(d) > max(diag(d))*1e-5); fprintf('rank(G)=%d, trace(W*G)=%f\n', rankG, trace(G*W)); if (rankG == 1) break end figure(1) % FIR impulse response h = G(:,1)/sqrt(G(1,1)); plot([0:N-1],h','o',[0:N-1],h','b:') xlabel('t'), ylabel('h(t)') pause(1) % check if problem was successfully solved disp(['Problem is ' cvx_status]) if ~strfind(cvx_status,'Solved') return end iteration = iteration + 1; end % compute the min attenuation in the stopband (convert to original vars) Ustop = delta2; fprintf(1,'Min attenuation in the stopband is %3.2f dB.\n',Ustop); %********************************************************************* % plotting routines %********************************************************************* % frequency response of the designed filter, where j = sqrt(-1) h = G(:,1)/sqrt(G(1,1)); H = [exp(-j*kron(w,[0:N-1]))]*h; figure(2) % magnitude plot(w,20*log10(abs(H)), ... [0 wpass],[delta delta],'r--', ... [0 wpass],[-delta -delta],'r--', ... [wstop pi],[Ustop Ustop],'r--') xlabel('w') ylabel('mag H(w) in dB') axis([0 pi -50 5]) %compare impulse response designed by conventional method figure(3) h = spectral_fact(r); %from CVX distribution, Examples subdirectory plot([0:N-1],h','o',[0:N-1],h','b:') xlabel('t'), ylabel('h(t)') figure(1) % FIR impulse response h = G(:,1)/sqrt(G(1,1)); plot([0:N-1],h','o',[0:N-1],h','b:') xlabel('t'), ylabel('h(t)')