# Filter design by convex iteration

(Difference between revisions)
 Revision as of 12:29, 24 August 2010 (edit)← Previous diff Current revision (17:12, 18 September 2016) (edit) (undo) (→Minimum peak time-domain response by Convex Optimization) (59 intermediate revisions not shown.) Line 1: Line 1: - $H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}$ + ==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 $90^o$ pulse will tip the net magnetization vector to the longitudinal axis. If the desired flip angle $\theta$ is small $(\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 $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 $A_N(z)$ and $B_N(z)$. + + ''i.e.'' $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, $B_N(z)$, where $z=e^{\jmath\,\gamma\,G_x\,\Delta\,t}$ + + '''3.''' + $A_N(z)$ and $B_N(z)$ are related by + $|A_N(z)| = \sqrt{ 1-B_N(z)\,B_N^\ast\,(z)}$. + If $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 $B_N(z)$, the unique minimum-phase $A_N(z)$ is + $A_N(z)=|A_N(z)|\,e^{\jmath\,\mathcal{H}(\mathrm{log}|A_N(z)|)}$, where $\mathcal{H}$ is the Hilbert transform operator. + + '''4.''' + Once $A_N\,(z)$ and $B_N\,(z)$ are found, RF waveform $B_1(t)$ can be computed by inverse SLR transform. + + + Note: upon applying this RF pulse, the resulting transverse magnetization is + $|M_{xy}(\omega)| = 2\sqrt{1-|B_N\,(\omega)|^2}\,|B_N\,(\omega)|\,M_o$ + where $z=e^{\jmath\,\omega}$ and $M_o$ is the initial magnetization. + If $|B_N\,(\omega)|$ represents a rectangular profile, $|M_{xy}(\omega)|$ will also have a rectangular profile. + + ==Minimum peak $B_1(t)$ RF waveform for saturation pulse== + Our goal here is to find a minimum peak magnitude $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 $B_N(z)$, in step 2, use Optimization: + + ==Minimum peak time-domain response by Convex Optimization== + $H(\omega) = h(0) + h(1)e^{-\jmath\,\omega} + ... + h(N-1)e^{-\jmath(N-1)\omega}$  where  $h(n)\in\mathbb{C}^N$ denotes impulse response.  where  $h(n)\in\mathbb{C}^N$ denotes impulse response. Line 6: Line 47: $[itex] \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}$ [/itex] - To minimize peak magnitude of $h\,$ , the problem becomes + To minimize peak amplitude of impulse response $h\,$ , the problem becomes $\begin{array}{cll} [itex]\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}$ \end{array}[/itex] - But this problem statement is nonconvex. + But this problem statement is nonconvex because $\frac{1}{1+\delta_1}\leq|H(\omega)|$ specifies a nonconvex set of feasible $h\,$. - So instead, a new vector + So instead, - $+ [itex]G=h(n)h(n)^{\rm H} - g \triangleq \left[ + =\!\left[\!\begin{array}{ccccccc} - \begin{array}{c} + |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(n) \\ + 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(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)^*\\ - \vdots \\ + :&:&:&.&:&:&:\\ - h(n-N) \\ + 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)^*\\ - \end{array} + 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)^*\\ - \right]\in\mathbb{C}^{N^2} + 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}[/itex] - is defined by concatenation of time-shifted versions of $h\,$ . + - + - Then + - + - $G\triangleq gg^\mathrm{H}=\,\left[\begin{array}{*{20}c} + - h(n)h(n)^{\rm H} & h(n)h(n-1)^{\rm H} & h(n)h(n-2)^{\rm H} & \cdots & h(n)h(n-N)^{\rm H}\\ + - h(n-1)h(n)^{\rm H} & h(n-1)h(n-1)^{\rm H} & h(n-1)h(n-2)^{\rm H} & \cdots & h(n-1)h(n-N)^{\rm H}\\ + - h(n-2)h(n)^{\rm H} & h(n-2)h(n-1)^{\rm H} & h(n-2)h(n-2)^{\rm H} & \ddots & h(n-2)h(n-N)^{\rm H}\\ + - \vdots & \vdots & \ddots & \ddots & \vdots\\ + - h(n-N)h(n)^{\rm H} & h(n-N)h(n-1)^{\rm H} & h(n-N)h(n-2)^{\rm H} & \cdots & h(n-N)h(n-N)^{\rm H} + - \end{array}\right]\in\mathbb{C}^{N^2\times N^2}$ + is a positive semidefinite rank 1 matrix. is a positive semidefinite rank 1 matrix. - Summing along each of $2N-1\,$ subdiagonals gives entries of the autocorrelation function $r\,$ of $h\,$ . + Summing along each of $2N-1\,$ subdiagonals gives entries of the autocorrelation function $r\,$ of $h\,$ where $r\,= r_{\rm re}+\jmath\,r_{\rm im}\,\in\mathbb{C}^N$  and $r(n)=r(-n)^*\,$. In particular, the main diagonal of $G\,$ holds squared entries of $h\,$ . In particular, the main diagonal of $G\,$ holds squared entries of $h\,$ . - Minimizing $\|h\|_\infty$ is therefore equivalent to minimizing $\|\textrm{diag}(G)\|_\infty$ . + Minimizing $\|h\|_\infty$ is therefore equivalent to minimizing $\|{\text diag}(G)\|_\infty$ . - Define $I_n\triangleq\,$... + Define $I_0= I\,$ and define $I_n\,$ as a square zero-matrix having vector $\mathbf{1}$ along the $n^{\rm{th}}\,$ superdiagonal when $n\,$ is positive or $\mathbf{1}$ along the $n^{\rm{th}}\,$ subdiagonal when $n\,$ is negative. - By spectral factorization, $R(\omega)=|H(\omega)|^2\,$ , + By spectral factorization, $R(\omega)=|H(\omega)|^2\,$,  an equivalent problem is: an equivalent problem is: $\begin{array}{cll} [itex]\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)=\frac{1}{n}\textrm{trace}(I_{n-N}G) &n=1\ldots N\\ + & r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\ - & r(n)=\frac{1}{n-N}\textrm{trace}(I_{n-N}G) &n=N+1\ldots 2N-1\\ + & G\succeq0\\ & G\succeq0\\ - & \textrm{rank}(G) = 1 + & {\text rank}(G) = 1 \end{array}$ \end{array}[/itex] Excepting the rank constraint, this problem is convex. Excepting the rank constraint, this problem is convex. + The technique of [[Convex Iteration]] (to find direction vector $W$) from Dattorro's book is applied: + + $\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 $\kappa\,$ is a constant chosen by cut-and-try or bisection. + + ==Matlab demonstration==

-                                                           N = 16;                                                                                                                                           +                                                                                                                      % CVX code by Almir Mutapcic in 2006.
-                                                           delta1 = 1.01;                                                                                                                                    +                                                                                                                      % Adapted in 2010 for impulse response peak-minimization by convex iteration by Christine Law.
-                                                           delta2 = 0.01;                                                                                                                                    +                                                                                                                      %
-                                                           w = linspace(0,pi,N);                                                                                                                             +                                                                                                                      % "FIR Filter Design via Spectral Factorization and Convex Optimization"
-                                                           DFT = real(fft(eye(2*N)));  % generate DFT matrix but only use the real components                                                                +                                                                                                                      % by S.-P. Wu, S. Boyd, and L. Vandenberghe
-                                                           % https://ccrma.stanford.edu/~jos/st/DFT_Matrix.html                                                                                              +                                                                                                                      %
+                                                                                                                      % 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

-                                                           cvx_begin                                                                                                                                         +                                                                                                                      %*********************************************************************
-                                                           variable r(2*N-1,1);                                                                                                                              +                                                                                                                      % optimization parameters
-                                                           variable G(N^2,N^2);                                                                                                                              +                                                                                                                      %*********************************************************************
+                                                                                                                      % rule-of-thumb discretization (from Cheney's Approximation Theory)
+                                                                                                                      m = 15*N;
+                                                                                                                      w = linspace(0,pi,m)';  % omega

-                                                           r = [0; r];   % otherwise r is not even                                                                                                           +                                                                                                                      % 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]))];
-                                                           minimize(max(diag(G)));                                                                                                                           +

-                                                           G2 = G(1:N,1:N);       % I don't really need to use the entire G.                                                                                 +                                                                                                                      % passband 0 <= w <= w_pass
-                                                           % So just use the 1st block: G2                                                                                                                   +                                                                                                                      ind_p = find((0 <= w) & (w <= wpass));    % passband
-                                                           for I=1:N                                                                                                                                         +                                                                                                                      Lp  = 10^(-delta/20);
-                                                           temp = diag(ones(1,I),-(N-I));                                                                                                                    +                                                                                                                      Up  = 10^(+delta/20);
-                                                           [i,j,k] = find(temp);                                                                                                                             +                                                                                                                      Ap  = A(ind_p,:);
-                                                           index = [i j ones(length(i),1)];                                                                                                                  +
-                                                           index = [index; [N, N, 0]];                                                                                                                       +                                                                                                                      % stopband (w_stop <= w)
-                                                           temp2 = spconvert(index);                                                                                                                         +                                                                                                                      ind_s = find((wstop <= w) & (w <= pi));   % stopband
-                                                           r(I) == trace(temp2*G2);                                                                                                                          +                                                                                                                      Sp  = 10^(delta2/20);
-                                                           end                                                                                                                                               +                                                                                                                      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

-                                                           for I=N+1:2*N-1                                                                                                                                   +                                                                                                                      %compute new direction vector
-                                                           temp = diag(ones(1,2*N-I),I-N);                                                                                                                   +                                                                                                                      [v,d,q] = svd(G);
-                                                           [i,j,k] = find(temp);                                                                                                                             +                                                                                                                      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));
-                                                           index = [i j ones(length(i),1)];                                                                                                                  +                                                                                                                      W = v(:,2:N)*v(:,2:N)';
-                                                           index = [index; [N, N, 0]];                                                                                                                       +                                                                                                                      rankG = sum(diag(d) > max(diag(d))*1e-5);
-                                                           temp2 = spconvert(index);                                                                                                                         +                                                                                                                      fprintf('rank(G)=%d,  trace(W*G)=%f\n', rankG, trace(G*W));
-                                                           r(I) == trace(temp2*G2);                                                                                                                          +
-                                                           end                                                                                                                                               +                                                                                                                      figure(1)
-                                                                                                                                                                                                             +                                                                                                                      % FIR impulse response
-                                                           R = fftshift(DFT*r);                                                                                                                              +                                                                                                                      h = v(:,1)*sqrt(d(1,1));
-                                                                                                                                                                                                             +                                                                                                                      plot([0:N-1],h','ob:')
-                                                           for I=1:N                 %%% R is of length of 2N                                                                                                +                                                                                                                      xlabel('t'), ylabel('h(t)')
-                                                           1/delta1^2 < R(I);                                                                                                                                +
-                                                           R(I) < delta1^2;                                                                                                                                  +
-                                                           end                                                                                                                                               +

-                                                           for I=N+1:2*N                                                                                                                                     +                                                                                                                      % monitor convergence to 0
-                                                           R(I) < delta2^2;                                                                                                                                  +                                                                                                                      convergence = [convergence trace(G*W)];
-                                                           end                                                                                                                                               +                                                                                                                      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

-                                                           R > 0;                                                                                                                                            +                                                                                                                      iteration = iteration + 1;
-                                                                                                                                                                                                             +                                                                                                                      end
-                                                                                                                                                                                                             +
-                                                           cvx_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); - But the above doesn't work. Debug: try the following + %********************************************************************* + % 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() ===

-                                                           N = 16;                                                                                                                                           +                                                                                                                      % Spectral factorization using Kolmogorov 1939 approach.
-                                                           delta1 = 1.01;                                                                                                                                    +                                                                                                                      % (code follows pp. 232-233, Signal Analysis, by A. Papoulis)
-                                                           delta2 = 0.01;                                                                                                                                    +                                                                                                                      %
-                                                           w = linspace(0,pi,N);                                                                                                                             +                                                                                                                      % Computes the minimum-phase impulse response which satisfies
-                                                           DFT = real(fft(eye(2*N)));                                                                                                                        +                                                                                                                      % 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

-                                                           cvx_begin                                                                                                                                         +                                                                                                                      function h = spectral_fact(r)
-                                                           variable r(2*N-1,1);                                                                                                                              +

-                                                           r = [0; r];                                                                                                                                       +                                                                                                                      % length of the impulse response sequence
+                                                                                                                      n = length(r);

-                                                           R = fftshift(DFT*r);                                                                                                                              +                                                                                                                      % over-sampling factor
-                                                                                                                                                                                                             +                                                                                                                      mult_factor = 100;        % should have mult_factor*(n) >> n
-                                                           for I=1:N                                                                                                                                         +                                                                                                                      m = mult_factor*n;
-                                                           1/delta1^2 <= R(I);                                                                                                                               +
-                                                           R(I) <= delta1^2;                                                                                                                                 +
-                                                           end                                                                                                                                               +
-                                                                                                                                                                                                             +
-                                                           for I=N+1:2*N                                                                                                                                     +
-                                                           R(I) <= delta2^2;                                                                                                                                 +
-                                                           end                                                                                                                                               +

-                                                           R > 0;                                                                                                                                            +                                                                                                                      % computation method:
-                                                                                                                                                                                                             +                                                                                                                      % H(exp(jTw)) = alpha(w) + j*phi(w)
-                                                           cvx_end                                                                                                                                           +                                                                                                                      % where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))
-
+ - The above just to see if there exists an [itex] r such that the problem is feasible. + % compute 1/2*ln(R(w)) - But it still doesn't work... + 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)); +

## 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\,$. $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)^*\\

&:&:&.&:&:&:\\

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