Filter design by convex iteration
From Wikimization
(Difference between revisions)
Line 53: | Line 53: | ||
Excepting the rank constraint, this problem is convex. | Excepting the rank constraint, this problem is convex. | ||
+ | ==demonstration in Matlab== | ||
<pre> | <pre> | ||
- | + | % Maximize stopband attenuation of a lowpass FIR filter (magnitude design) | |
- | + | % "FIR Filter Design via Spectral Factorization and Convex Optimization" | |
- | + | % by S.-P. Wu, S. Boyd, and L. Vandenberghe | |
- | w = | + | % |
- | + | % Designs an FIR lowpass filter using spectral factorization method where: | |
- | + | % - minimize maximum stopband attenuation | |
+ | % - constraint on maximum passband ripple | ||
+ | % | ||
+ | % minimize max |H(w)| for w in stopband | ||
+ | % s.t. 1/delta <= |H(w)| <= delta for w in passband | ||
+ | % | ||
+ | % We 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. | ||
+ | % This is a convex problem (can be formulated as an SDP after sampling). | ||
+ | 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; | ||
- | + | %******************************************************************** | |
- | N = | + | % 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)') | ||
+ | </pre> |
Revision as of 02:59, 25 August 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
% Maximize stopband attenuation of a lowpass FIR filter (magnitude design) % "FIR Filter Design via Spectral Factorization and Convex Optimization" % by S.-P. Wu, S. Boyd, and L. Vandenberghe % % Designs an FIR lowpass filter using spectral factorization method where: % - minimize maximum stopband attenuation % - constraint on maximum passband ripple % % minimize max |H(w)| for w in stopband % s.t. 1/delta <= |H(w)| <= delta for w in passband % % We 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. % This is a convex problem (can be formulated as an SDP after sampling). 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)')