# Filter design by convex iteration

(Difference between revisions)
 Revision as of 22:47, 24 August 2010 (edit)← Previous diff Revision as of 02:59, 25 August 2010 (edit) (undo)Next diff → Line 53: Line 53: Excepting the rank constraint, this problem is convex. Excepting the rank constraint, this problem is convex. + ==demonstration in Matlab==

-                                                           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,:);
-
+ + % 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; -
+                                                              %********************************************************************
-                                                           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)), ...
-
+ [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 $r$ 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)') +

## Revision as of 02: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.

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