Filter design by convex iteration

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Current revision (17:12, 18 September 2016) (edit) (undo)
(Minimum peak time-domain response by Convex Optimization)
 
(48 intermediate revisions not shown.)
Line 1: Line 1:
-
'''RF pulse design with small flip angle'''
+
==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 <math>90^o</math> pulse will tip the net magnetization vector to the longitudinal axis. If the desired flip angle <math>\theta</math> is small <math>(\sin(\theta)\simeq\theta)</math>, then the time domain RF waveform and the pulse response can be approximated by Fourier Transform.
-
Flip angle is the amount of devriation one applies to net magnetization vector from its transverse axis. A <math>90^o</math> pulse will tip the net magnetization vector to the longitudinal axis. If the desired flip angle <math>\theta</math> is small <math>(\sin(\theta)\simeq\theta)</math>, then the time domain RF waveform and the pulse response can be approximated by Fourier Transform.
+
-
 
+
-
 
+
-
'''RF pulse design with large flip angle'''
+
 +
==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.
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 transverse magnetization pulse profile and longitudinal 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>
Line 14: Line 11:
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.
-
 
+
==A typical design procedure involves the following steps:==
-
'''A typical design procedure involves the following steps:'''
+
'''1.'''
-
 
+
-
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.
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.
+
'''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>
+
Use Parks-McClellan algorithm to come up with a waveform, <math>B_N(z)</math>, where <math>z=e^{\jmath\,\gamma\,G_x\,\Delta\,t}</math>
-
3.
+
'''3.'''
<math>A_N(z)</math> and <math>B_N(z)</math> are related by
<math>A_N(z)</math> and <math>B_N(z)</math> are related by
<math>|A_N(z)| = \sqrt{ 1-B_N(z)\,B_N^\ast\,(z)}</math>.
<math>|A_N(z)| = \sqrt{ 1-B_N(z)\,B_N^\ast\,(z)}</math>.
-
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.
+
If <math>A_N(z)</math> 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 <math>B_N(z)</math>, the unique minimum-phase <math>A_N(z)</math> is
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.
+
<math>A_N(z)=|A_N(z)|\,e^{\jmath\,\mathcal{H}(\mathrm{log}|A_N(z)|)}</math>, where <math>\mathcal{H}</math> is the Hilbert transform operator.
-
4.
+
'''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.
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 transverse 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^{\jmath\,\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 have 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.&nbsp; Since the application is for saturation, phase of the profile response is not important. But the passband should be flat and transition region sharp.
-
'''Minimum peak <math>B_1(t)</math> RF waveform for saturation pulse'''
+
Plan: Follow the four steps above except instead of using Parks-McClellan to find <math>B_N(z)</math>, in step 2, use Optimization:
-
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.
+
==Minimum peak time-domain response by Convex Optimization==
-
 
+
<math>H(\omega) = h(0) + h(1)e^{-\jmath\,\omega} + ... + h(N-1)e^{-\jmath(N-1)\omega}</math>
-
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.
+
-
--------------------------------------------
+
-
 
+
-
 
+
-
<math>H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}</math>
+
&nbsp;where&nbsp; <math>h(n)\in\mathbb{C}^N</math> denotes impulse response.
&nbsp;where&nbsp; <math>h(n)\in\mathbb{C}^N</math> denotes impulse response.
Line 55: Line 47:
<math>
<math>
\begin{array}{ll}
\begin{array}{ll}
-
\frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\
+
\frac{1}{1+\delta_1}\leq|H(\omega)|\leq1+\delta_1\,, &\omega\in[0,\omega_p]\\
|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
\end{array}
\end{array}
</math>
</math>
-
To minimize peak magnitude of <math>h\,</math>&nbsp;, the problem becomes
+
To minimize peak amplitude of impulse response <math>h\,</math>&nbsp;, the problem becomes
<math>\begin{array}{cll}
<math>\begin{array}{cll}
\textrm{minimize} &\|h\|_\infty\\
\textrm{minimize} &\|h\|_\infty\\
-
\textrm{subject~to} &\frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\
+
\textrm{subject~to} &\frac{1}{1+\delta_1}\leq|H(\omega)|\leq1+\delta_1\,, &\omega\in[0,\omega_p]\\
&|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
&|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
\end{array}</math>
\end{array}</math>
-
But this problem statement is nonconvex.
+
But this problem statement is nonconvex because <math>\frac{1}{1+\delta_1}\leq|H(\omega)|</math> specifies a nonconvex set of feasible <math>h\,</math>.
So instead,
So instead,
-
<math>G\triangleq h(n)h(n)^{\rm H} \in\mathbb{C}^{N\times N}</math>
+
<math>G=h(n)h(n)^{\rm H}
 +
=\!\left[\!\begin{array}{ccccccc}
 +
|h(0)|^2&h(0)h(1)^*&h(0)h(2)^*&...&h(0)h(N\!-\!3)^*&h(0)h(N\!-\!2)^*&h(0)h(N\!-\!1)^*\\
 +
h(1)h(0)^*&|h(1)|^2&h(1)h(2)^*&...&h(1)h(N\!-\!3)^*&h(1)h(N\!-\!2)^*&h(1)h(N\!-\!1)^*\\
 +
h(2)h(0)^*&h(2)h(1)^*&|h(2)|^2&...&h(2)h(N\!-\!3)^*&h(2)h(N\!-\!2)^*&h(2)h(N\!-\!1)^*\\
 +
:&:&:&.&:&:&:\\
 +
h(N\!-\!3)h(0)^*&h(N\!-\!3)h(1)^*&h(N\!-\!3)h(2)^*&...&|h(N\!-\!3)|^2&h(N\!-\!3)h(N\!-\!2)^*&h(N\!-\!3)h(N\!-\!1)^*\\
 +
h(N\!-\!2)h(0)^*&h(N\!-\!2)h(1)^*&h(N\!-\!2)h(2)^*&...&h(N\!-\!2)h(N\!-\!3)^*&|h(N\!-\!2)|^2&h(N\!-\!2)h(N\!-\!1)^*\\
 +
h(N\!-\!1)h(0)^*&h(N\!-\!1)h(1)^*&h(N\!-\!1)h(2)^*&...&h(N\!-\!1)h(N\!-\!3)^*&h(N\!-\!1)h(N\!-\!2)^*&|h(N\!-\!1)|^2
 +
\end{array}\!\right]\!
 +
\in\mathbb{C}^{N\times N}</math>
is a positive semidefinite rank 1 matrix.
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\,= r_{\rm re}+\jmath\,r_{\rm im}\,\in\mathbb{C}^N</math>&nbsp; and <math>r(n)=r(-n)^*\,</math>.
In particular, the main diagonal of <math>G\,</math> holds squared entries of <math>h\,</math>&nbsp;.
In particular, the main diagonal of <math>G\,</math> holds squared entries of <math>h\,</math>&nbsp;.
-
Minimizing <math>\|h\|_\infty</math> is therefore equivalent to minimizing <math>\|\textrm{diag}(G)\|_\infty</math>&nbsp;.
+
Minimizing <math>\|h\|_\infty</math> is therefore equivalent to minimizing <math>\|{\text diag}(G)\|_\infty</math>&nbsp;.
-
Define <math>I_0\triangleq I\,</math> and define <math>I_n\,</math> as a zero matrix having vector <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> superdiagonal when <math>n\,</math> is positive or <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> subdiagonal when <math>n\,</math> is negative.
+
Define <math>I_0= I\,</math> and define <math>I_n\,</math> as a square zero-matrix having vector <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> superdiagonal when <math>n\,</math> is positive or <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> subdiagonal when <math>n\,</math> is negative.
-
By spectral factorization, <math>R(\omega)=|H(\omega)|^2\,</math>&nbsp;,
+
By spectral factorization, <math>R(\omega)=|H(\omega)|^2\,</math>,&nbsp;
an equivalent problem is:
an equivalent problem is:
<math>\begin{array}{cll}
<math>\begin{array}{cll}
-
\textrm{minimize}_{G\,,\,r}&\|\textrm{diag}(G)\|_\infty\\
+
\textrm{minimize}\limits_{G\,,\,r}&\|{\text diag}(G)\|_\infty\\
\textrm{subject~to}
\textrm{subject~to}
-
& R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2r(n)\cos(\omega n)\\
+
& R(\omega) = r_{\rm re}(0)+\!\sum\limits_{n=1}^{N-1}2(r_{\rm re}(n)\cos(\omega n)+r_{\rm im}(n)\sin(\omega n))\\
-
&\frac{1}{\delta_1^2}\leq R(\omega)\leq\delta_1^2 &\omega\in[0,\omega_p]\\
+
&\frac{1}{(1+\delta_1)^2}\leq R(\omega)\leq(1+\delta_1)^2 &\omega\in[-\omega_p,\omega_p]\\
-
& R(\omega)\leq\delta_2^2 & \omega\in[\omega_s\,,\pi]\\
+
& R(\omega)\leq\delta_2^2 & \pm\omega\in[\omega_s\,,\pi]\\
-
& R(\omega)\geq0 & \omega\in[0,\pi]\\
+
& R(\omega)\geq0 & \omega\in[-\pi,\pi]\\
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
& G\succeq0\\
& G\succeq0\\
-
& \textrm{rank}(G) = 1
+
& {\text rank}(G) = 1
\end{array}</math>
\end{array}</math>
Excepting the rank constraint, this problem is convex.
Excepting the rank constraint, this problem is convex.
-
==demonstration in Matlab==
+
The technique of [[Convex Iteration]] (to find direction vector <math>W</math>) from Dattorro's book is applied:
 +
 
 +
<math>\begin{array}{cll}
 +
\textrm{minimize}\limits_{G\,,\,r}&\langle G\,,W\rangle\\
 +
\textrm{subject~to}
 +
& \max{\text diag}(G)\leq\kappa^2\\
 +
& R(\omega) = r_{\rm re}(0)+\!\sum\limits_{n=1}^{N-1}2(r_{\rm re}(n)\cos(\omega n)+r_{\rm im}(n)\sin(\omega n))\\
 +
&\frac{1}{(1+\delta_1)^2}\leq R(\omega)\leq(1+\delta_1)^2 &\omega\in[-\omega_p,\omega_p]\\
 +
& R(\omega)\leq\delta_2^2 & \pm\omega\in[\omega_s\,,\pi]\\
 +
& R(\omega)\geq0 & \omega\in[-\pi,\pi]\\
 +
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
 +
& G\succeq0
 +
\end{array}</math>
 +
 
 +
where peak impulse response amplitude <math>\kappa\,</math> is a constant chosen by cut-and-try or bisection.
 +
 
 +
==Matlab demonstration==
<pre>
<pre>
% CVX code by Almir Mutapcic in 2006.
% CVX code by Almir Mutapcic in 2006.
-
% Adapted in 2010 for impulse response peak-minimization by convex iteration.
+
% 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
%
%
-
% S.-P. Wu, S. Boyd, and L. Vandenberghe,
+
% Designs an FIR lowpass filter using spectral factorization method with
-
% "FIR Filter Design via Spectral Factorization and Convex Optimization":
+
% constraint on maximum passband ripple and stopband attenuation:
%
%
% 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
%
%
-
% Change variables via spectral factorization method and get:
+
% We 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 119: Line 140:
% R(w) >= 0 for all w
% R(w) >= 0 for all w
%
%
-
% where R(w) is squared magnitude of frequency response
+
% where R(w) is squared magnitude 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.
-
% delta2 is specified stopband attenuation.
+
% This is a convex problem (can be formulated as an SDP after sampling).
-
 
+
clear all, clc, close all, fclose('all');
clear all, clc, close all, fclose('all');
rand('twister',sum(100*clock));
rand('twister',sum(100*clock));
Line 133: Line 153:
% number of FIR coefficients (including zeroth)
% number of FIR coefficients (including zeroth)
N = 32;
N = 32;
-
 
+
wpass = 0.12*pi; % end of passband
-
wpass = 0.12*pi; % end of the passband
+
wstop = 0.20*pi; % start of stopband
-
wstop = 0.20*pi; % start of the stopband
+
delta0_wpass = 0.125;
-
delta = 1; % maximum passband ripple in dB (+/- around 0 dB)
+
delta0_wstop = 0.125;
-
delta2 = -20; % stopband attenuation desired in dB
+
delta = 20*log10(1 + delta0_wpass) % maximum passband ripple in dB (+/- around 0 dB)
 +
delta2 = 20*log10(delta0_wstop) % stopband attenuation desired in dB
%*********************************************************************
%*********************************************************************
Line 151: Line 172:
% passband 0 <= w <= w_pass
% passband 0 <= w <= w_pass
-
ind = find((0 <= w) & (w <= wpass)); % passband
+
ind_p = find((0 <= w) & (w <= wpass)); % passband
-
Lp = 10^(-delta/20)*ones(length(ind),1);
+
Lp = 10^(-delta/20);
-
Up = 10^(+delta/20)*ones(length(ind),1);
+
Up = 10^(+delta/20);
-
Ap = A(ind,:);
+
Ap = A(ind_p,:);
% stopband (w_stop <= w)
% stopband (w_stop <= w)
-
ind = find((wstop <= w) & (w <= pi)); % stopband
+
ind_s = find((wstop <= w) & (w <= pi)); % stopband
-
Sp = 10^(delta2/20)*ones(length(ind),1);
+
Sp = 10^(delta2/20);
-
As = A(ind,:);
+
As = A(ind_s,:);
 +
 
 +
%remove redundant contraints
 +
ind_nr = setdiff(1:m,ind_p); % fullband less passband
 +
Anr = A(ind_nr,:);
% make I matrices
% make I matrices
Line 171: Line 196:
%initial direction vector
%initial direction vector
W = eye(N);
W = eye(N);
- 
-
% peak impulse response
 
-
peak = 0.125;
 
%********************************************************************
%********************************************************************
% optimization
% optimization
%********************************************************************
%********************************************************************
-
% formulate and solve the magnitude design problem
+
convergence = [];
iteration = 1;
iteration = 1;
while 1
while 1
tic
tic
-
fprintf('\nMinimizing impulse response peak: iteration %d\n', iteration);
+
fprintf('\niteration %d\n', iteration);
cvx_quiet(true);
cvx_quiet(true);
-
% cvx_solver('sdpt3');
+
cvx_solver('sdpt3');
 +
cvx_precision('best');
cvx_begin
cvx_begin
variable r(N,1);
variable r(N,1);
variable G(N,N) symmetric;
variable G(N,N) symmetric;
minimize(trace(W'*G));
minimize(trace(W'*G));
 +
%peak impulse response by cut and try
 +
max(diag(G)) <= 0.1053^2;
% passband constraints
% passband constraints
Ap*r >= Lp.^2;
Ap*r >= Lp.^2;
Line 195: Line 220:
As*r <= Sp.^2;
As*r <= Sp.^2;
% nonnegative-real constraint
% nonnegative-real constraint
-
A*r >= 0;
+
Anr*r >= 0;
% relate r to h
% relate r to h
r == B*vect(G);
r == B*vect(G);
G == semidefinite(N);
G == semidefinite(N);
-
%minimize peak of h
 
-
diag(G) <= peak^2;
 
cvx_end
cvx_end
toc
toc
-
%computer new direction vector
+
%compute new direction vector
[v,d,q] = svd(G);
[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));
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)*v(:,2:N)';
-
% W = v(:,2:N)*diag(rand(N-1,1))*v(:,2:N)';
 
rankG = sum(diag(d) > max(diag(d))*1e-5);
rankG = sum(diag(d) > max(diag(d))*1e-5);
fprintf('rank(G)=%d, trace(W*G)=%f\n', rankG, trace(G*W));
fprintf('rank(G)=%d, trace(W*G)=%f\n', rankG, trace(G*W));
-
if (rankG == 1)
 
-
break
 
-
end
 
figure(1)
figure(1)
% FIR impulse response
% FIR impulse response
-
h = G(:,1)/sqrt(G(1,1));
+
h = v(:,1)*sqrt(d(1,1));
-
plot([0:N-1],h','o',[0:N-1],h','b:')
+
plot([0:N-1],h','ob:')
xlabel('t'), ylabel('h(t)')
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)
pause(1)
% check if problem was successfully solved
% check if problem was successfully solved
disp(['Problem is ' cvx_status])
disp(['Problem is ' cvx_status])
 +
if (rankG == 1)
 +
break
 +
end
if ~strfind(cvx_status,'Solved')
if ~strfind(cvx_status,'Solved')
-
return
+
fprintf(2,'Excuse me.\n')
end
end
Line 239: Line 269:
%*********************************************************************
%*********************************************************************
% frequency response of the designed filter, where j = sqrt(-1)
% 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;
H = [exp(-j*kron(w,[0:N-1]))]*h;
-
figure(2)
+
figure(2);
 +
subplot(121)
% magnitude
% magnitude
plot(w,20*log10(abs(H)), ...
plot(w,20*log10(abs(H)), ...
Line 251: Line 281:
ylabel('mag H(w) in dB')
ylabel('mag H(w) in dB')
axis([0 pi -50 5])
axis([0 pi -50 5])
 +
title(sprintf('N=%d, w_p(pi)=%3.2f, w_s(pi)=%3.2f, delta=%3.2f', N, wpass/pi, wstop/pi, delta));
%compare impulse response designed by conventional method
%compare impulse response designed by conventional method
-
figure(3)
+
subplot(122)
-
h = spectral_fact(r); %from CVX distribution, Examples subdirectory
+
h_sp = spectral_fact(r); %from CVX distribution, Examples subdirectory
-
plot([0:N-1],h','o',[0:N-1],h','b:')
+
plot([0:N-1],h_sp','+r--');
-
xlabel('t'), ylabel('h(t)')
+
hold on;
-
 
+
plot([0:N-1],h(end:-1:1)','ob:');
 +
legend('conventional','optimal');
 +
xlabel('t'), ylabel('h(t)'); grid
 +
title(sprintf('h_{max} conventional=%3.4f, h_{max} optimal=%3.4f',max(abs(h_sp)),max(abs(h))));
 +
set(gcf,'Outerposition',[300 300 256*4 256*2])
 +
figure(1)
figure(1)
% FIR impulse response
% FIR impulse response
-
h = G(:,1)/sqrt(G(1,1));
+
plot([0:N-1],h','ob:');
-
plot([0:N-1],h','o',[0:N-1],h','b:')
+
xlabel('t'), ylabel('h(t)')
xlabel('t'), ylabel('h(t)')
 +
</pre>
 +
 +
=== spectral_fact() ===
 +
<pre>
 +
% Spectral factorization using Kolmogorov 1939 approach.
 +
% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)
 +
%
 +
% Computes the minimum-phase impulse response which satisfies
 +
% given auto-correlation.
 +
%
 +
% Input:
 +
% r: top-half of the auto-correlation coefficients
 +
% starts from 0th element to end of the auto-corelation
 +
% should be passed in as a column vector
 +
% Output
 +
% h: impulse response that gives the desired auto-correlation
 +
 +
function h = spectral_fact(r)
 +
 +
% length of the impulse response sequence
 +
n = length(r);
 +
 +
% over-sampling factor
 +
mult_factor = 100; % should have mult_factor*(n) >> n
 +
m = mult_factor*n;
 +
 +
% computation method:
 +
% H(exp(jTw)) = alpha(w) + j*phi(w)
 +
% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))
 +
 +
% compute 1/2*ln(R(w))
 +
w = 2*pi*[0:m-1]/m;
 +
R = [ ones(m,1) 2*cos(kron(w',[1:n-1])) ]*r;
 +
alpha = 1/2*log(R);
 +
 +
% find the Hilbert transform
 +
alphatmp = fft(alpha);
 +
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
 +
alphatmp(1) = 0;
 +
alphatmp(floor(m/2)+1) = 0;
 +
phi = real(ifft(j*alphatmp));
 +
 +
% now retrieve the original sampling
 +
index = find(rem([0:m-1],mult_factor)==0);
 +
alpha1 = alpha(index);
 +
phi1 = phi(index);
 +
 +
% compute the impulse response (inverse Fourier transform)
 +
h = real(ifft(exp(alpha1+j*phi1),n));
</pre>
</pre>

Current revision

Contents

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^{\jmath\,\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^{\jmath\,\mathcal{H}(\mathrm{log}|A_N(z)|)}, where LaTeX: \mathcal{H} is the Hilbert transform operator.

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^{\jmath\,\omega} and LaTeX: M_o is the initial magnetization. If LaTeX: |B_N\,(\omega)| represents a rectangular profile, LaTeX: |M_{xy}(\omega)| will also have 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, phase of the profile response is not important. But the passband should be flat and transition region sharp.

Plan: Follow the four steps above except instead of using Parks-McClellan to find LaTeX: B_N(z), in step 2, use Optimization:

Minimum peak time-domain response by Convex Optimization

LaTeX: H(\omega) = h(0) + h(1)e^{-\jmath\,\omega} + ... + h(N-1)e^{-\jmath(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}{1+\delta_1}\leq|H(\omega)|\leq1+\delta_1\,, &\omega\in[0,\omega_p]\\
|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi]
\end{array}

To minimize peak amplitude of impulse response LaTeX: h\, , the problem becomes

LaTeX: \begin{array}{cll}
\textrm{minimize} &\|h\|_\infty\\
\textrm{subject~to} &\frac{1}{1+\delta_1}\leq|H(\omega)|\leq1+\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 because LaTeX: \frac{1}{1+\delta_1}\leq|H(\omega)| specifies a nonconvex set of feasible LaTeX: h\,.

So instead,

LaTeX: G=h(n)h(n)^{\rm H} 
=\!\left[\!\begin{array}{ccccccc}
|h(0)|^2&h(0)h(1)^*&h(0)h(2)^*&...&h(0)h(N\!-\!3)^*&h(0)h(N\!-\!2)^*&h(0)h(N\!-\!1)^*\\
h(1)h(0)^*&|h(1)|^2&h(1)h(2)^*&...&h(1)h(N\!-\!3)^*&h(1)h(N\!-\!2)^*&h(1)h(N\!-\!1)^*\\
h(2)h(0)^*&h(2)h(1)^*&|h(2)|^2&...&h(2)h(N\!-\!3)^*&h(2)h(N\!-\!2)^*&h(2)h(N\!-\!1)^*\\
</p>
<dl><dd>&:&:&.&:&:&:\\
</dd></dl>
<p>h(N\!-\!3)h(0)^*&h(N\!-\!3)h(1)^*&h(N\!-\!3)h(2)^*&...&|h(N\!-\!3)|^2&h(N\!-\!3)h(N\!-\!2)^*&h(N\!-\!3)h(N\!-\!1)^*\\
h(N\!-\!2)h(0)^*&h(N\!-\!2)h(1)^*&h(N\!-\!2)h(2)^*&...&h(N\!-\!2)h(N\!-\!3)^*&|h(N\!-\!2)|^2&h(N\!-\!2)h(N\!-\!1)^*\\
h(N\!-\!1)h(0)^*&h(N\!-\!1)h(1)^*&h(N\!-\!1)h(2)^*&...&h(N\!-\!1)h(N\!-\!3)^*&h(N\!-\!1)h(N\!-\!2)^*&|h(N\!-\!1)|^2
\end{array}\!\right]\!
\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\,= r_{\rm re}+\jmath\,r_{\rm im}\,\in\mathbb{C}^N  and LaTeX: r(n)=r(-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: \|{\text diag}(G)\|_\infty .


Define LaTeX: I_0= I\, and define LaTeX: I_n\, as a square 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}\limits_{G\,,\,r}&\|{\text diag}(G)\|_\infty\\
\textrm{subject~to} 
& R(\omega) = r_{\rm re}(0)+\!\sum\limits_{n=1}^{N-1}2(r_{\rm re}(n)\cos(\omega n)+r_{\rm im}(n)\sin(\omega n))\\
&\frac{1}{(1+\delta_1)^2}\leq R(\omega)\leq(1+\delta_1)^2 &\omega\in[-\omega_p,\omega_p]\\
& R(\omega)\leq\delta_2^2 & \pm\omega\in[\omega_s\,,\pi]\\
& R(\omega)\geq0 & \omega\in[-\pi,\pi]\\
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
& G\succeq0\\
& {\text 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:

LaTeX: \begin{array}{cll}
\textrm{minimize}\limits_{G\,,\,r}&\langle G\,,W\rangle\\
\textrm{subject~to}
& \max{\text diag}(G)\leq\kappa^2\\
& R(\omega) = r_{\rm re}(0)+\!\sum\limits_{n=1}^{N-1}2(r_{\rm re}(n)\cos(\omega n)+r_{\rm im}(n)\sin(\omega n))\\
&\frac{1}{(1+\delta_1)^2}\leq R(\omega)\leq(1+\delta_1)^2 &\omega\in[-\omega_p,\omega_p]\\
& R(\omega)\leq\delta_2^2 & \pm\omega\in[\omega_s\,,\pi]\\
& R(\omega)\geq0 & \omega\in[-\pi,\pi]\\
& r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\
& G\succeq0
\end{array}

where peak impulse response amplitude LaTeX: \kappa\, is a constant chosen by cut-and-try or bisection.

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;
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);
Up  = 10^(+delta/20);
Ap  = A(ind_p,:);

% stopband (w_stop <= w)
ind_s = find((wstop <= w) & (w <= pi));   % stopband
Sp  = 10^(delta2/20);
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('\niteration %d\n', iteration);
   cvx_quiet(true);
   cvx_solver('sdpt3');
   cvx_precision('best');
   cvx_begin
      variable r(N,1);
      variable G(N,N) symmetric;
      minimize(trace(W'*G));
      %peak impulse response by cut and try
      max(diag(G)) <= 0.1053^2;
      % 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)';
   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 = v(:,1)*sqrt(d(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 = [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=%d, w_p(pi)=%3.2f, w_s(pi)=%3.2f, delta=%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','+r--');
hold on;
plot([0:N-1],h(end:-1:1)','ob:');
legend('conventional','optimal');
xlabel('t'), ylabel('h(t)'); grid
title(sprintf('h_{max} conventional=%3.4f, h_{max} optimal=%3.4f',max(abs(h_sp)),max(abs(h))));
set(gcf,'Outerposition',[300 300 256*4 256*2])
 
figure(1)
% FIR impulse response
plot([0:N-1],h','ob:'); 
xlabel('t'), ylabel('h(t)')

spectral_fact()

% Spectral factorization using Kolmogorov 1939 approach.
% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)
%
% Computes the minimum-phase impulse response which satisfies
% given auto-correlation.
%
% Input:
%   r: top-half of the auto-correlation coefficients
%      starts from 0th element to end of the auto-corelation
%      should be passed in as a column vector
% Output
%   h: impulse response that gives the desired auto-correlation

function h = spectral_fact(r)

% length of the impulse response sequence
n = length(r);

% over-sampling factor
mult_factor = 100;        % should have mult_factor*(n) >> n 
m = mult_factor*n;

% computation method:
% H(exp(jTw)) = alpha(w) + j*phi(w)
% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))

% compute 1/2*ln(R(w))
w = 2*pi*[0:m-1]/m;
R = [ ones(m,1) 2*cos(kron(w',[1:n-1])) ]*r;
alpha = 1/2*log(R);

% find the Hilbert transform 
alphatmp = fft(alpha);
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
alphatmp(1) = 0;
alphatmp(floor(m/2)+1) = 0;
phi = real(ifft(j*alphatmp));

% now retrieve the original sampling 
index  = find(rem([0:m-1],mult_factor)==0);
alpha1 = alpha(index);
phi1   = phi(index);

% compute the impulse response (inverse Fourier transform)
h = real(ifft(exp(alpha1+j*phi1),n));
Personal tools