# Filter design by convex iteration

### From Wikimization

(→demonstration in Matlab) |
|||

Line 1: | Line 1: | ||

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

+ | |||

+ | 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 that will produce the desired transversed magnetization pulse profile and lognitudinal magnetization pulse profile. SLR transform relates the desired profile to two complex polynomials. 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, <math>B_N(z)</math>, where <math>z=e^{j\gamma\,G_x\,\delta\,t}</math>. | ||

+ | |||

+ | Forward SLR transform relates time domain RF waveform to the frequency domain pulse response. Inverse SLR transform is the converse. | ||

+ | |||

<math>H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}</math> | <math>H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}</math> | ||

where <math>h(n)\in\mathbb{C}^N</math> denotes impulse response. | where <math>h(n)\in\mathbb{C}^N</math> denotes impulse response. |

## Revision as of 12:55, 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 pulse will tip the net magnetization vector to the longitudinal axis. If the desired flip angle is small , 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 that will produce the desired transversed magnetization pulse profile and lognitudinal magnetization pulse profile. SLR transform relates the desired profile to two complex polynomials. 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, , where .

Forward SLR transform relates time domain RF waveform to the frequency domain pulse response. Inverse SLR transform is the converse.

where denotes impulse response.

For a low pass filter, frequency domain specifications are:

To minimize peak magnitude of , the problem becomes

But this problem statement is nonconvex.

So instead,

is a positive semidefinite rank 1 matrix.

Summing along each of subdiagonals gives entries of the autocorrelation function of .

In particular, the main diagonal of holds squared entries of .

Minimizing is therefore equivalent to minimizing .

Define and define as a zero matrix having vector along the superdiagonal when is positive or along the subdiagonal when is negative.

By spectral factorization, , an equivalent problem is:

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