Filter design by convex iteration

From Wikimization

(Difference between revisions)

Revision as of 07:51, 1 November 2010

1 RF pulse design with small flip angle
2 RF pulse design with large flip angle
3 A typical design procedure involves the following steps:
4 Minimum peak RF waveform for saturation pulse
5 Minimum peak time-domain response by Optimization
6 Matlab demonstration

RF pulse design with small flip angle

Flip angle is the amount of deviation one applies to net magnetization vector from its transverse axis. A $LaTeX: 90^o$ pulse will tip the net magnetization vector to the longitudinal axis. If the desired flip angle $LaTeX: \theta$ is small $LaTeX: (\sin(\theta)\simeq\theta)$ , then the time domain RF waveform and the pulse response can be approximated by Fourier Transform.

RF pulse design with large flip angle

For design RF pulses with large flip angles, a technique called Shinnar-Le Roux (SLR) transform is often used. SLR transform relates the desired pulse response to the design of FIR filters.

The goal of RF pulse design is: find the time domain waveform $LaTeX: B_1(t)$ that will produce the desired transverse magnetization pulse profile and longitudinal magnetization pulse profile. SLR transform relates the RF waveform to two complex polynomials $LaTeX: A_N(z)$ and $LaTeX: B_N(z)$ .

i.e. $LaTeX: B_1(t)\begin{array}{c}\mbox{SLR}\\\Leftrightarrow \end{array}\, A_N(z),\, B_N(z)$

Thus, the task of filter design becomes finding these two polynomials, which can be done via FIR filter design methods.

A typical design procedure involves the following steps:

1. Establish a set of design parameters: e.g. waveform duration, pulse bandwidth, passband and stopband ripples. These parameters are converted to their FIR filter design counterparts.

2. Use Parks-McClellan algorithm to come up with a waveform, $LaTeX: B_N(z)$ , where $LaTeX: z=e^{j\gamma\,G_x\,\Delta\,t}$

3. $LaTeX: A_N(z)$ and $LaTeX: B_N(z)$ are related by $LaTeX: |A_N(z)| = \sqrt{ 1-B_N(z)\,B_N^\ast\,(z)}$ . If $LaTeX: A_N(z)$ is chosen to be a minimum-phase polynomial, then it is also an analytic signal. Analytic signals have the property that their log-magnitude and phase are a Hilbert transform pair. For the given $LaTeX: B_N(z)$ , the unique minimum-phase $LaTeX: A_N(z)$ is $LaTeX: A_N(z)=|A_N(z)|\,e^{i\mathcal{H}(\mathrm{log}|A_N(z)|)}$ , where $LaTeX: \mathcal{H}$ is the Hilbert transform opetator.

4. Once $LaTeX: A_N\,(z)$ and $LaTeX: B_N\,(z)$ are found, RF waveform $LaTeX: B_1(t)$ can be computed by inverse SLR transform.

Note: upon applying this RF pulse, the resulting transverse magnetization is $LaTeX: |M_{xy}(\omega)| = 2\sqrt{1-|B_N\,(\omega)|^2}\,|B_N\,(\omega)|\,M_o$ where $LaTeX: z=e^{j\omega}$ and $LaTeX: M_o$ is the initial magnetization. If $LaTeX: |B_N\,(\omega)|$ represents a rectangular profile, $LaTeX: |M_{xy}(\omega)|$ will also has a rectangular profile.

Minimum peak $LaTeX: B_1(t)$ RF waveform for saturation pulse

Our goal here is to find a minimum peak magnitude $LaTeX: B_1(t)$ for a given magnitude profile response. Since the application is for saturation, the phase of the profile response is not important but the passband should be flat and transition region should be sharp.

Plan: follow the above 4 steps. Except in step 2, instead of using Parks-McClellan algorithm to find $LaTeX: B_N(z)$ , use optimization:

Minimum peak time-domain response by Optimization

$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\,$ , where $LaTeX: r\,\triangleq r_r+ir_i\,\,\in\mathbb{C}^N$ .

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}2(r_r(n)\cos(\omega n)+r_i(n)\sin(\omega n))\\ &(1-\delta_1)^2\leq R(\omega)\leq(1+\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. The technique of convex iteration (to find direction vector $LaTeX: W\,$ ) from Dattorro's book is applied. Regularization weight $LaTeX: \lambda\,$ is chosen by cut-and-try as in the figure.

$LaTeX: \begin{array}{cll} \textrm{minimize}_{G\,,\,r}&\langle G\,,W\rangle\lambda+\|\textrm{diag}(G)\|_\infty\\ \textrm{subject~to} & R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2(r_r(n)\cos(\omega n)+r_i(n)\sin(\omega n))\\ &(1-\delta_1)^2\leq R(\omega)\leq(1+\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 \end{array}$

Matlab demonstration

% CVX code by Almir Mutapcic in 2006.  
% Adapted in 2010 for impulse response peak-minimization by convex iteration by Christine Law.
%
% "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 with
% constraint on maximum passband ripple and stopband attenuation:
%
%   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 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;

lambda = 0.5;  %weighting for convex iteration. Want this as small as possible.  
               %If no convergence, make larger. See Optimization book by Dattorro ch.4.4.1.3

wpass = 0.12*pi;   % end of passband
wstop = 0.20*pi;   % start of stopband
delta0_wpass = 0.125; 
delta0_wstop = 0.125; 
delta  = 20*log10(1 + delta0_wpass)  % maximum passband ripple in dB (+/- around 0 dB)
delta2 = 20*log10(delta0_wstop)      % 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_p = find((0 <= w) & (w <= wpass));    % passband
Lp  = 10^(-delta/20)*ones(length(ind_p),1);
Up  = 10^(+delta/20)*ones(length(ind_p),1);
Ap  = A(ind_p,:);

% stopband (w_stop <= w)
ind_s = find((wstop <= w) & (w <= pi));   % stopband
Sp  = 10^(delta2/20)*ones(length(ind_s),1);
As  = A(ind_s,:);

%remove redundant contraints
ind_nr = setdiff(1:m,ind_p);   % fullband less passband
Anr = A(ind_nr,:);

% 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);

%********************************************************************
% optimization
%********************************************************************
convergence = [];
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)*lambda + max(diag(G)));
      % passband constraints
      Ap*r >= Lp.^2;
      Ap*r <= Up.^2;
      %stopband constraint
      As*r <= Sp.^2;
      % nonnegative-real constraint
      Anr*r >= 0;
      % relate r to h
      r == B*vect(G);
      G == semidefinite(N);
   cvx_end
   toc
   
   %compute 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)';  %Use this if convex iteration stalling alot.  
                                               %Stalling will also occur if lambda too small.
   rankG = sum(diag(d) > max(diag(d))*1e-5);
   fprintf('rank(G)=%d,  trace(W*G)=%f\n', rankG, trace(G*W)); 
   
   figure(1)
   % FIR impulse response
   h = G(:,1)/sqrt(G(1,1));
   plot([0:N-1],h','ob:')
   xlabel('t'), ylabel('h(t)')

   % monitor convergence to 0
   convergence = [convergence trace(G*W)];
   if iteration > 1
      figure(4)
      plot(convergence)
      set(gcf,'position',[70 200 256 256])
   end
   pause(1)
   
   % check if problem was successfully solved
   disp(['Problem is ' cvx_status])
   if (rankG == 1)
      break
   end
   if ~strfind(cvx_status,'Solved')
       fprintf(2,'Excuse me.\n')
   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);
subplot(121),
% 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])
title(sprintf('N, wp(pi), ws(pi), delta = %d  %3.2f  %3.2f  %3.2f',N, wpass/pi, wstop/pi, delta));

%compare impulse response designed by conventional method
subplot(122),
h_sp = spectral_fact(r);  %from CVX distribution, Examples subdirectory
plot([0:N-1],h_sp','ob:');
hold on;
h = G(:,1)/sqrt(G(1,1));
plot([0:N-1],h','+r--');
legend('conventional','optimal');
xlabel('t'), ylabel('h(t)'); grid
title(sprintf('h max conventional and optimal; lambda = %3.4f  %3.4f  %3.2f',max(abs(h_sp)),max(abs(h)),lambda));
set(gcf,'Outerposition',[300 300 256*4 256*2])
 
figure(1)
% FIR impulse response
h = G(:,1)/sqrt(G(1,1));
plot([0:N-1],h','ob:'); 
xlabel('t'), ylabel('h(t)')

Retrieved from "http://www.convexoptimization.com/wikimization/index.php/Filter_design_by_convex_iteration"

Filter design by convex iteration

From Wikimization

Revision as of 07:51, 1 November 2010

Contents

RF pulse design with small flip angle

RF pulse design with large flip angle

A typical design procedure involves the following steps:

Minimum peak $LaTeX: B_1(t)$ RF waveform for saturation pulse

Minimum peak time-domain response by Optimization

Matlab demonstration

Views

Personal tools

Navigation

Search

Toolbox

@@ Line 68: / Line 68: @@
 is a positive semidefinite rank 1 matrix.
-Summing along each of <math>2N-1\,</math> subdiagonals gives entries of the autocorrelation function <math>r\,</math> of <math>h\,</math>&nbsp;.
+Summing along each of <math>2N-1\,</math> subdiagonals gives entries of the autocorrelation function <math>r\,</math> of <math>h\,</math>, where <math>r\,\triangleq r_r+ir_i\,\,\in\mathbb{C}^N</math>&nbsp;.
 In particular, the main diagonal of <math>G\,</math> holds squared entries of <math>h\,</math>&nbsp;.
@@ Line 83: / Line 83: @@
 \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)\\
+& R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2(r_r(n)\cos(\omega n)+r_i(n)\sin(\omega n))\\
-&\frac{1}{\delta_1^2}\leq R(\omega)\leq\delta_1^2 &\omega\in[0,\omega_p]\\
+&(1-\delta_1)^2\leq R(\omega)\leq(1+\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]\\
@@ Line 100: / Line 100: @@
 \textrm{minimize}_{G\,,\,r}&\langle G\,,W\rangle\lambda+\|\textrm{diag}(G)\|_\infty\\
 \textrm{subject~to}
-& R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2r(n)\cos(\omega n)\\
+& R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2(r_r(n)\cos(\omega n)+r_i(n)\sin(\omega n))\\
-&\frac{1}{\delta_1^2}\leq R(\omega)\leq\delta_1^2 &\omega\in[0,\omega_p]\\
+&(1-\delta_1)^2\leq R(\omega)\leq(1+\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]\\

Filter design by convex iteration

From Wikimization

Revision as of 07:51, 1 November 2010

Contents

RF pulse design with small flip angle

RF pulse design with large flip angle

A typical design procedure involves the following steps:

Minimum peak RF waveform for saturation pulse

Minimum peak time-domain response by Optimization

Matlab demonstration

Views

Personal tools

Navigation

Search

Toolbox

Minimum peak $LaTeX: B_1(t)$ RF waveform for saturation pulse