function [rf, rf_g] = design_dualband_pulse(N, pw, bw1, bw2, foff2, d1, d2);
% DESIGN_DUALBAND_PULSE Creates a maximum-phase dual-band RF suppression 
%        pulse using the complex Parks-McClellan filter design algorithm
%    [rf, rf_g] = design_dualband_pulse(N, pw, bw1, bw2, foff2, [d1, d2])
%
%    Uses the complex Parks-McClellan filter design algorithm (see CFIRPM) 
%    to design a dual-band RF pulse.  The input bandwidths are the exact
%    suppression bandwidths,  not the full-width, half-maximum bandwidths.
%    The resulting pulse is maximum phase.
%
%    A filter order of N > 70 is recommended.
% Inputs:
%    N - number of points in pulse / filter order
%    pw - pulsewidth in seconds
%    bw1 - Bandwidth of on-resonance suppression band
%    bw2 - Bandwidth of second suppression band (eg for fat sat)
%    foff2 - frequency offset of second suppression band, in Hz
%    d1 (optional) - ripple in suppression bands (default is .01)
%    d2 (optional) - ripple elsewhere (default is .01)
% Outputs:
%    rf - normalized rf pulse (sums to flip angle)
%    rf_g - rf pulse in Gauss
%
% NOTE: fmp.m, b2rf.m, b2a.m, ab2rf.m and rfscale.m must be installed to 
% run this function.
% They are available at http://www-mrsrl.stanford.edu/~peder/longt2supp/
%
% Written by Peder Larson, 12/13/2005
% (c) Board of Trustees, Leland Stanford Junior University

if (nargin < 6) % no optional arguments supplied
  d1 = .01;
  d2 = .01;
end

% Adjust ripples for suppression pulse
d1 = d1/2;
d2 = sqrt(d2);

% normalized frequencies
bw1n = bw1/2*pw;
bw2n = bw2/2*pw;
foff2n = foff2*pw;

% for Max/Min phase pulse
d1m = 2*d1; d2m = d2^2/2;
di = 0.5*dinf(d1m,d2m);

% transition width
bw_trans = di/pw;
% normalized transition width
bw_t = bw_trans*pw;

if (foff2n < 0)
  f = [-N/2 foff2n-bw2n-bw_t foff2n-bw2n foff2n+bw2n foff2n+bw2n+bw_t ...
       -bw1n-bw_t -bw1n  bw1n  bw1n+bw_t N/2]/(N/2);
else
  f = [-N/2 -bw1n-bw_t -bw1n  bw1n  bw1n+bw_t foff2n-bw2n-bw_t foff2n-bw2n ...
       foff2n+bw2n foff2n+bw2n+bw_t N/2]/(N/2);
end

if any(diff(f) <= 0),
   error(['Pulse with given parameters is not feasible.\n' ...
          'Try decreasing the bandwidths, increasing the pulse duration\n' ...
          'or increasing the ripple values'],1)
end

a = [0    0                1           1           0 ...
     0          1      1   0         0];

% Creates linear phase filter
Wpass = d2m/d1m;
w = [1  Wpass 1 Wpass 1];
h_lin =  cfirpm((N-1)*2,f, a, w);
% Converts to maximum phase filter - this compromises the frequency response
h = fmp(h_lin(end:-1:1)); % h is of length N

rf = b2rf(h*sqrt(1/2));
% scale to pi/2
rf = pi/2*rf/sum(rf);
rf_g = rfscaleg(rf, pw*1e3);