# Filter design by convex iteration

## 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 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^{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 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^{-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} =\!\left[\!\begin{array}{ccccccc} |h(0)|^2&h(0)h(1)^*&h(0)h(2)^*&\cdots&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)^*&\cdots&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&\cdots&h(2)h(N\!-\!3)^*&h(2)h(N\!-\!2)^*&h(2)h(N\!-\!1)^*\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\ h(N\!-\!3)h(0)^*&h(N\!-\!3)h(1)^*&h(N\!-\!3)h(2)^*&\cdots&|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)^*&\cdots&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)^*&\cdots&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\,\triangleq r_{\rm re}+ir_{\rm im}\,\,\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 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}_{G\,,\,r}&\|\textrm{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))\\ &(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\\ & \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_{\rm re}(0)+\!\sum\limits_{n=1}^{N-1}2(r_{\rm re}(n)\cos(\omega n)+r_{\rm im}(n)\sin(\omega n))\\ &(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}$

## 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);
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('\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)';
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 = 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)')


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