Filter design by convex iteration

From Wikimization

(Difference between revisions)
Jump to: navigation, search
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>
-
N = 16;
+
% Maximize stopband attenuation of a lowpass FIR filter (magnitude design)
-
delta1 = 1.01;
+
% "FIR Filter Design via Spectral Factorization and Convex Optimization"
-
delta2 = 0.01;
+
% by S.-P. Wu, S. Boyd, and L. Vandenberghe
-
w = linspace(0,pi,N);
+
%
-
DFT = real(fft(eye(2*N))); % generate DFT matrix but only use the real components
+
% Designs an FIR lowpass filter using spectral factorization method where:
-
% https://ccrma.stanford.edu/~jos/st/DFT_Matrix.html
+
% - 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;
-
cvx_begin
+
wpass = 0.12*pi; % end of the passband
-
variable r(2*N-1,1);
+
wstop = 0.20*pi; % start of the stopband
-
variable G(N^2,N^2);
+
delta = 1; % maximum passband ripple in dB (+/- around 0 dB)
 +
delta2 = -20; % stopband attenuation desired in dB
-
r = [0; r]; % otherwise r is not even
+
%*********************************************************************
-
+
% optimization parameters
-
+
%*********************************************************************
-
minimize(max(diag(G)));
+
% rule-of-thumb discretization (from Cheney's Approximation Theory)
 +
m = 15*N;
 +
w = linspace(0,pi,m)'; % omega
-
G2 = G(1:N,1:N); % I don't really need to use the entire G.
+
% A is the matrix used to compute the power spectrum
-
% So just use the 1st block: G2
+
% A(w,:) = [1 2*cos(w) 2*cos(2*w) ... 2*cos(N*w)]
-
for I=1:N
+
A = [ones(m,1) 2*cos(kron(w,[1:N-1]))];
-
temp = diag(ones(1,I),-(N-I));
+
-
[i,j,k] = find(temp);
+
-
index = [i j ones(length(i),1)];
+
-
index = [index; [N, N, 0]];
+
-
temp2 = spconvert(index);
+
-
r(I) == trace(temp2*G2);
+
-
end
+
-
+
-
for I=N+1:2*N-1
+
-
temp = diag(ones(1,2*N-I),I-N);
+
-
[i,j,k] = find(temp);
+
-
index = [i j ones(length(i),1)];
+
-
index = [index; [N, N, 0]];
+
-
temp2 = spconvert(index);
+
-
r(I) == trace(temp2*G2);
+
-
end
+
-
+
-
R = fftshift(DFT*r);
+
-
+
-
for I=1:N %%% R is of length of 2N
+
-
1/delta1^2 < R(I);
+
-
R(I) < delta1^2;
+
-
end
+
-
for I=N+1:2*N
+
% passband 0 <= w <= w_pass
-
R(I) < delta2^2;
+
ind = find((0 <= w) & (w <= wpass)); % passband
-
end
+
Lp = 10^(-delta/20)*ones(length(ind),1);
 +
Up = 10^(+delta/20)*ones(length(ind),1);
 +
Ap = A(ind,:);
-
R > 0;
+
% stopband (w_stop <= w)
-
+
ind = find((wstop <= w) & (w <= pi)); % stopband
-
+
Sp = 10^(delta2/20)*ones(length(ind),1);
-
cvx_end
+
As = A(ind,:);
-
</pre>
+
 +
% 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
-
But the above doesn't work.
+
%initial direction vector
 +
W = eye(N);
-
Debug: try the following:
+
% peak impulse response
 +
peak = 0.125;
-
<pre>
+
%********************************************************************
-
N = 16;
+
% optimization
-
delta1 = 1.01;
+
%********************************************************************
-
delta2 = 0.01;
+
% formulate and solve the magnitude design problem
-
w = linspace(0,pi,N);
+
iteration = 1;
-
DFT = real(fft(eye(2*N)));
+
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
-
cvx_begin
+
iteration = iteration + 1;
-
variable r(2*N-1,1);
+
end
-
r = [0; r];
+
% 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);
-
R = fftshift(DFT*r);
+
%*********************************************************************
-
+
% plotting routines
-
for I=1:N
+
%*********************************************************************
-
1/delta1^2 <= R(I);
+
% frequency response of the designed filter, where j = sqrt(-1)
-
R(I) <= delta1^2;
+
h = G(:,1)/sqrt(G(1,1));
-
end
+
H = [exp(-j*kron(w,[0:N-1]))]*h;
-
+
-
for I=N+1:2*N
+
-
R(I) <= delta2^2;
+
-
end
+
-
R > 0;
+
figure(2)
-
+
% magnitude
-
cvx_end
+
plot(w,20*log10(abs(H)), ...
-
</pre>
+
[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])
-
The above just to see if there exists an <math> r </math> such that the problem is feasible.
+
%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)')
-
But it still doesn't work...
+
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 01:59, 25 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

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