# Filter design by convex iteration

### From Wikimization

(→Matlab demonstration) |
|||

Line 267: | Line 267: | ||

%********************************************************************* | %********************************************************************* | ||

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

Line 285: | Line 284: | ||

subplot(122), | subplot(122), | ||

h_sp = spectral_fact(r); %from CVX distribution, Examples subdirectory | h_sp = spectral_fact(r); %from CVX distribution, Examples subdirectory | ||

- | plot([0:N-1],h_sp',' | + | plot([0:N-1],h_sp','+r--'); |

hold on; | hold on; | ||

h = G(:,1)/sqrt(G(1,1)); | h = G(:,1)/sqrt(G(1,1)); | ||

- | plot([0:N-1],-h(end:-1:1)',' | + | plot([0:N-1],-h(end:-1:1)','ob:'); |

legend('conventional','optimal'); | legend('conventional','optimal'); | ||

xlabel('t'), ylabel('h(t)'); grid | xlabel('t'), ylabel('h(t)'); grid |

## Revision as of 01:31, 10 February 2011

## 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 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 transverse magnetization pulse profile and longitudinal magnetization pulse profile. SLR transform relates the RF waveform to two complex polynomials and .

*i.e.*

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

**3.**
and are related by
.
If 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 , the unique minimum-phase is
, where is the Hilbert transform operator.

**4.**
Once and are found, RF waveform can be computed by inverse SLR transform.

Note: upon applying this RF pulse, the resulting transverse magnetization is
where and is the initial magnetization.
If represents a rectangular profile, will also have a rectangular profile.

## Minimum peak RF waveform for saturation pulse

Our goal here is to find a minimum peak magnitude 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 , in step 2, use Optimization:

## Minimum peak time-domain response by Convex Optimization

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 , where .

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

Minimizing is therefore equivalent to minimizing .

Define and define as a square 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. The technique of Convex Iteration (to find direction vector ) from Dattorro's book is applied. Regularization weight is chosen by cut-and-try as in the figure.

## 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; h = G(:,1)/sqrt(G(1,1)); 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 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));