Filter design by convex iteration

From Wikimization

(Difference between revisions)
Jump to: navigation, search
(demonstration in Matlab)
Line 55: Line 55:
==demonstration in Matlab==
==demonstration in Matlab==
<pre>
<pre>
-
% Maximize stopband attenuation of a lowpass FIR filter (magnitude design)
+
% CVX by Almir Mutapcic in 2006.
-
% "FIR Filter Design via Spectral Factorization and Convex Optimization"
+
% Adapted in 2010 for impulse response peak-minimization by convex iteration.
-
% by S.-P. Wu, S. Boyd, and L. Vandenberghe
+
%
%
-
% Designs an FIR lowpass filter using spectral factorization method where:
+
% S.-P. Wu, S. Boyd, and L. Vandenberghe,
-
% - minimize maximum stopband attenuation
+
% "FIR Filter Design via Spectral Factorization and Convex Optimization":
-
% - constraint on maximum passband ripple
+
%
%
% minimize max |H(w)| for w in stopband
% minimize max |H(w)| for w in stopband
% s.t. 1/delta <= |H(w)| <= delta for w in passband
% s.t. 1/delta <= |H(w)| <= delta for w in passband
%
%
-
% We change variables via spectral factorization method and get:
+
% Change variables via spectral factorization method and get:
%
%
% minimize max R(w) for w in stopband
% minimize max R(w) for w in stopband
Line 74: Line 72:
% where R(w) is squared magnitude of frequency response
% where R(w) is squared magnitude of frequency response
% (and Fourier transform of autocorrelation coefficients r).
% (and Fourier transform of autocorrelation coefficients r).
-
% Variables are coeffients r and G = hh' where h is impulse response.
+
% Variables are coeffients r and G := hh' where h is impulse response.
% delta is allowed passband ripple.
% delta is allowed passband ripple.
-
% This is a convex problem (can be formulated as an SDP after sampling).
+
% delta2 is specified stopband attenuation.
 +
 
clear all, clc, close all, fclose('all');
clear all, clc, close all, fclose('all');
rand('twister',sum(100*clock));
rand('twister',sum(100*clock));

Revision as of 00:58, 26 August 2010

LaTeX: H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}  where  LaTeX: h(n)\in\mathbb{C}^N denotes impulse response.

For a low pass filter, frequency domain specifications are:

LaTeX: 
\begin{array}{ll}
\frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\
|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
\end{array}

To minimize peak magnitude of LaTeX: h\, , the problem becomes

LaTeX: \begin{array}{cll}
\textrm{minimize} &\|h\|_\infty\\
\textrm{subject~to} &\frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\
&|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
\end{array}

But this problem statement is nonconvex.

So instead,

LaTeX: G\triangleq h(n)h(n)^{\rm H} \in\mathbb{C}^{N\times N}

is a positive semidefinite rank 1 matrix.

Summing along each of LaTeX: 2N-1\, subdiagonals gives entries of the autocorrelation function LaTeX: r\, of LaTeX: h\, .

In particular, the main diagonal of LaTeX: G\, holds squared entries of LaTeX: h\, .

Minimizing LaTeX: \|h\|_\infty is therefore equivalent to minimizing LaTeX: \|\textrm{diag}(G)\|_\infty .


Define LaTeX: I_0\triangleq I\, and define LaTeX: I_n\, as a zero matrix having vector LaTeX: \mathbf{1} along the LaTeX: n^{\rm{th}}\, superdiagonal when LaTeX: n\, is positive or LaTeX: \mathbf{1} along the LaTeX: n^{\rm{th}}\, subdiagonal when LaTeX: n\, is negative.

By spectral factorization, LaTeX: R(\omega)=|H(\omega)|^2\, , an equivalent problem is:

LaTeX: \begin{array}{cll}
\textrm{minimize}_{G\,,\,r}&\|\textrm{diag}(G)\|_\infty\\
\textrm{subject~to} 
& R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2r(n)\cos(\omega n)\\
&\frac{1}{\delta_1^2}\leq R(\omega)\leq\delta_1^2 &\omega\in[0,\omega_p]\\
& R(\omega)\leq\delta_2^2 & \omega\in[\omega_s\,,\pi]\\
& R(\omega)\geq0 & \omega\in[0,\pi]\\
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
& G\succeq0\\
& \textrm{rank}(G) = 1
\end{array}

Excepting the rank constraint, this problem is convex.

demonstration in Matlab

% CVX 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)')
Personal tools