Filter design by convex iteration

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Line 8: Line 8:
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.
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 <math> B_1(t)</math> that will produce the desired transversed magnetization pulse profile and lognitudinal magnetization pulse profile. SLR transform relates the RF waveform to two complex polynomials <math>A_N\,(z)</math> and <math>B_N\,(z)</math>.
+
The goal of RF pulse design is find the time domain waveform <math> B_1(t)</math> that will produce the desired transversed magnetization pulse profile and lognitudinal magnetization pulse profile. SLR transform relates the RF waveform to two complex polynomials <math>A_N(z)</math> and <math>B_N(z)</math>.
-
''i.e.'' <math> B_1(t)\begin{array}{c}\mbox{SLR}\\\Leftrightarrow \end{array}\, A_N\,(z),\, B_N\,(z)</math>
+
''i.e.'' <math> B_1(t)\begin{array}{c}\mbox{SLR}\\\Leftrightarrow \end{array}\, A_N(z),\, B_N(z)</math>
Thus, the task of filter design becomes finding these two polynomials, which can be done via FIR filter design methods.
Thus, the task of filter design becomes finding these two polynomials, which can be done via FIR filter design methods.
Line 17: Line 17:
'''A typical design procedure involves the following steps:'''
'''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.
+
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, <math>B_N(z)</math>, where <math>z=e^{j\gamma\,G_x\,\Delta\,t}</math>
+
2.
 +
Use Parks-McClellan algorithm to come up with a waveform, <math>B_N(z)</math>, where <math>z=e^{j\gamma\,G_x\,\Delta\,t}</math>
-
3. <math>A_N\,(z)</math> and <math>B_N\,(z)</math> are related by
+
3.
-
<math>|A_N\,(z)| = \sqrt{ 1-B_N\,(z)\,B_N^\ast\,(z)}</math>.
+
<math>A_N(z)</math> and <math>B_N(z)</math> are related by
-
If <math>A_N\,(z)</math> is choosen 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.
+
<math>|A_N(z)| = \sqrt{ 1-B_N(z)\,B_N^\ast\,(z)}</math>.
-
For the given <math>B_N\,(z)</math>, the unique minimum-phase <math>A_N\,(z)</math> is
+
If <math>A_N(z)</math> is choosen 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.
-
<math>A_N\,(z)=|A_N\,(z)|\,e^{i\mathcal{H}\displaystyle(\mbox{log}|A_N\,(z)|\displaystyle)}</math>, where <math>\mathcal{H}</math> is the Hilbert transform opetator.
+
For the given <math>B_N(z)</math>, the unique minimum-phase <math>A_N(z)</math> is
 +
<math>A_N(z)=|A_N(z)|\,e^{i\mathcal{H}\displaystyle(\mbox{log}|A_N(z)|\displaystyle)}</math>, where <math>\mathcal{H}</math> is the Hilbert transform opetator.
 +
 
 +
4.
 +
Once <math>A_N\,(z)</math> and <math>B_N\,(z)</math> are found, RF waveform <math> B_1(t)</math> can be computed by inverse SLR transform.
-
4. Once <math>A_N\,(z)</math> and <math>B_N\,(z)</math> are found, RF waveform <math> B_1(t)</math> can be computed by inverse SLR transform.
 
Note: upon applying this RF pulse, the resulting transerve magnetization is
Note: upon applying this RF pulse, the resulting transerve magnetization is
-
<math>|M_{xy}(\omega)| = 2\sqrt{1-B_N\,(\omega)|^2}|B_N\,(\omega)|\,M_o</math>
+
<math>|M_{xy}(\omega)| = 2\sqrt{1-B_N\,(\omega)|^2}\,|B_N\,(\omega)|\,M_o</math>
where <math> z=e^{j\omega}</math> and <math>M_o</math> is the initial magnetization.
where <math> z=e^{j\omega}</math> and <math>M_o</math> is the initial magnetization.
If <math>|B_N\,(\omega)|</math> represents a rectangular profile, <math>|M_{xy}(\omega)|</math> will also has a rectangular profile.
If <math>|B_N\,(\omega)|</math> represents a rectangular profile, <math>|M_{xy}(\omega)|</math> will also has a rectangular profile.
 +
 +
 +
'''Minimum peak <math>B_1(t)</math> RF waveform for saturation pulse'''
 +
 +
Our goal here is to find a minimum peak magnitude <math>B_1(t)</math> 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 <math>B_N(z)</math>, use optimization to

Revision as of 14:10, 22 October 2010

RF pulse design with small flip angle

Flip angle is the amount of devriation 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 transversed magnetization pulse profile and lognitudinal 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 choosen 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}\displaystyle(\mbox{log}|A_N(z)|\displaystyle)}, 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 transerve 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 to


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 code 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