Wave physics, nonlinear propagation and harmonic imaging

%% m15_wave_fdtd_interface.m
% Purpose : 1D linear acoustic FDTD (staggered pressure / particle velocity, leapfrog)
%           launched at an impedance step and through a quarter-wave layer.
%           Measured reflected / transmitted ENERGY fractions are compared with
%           R = (Z2-Z1)/(Z2+Z1) and 1 - R^2 (band-averaged for the layer case).
%           A second cell computes oblique-incidence energy partition at a
%           soft-tissue / skull-bone interface, including mode conversion to shear.
% Change  : the cases table (c, rho of each medium, layer), f0, FBW, PPW.
% Expect  : fat->muscle  R^2 = 0.0086 (R = 0.093, -20.7 dB)
%           tissue->skull R^2 = 0.288 ; tissue->air-like R^2 = 0.997
%           PZT (30 MRayl)->tissue R^2 = 0.81 bare, 0.18 band-averaged with a lambda/4 layer of
%           sqrt(Z1 Z2) (zero at f0 only; 60 % pulse)
%           skull: critical angle 33.1 deg, > 90 % of the energy enters as shear at 40-70 deg
% Runtime : about 3 s in Octave.
clear; close all; rng(1);

%% Parameters
f0  = 3e6;          % centre frequency (Hz)
FBW = 0.6;          % -6 dB fractional bandwidth of the Gaussian pulse
PPW = 24;           % points per wavelength at f0 in the slowest medium
sf  = FBW*f0/(2*sqrt(2*log(2)));   % spectral std (Hz)
st  = 1/(2*pi*sf);                 % envelope std (s)
pulse = @(t) exp(-t.^2/(2*st^2)).*sin(2*pi*f0*t);

% name, c1, rho1, c2, rho2, layer thickness (m), c_layer, rho_layer
cases = { 'fat -> muscle',            1450,  950, 1580, 1050, 0, 0, 0;
          'soft tissue -> skull',     1540, 1050, 2820, 1900, 0, 0, 0;
          'soft tissue -> air-like',  1540, 1050, 1000,  1.2, 0, 0, 0;
          'PZT -> tissue, bare',      4000, 7500, 1540, 1050, 0, 0, 0;
          'PZT -> tissue, lambda/4',  4000, 7500, 1540, 1050, 2500/(4*f0), 2500, NaN };
Z1Z2 = 4000*7500*1540*1050; cases{5,8} = sqrt(Z1Z2)/2500;   % layer impedance sqrt(Z1 Z2)

fprintf('%-26s %8s %8s %8s %8s %9s %9s\n', 'case', 'R th', 'R meas', 'T th', 'T meas', 'Er th', 'Er meas');
snap = {};
for ic = 1:size(cases, 1)
  c1 = cases{ic,2}; r1 = cases{ic,3}; c2 = cases{ic,4}; r2 = cases{ic,5};
  dL = cases{ic,6}; cL = cases{ic,7}; rL = cases{ic,8};
  cmin = min([c1 c2]); cmax = max([c1 c2]); if dL > 0, cmin = min(cmin, cL); cmax = max(cmax, cL); end
  dx = cmin/f0/PPW; if dL > 0, dx = dL/ceil(dL/dx); end   % layer = whole number of cells
  dt = 0.95*dx/cmax;
  xI = (round(15e-3/dx) - 0.5)*dx;     % interface position (m), midway between p nodes
  x0 = xI - 9e-3;                      % initial pulse centre
  xA = xI - 4e-3; xB = xI + dL + 2e-3; % receivers in medium 1 and medium 2
  tEnd = (xA - x0 + 8e-3)/c1 + 3e-6;   % reflected pulse (and layer ring-down) past receiver A
  xEnd = xB + c2*tEnd + 2e-3;
  x  = (0:dx:xEnd)';  nx = numel(x);  xh = x + dx/2;      % p at x, u at x + dx/2
  c = c2*ones(nx,1); rho = r2*ones(nx,1);
  c(x < xI) = c1; rho(x < xI) = r1;
  if dL > 0, in = (x >= xI) & (x < xI + dL); c(in) = cL; rho(in) = rL; end
  K = rho.*c.^2;                                  % bulk modulus at p nodes
  rh = rho; rh(1:end-1) = 0.5*(rho(1:end-1) + rho(2:end));   % density at u nodes
  % initial right-going pulse: p(x,0) = s(-(x-x0)/c1); u = p/Z1 half a time step later
  p = pulse(-(x - x0)/c1);                 p(x >= xI) = 0;
  u = pulse(-(xh - x0)/c1 + dt/2)/(r1*c1); u(xh >= xI) = 0;   % u lives at t = +dt/2 (leapfrog)
  nt = ceil(tEnd/dt); iA = round(xA/dx) + 1; iB = round(xB/dx) + 1;
  pA = zeros(nt,1); pB = zeros(nt,1);
  for n = 1:nt
    p(2:end) = p(2:end) - dt*K(2:end).*(u(2:end) - u(1:end-1))/dx;          % dp/dt = -K du/dx
    u(1:end-1) = u(1:end-1) - dt./rh(1:end-1).*(p(2:end) - p(1:end-1))/dx;  % du/dt = -(1/rho) dp/dx
    pA(n) = p(iA); pB(n) = p(iB);
    if n == round((xI - x0 + 3e-3)/c1/dt), snap{ic} = [x p]; end %#ok<SAGROW>
  end
  t = (1:nt)'*dt; tm = (xA - x0 + 4e-3)/c1;          % split incident / reflected at receiver A
  inc = pA.*(t < tm); ref = pA.*(t >= tm);
  % pressure coefficients by projection onto the incident waveform (shifted to each arrival)
  shiftTo = @(sig, d) interp1(t, sig, t - d, 'linear', 0);
  incR = shiftTo(inc, 8e-3/c1);                              % incident shape at the reflected arrival
  incT = shiftTo(inc, (xI - xA)/c1 + dL*(dL > 0)/max(cL,1) + (xB - xI - dL)/c2);
  % allow a small residual shift (the material jump sits between grid nodes): best-aligned projection
  proj = @(sig, ref0) max(arrayfun(@(d) abs(sum(sig.*shiftTo(ref0, d))), (-40:40)*dt/4)) / sum(ref0.^2);
  Rm = sign(sum(ref.*incR))*proj(ref, incR); Tm = proj(pB, incT);
  Erm = sum(ref.^2)/sum(inc.^2);                             % reflected energy flux ratio (same medium)
  Rth = (r2*c2 - r1*c1)/(r2*c2 + r1*c1); Tth = 1 + Rth;
  if dL == 0
    Erth = Rth^2;
  else  % band-averaged energy reflection of the layer, weighted by the pulse power spectrum
    f = linspace(f0 - 4*sf, f0 + 4*sf, 801); f = f(f > 0); w = exp(-(f - f0).^2/sf^2);
    Za = r1*c1; Zb = r2*c2; Zl = rL*cL; kd = 2*pi*f*dL/cL;
    Zin = Zl*(Zb + 1i*Zl*tan(kd))./(Zl + 1i*Zb*tan(kd));    % input impedance seen from medium 1
    Rf = (Zin - Za)./(Zin + Za); Erth = sum(w.*abs(Rf).^2)/sum(w);
    Rth = NaN; Tth = NaN; Rm = NaN; Tm = NaN;                % waveform changes: compare energy only
  end
  fprintf('%-26s %8.4f %8.4f %8.4f %8.4f %9.4f %9.4f\n', cases{ic,1}, Rth, Rm, Tth, Tm, Erth, Erm);
end
fprintf('(T = 1 + R is a pressure ratio and can exceed 1; the energy fractions are R^2 and 1 - R^2.)\n');

%% Figure: snapshots
figure('Position', [50 50 900 640]);
for ic = 1:size(cases, 1)
  subplot(size(cases,1), 1, ic); s = snap{ic};
  plot(s(:,1)*1e3, s(:,2), 'b'); hold on; yl = [-1.6 1.6]; ylim(yl);
  plot([15 15], yl, 'k--'); if cases{ic,6} > 0, plot([15 15]+cases{ic,6}*1e3, yl, 'k:'); end
  xlim([0 35]); ylabel('p / p_{inc}'); title(cases{ic,1});
end
xlabel('x (mm)');

%% Oblique incidence: fluid (tissue) -> isotropic solid (skull), energy partition
% Potentials: fluid phi1 = e^{i(kx x + kz1 z)} + R e^{i(kx x - kz1 z)};
% solid phi = A e^{i(kx x + kzL z)}, psi = B e^{i(kx x + kzS z)}. Boundary conditions at z = 0:
% u_z continuous, sigma_zz = -p, sigma_xz = 0.
r1 = 1050; c1 = 1540; r2 = 1900; cLb = 2820; cSb = 1500;   % White, Clement & Hynynen 2006 speeds (1 MHz)
lam = r2*(cLb^2 - 2*cSb^2); mu = r2*cSb^2; w = 1;
th = (0:1:89)*pi/180; ER = zeros(size(th)); EL = ER; ES = ER;
for i = 1:numel(th)
  k1 = w/c1; kL = w/cLb; kS = w/cSb; kx = k1*sin(th(i));
  kz1 = sqrt(k1^2 - kx^2); kzL = sqrt(kL^2 - kx^2 + 0i); kzS = sqrt(kS^2 - kx^2 + 0i);
  if imag(kzL) < 0, kzL = -kzL; end
  if imag(kzS) < 0, kzS = -kzS; end
  M = [1i*kz1, 1i*kzL, 1i*kx;
       r1*w^2, -lam*kL^2 - 2*mu*kzL^2, -2*mu*kx*kzS;
       0, -2*kx*kzL, kzS^2 - kx^2];
  v = M \ [1i*kz1; -r1*w^2; 0];
  ER(i) = abs(v(1))^2;
  EL(i) = r2*cLb*kL^2*real(kzL/kL)/(r1*c1*k1^2*cos(th(i)))*abs(v(2))^2;
  ES(i) = r2*cSb*kS^2*real(kzS/kS)/(r1*c1*k1^2*cos(th(i)))*abs(v(3))^2;
end
fprintf('\nTissue -> skull (cL %d, cS %d m/s): longitudinal critical angle %.1f deg\n', cLb, cSb, asind(c1/cLb));
fprintf(' angle  E_refl  E_long  E_shear   sum\n');
for a = [0 20 30 40 50 60 70 80]
  i = a + 1; fprintf('%5d  %6.3f  %6.3f  %6.3f  %6.4f\n', a, ER(i), EL(i), ES(i), ER(i) + EL(i) + ES(i));
end
figure('Position', [80 80 640 360]);
plot(th*180/pi, ER, 'k', th*180/pi, EL, 'b', th*180/pi, ES, 'r', 'LineWidth', 1.4);
xlabel('incidence angle (deg)'); ylabel('energy fraction'); grid on;
legend('reflected', 'longitudinal in bone', 'shear in bone', 'Location', 'east');
title('Soft tissue -> skull: mode conversion');

%% m15_attenuation_dispersion.m
% Purpose : power-law attenuation alpha(f) = alpha0 f^y, the dispersion that causality
%           (Kramers-Kronig, nearly local form) attaches to it, the downshift of a
%           Gaussian pulse, and the pulse after propagation (FFT model).
% Change  : alpha0 (dB/cm/MHz), y, depth z, f0, FBW.
% Expect  : y = 1, alpha0 = 0.5: c(10 MHz) - c(1 MHz) = 3.19 m/s (1.38 m/s per neper of ln f)
%           y = 1.1: 3.56 m/s ; y = 1.5: 4.71 m/s
%           5 MHz, 60 % pulse, 0.5 dB/cm/MHz: centroid shift 0.467 MHz per 5 cm one way
%           (sigma_f^2 * alpha0' * z), 0.93 MHz for a 5 cm round trip
% Runtime : under 1 s.
clear; close all; rng(1);

%% Dispersion from Kramers-Kronig (Szabo 1994 / Waters et al. 2000 nearly local form)
c0 = 1540;                 % phase speed at the reference frequency f_ref (m/s)
fref = 1e6;                % reference frequency (Hz)
f = linspace(0.5e6, 15e6, 600);
dBperNp = 20*log10(exp(1));          % 8.686 dB per neper
cKK = @(a0dB, y, fq) 1./( 1/c0 + ...
     (abs(y-1) < 1e-6) .* ( -2/pi*(a0dB/dBperNp*100/(2*pi*1e6)^y) .* log(fq/fref) ) + ...
     (abs(y-1) >= 1e-6) .* ( (a0dB/dBperNp*100/(2*pi*1e6)^y) .* tan(pi*y/2) .* ((2*pi*fq).^(y-1) - (2*pi*fref).^(y-1)) ) );
cases = [0.5 1.0; 0.5 1.1; 0.5 1.5; 1.0 1.0; 0.5 1.0001];
fprintf('alpha0 (dB/cm/MHz)   y     c(5)-c(1)   c(10)-c(1)  (m/s)\n');
for i = 1:size(cases,1)
  a0 = cases(i,1); y = cases(i,2);
  fprintf('   %5.2f            %6.4f   %7.2f     %7.2f\n', a0, y, cKK(a0,y,5e6) - c0, cKK(a0,y,10e6) - c0);
end
fprintf('y -> 1 limit check: the y = 1.0001 row must equal the y = 1 row.\n');

%% Gaussian pulse downshift: amplitude spectrum exp(-(f-f0)^2/(2 sf^2)) * exp(-alpha0' f z)
f0 = 5e6; FBW = 0.6; sf = FBW*f0/(2*sqrt(2*log(2)));       % sf = 1.274 MHz
a0 = 0.5; a0p = a0/dBperNp*100/1e6;                          % Np/(m*Hz)
fg = linspace(0, 15e6, 3001);
fprintf('\n5 MHz, 60 %% pulse, alpha0 = 0.5 dB/cm/MHz, y = 1 (one-way):\n');
fprintf(' z (cm)  centroid (MHz)  theory f0 - sf^2 a0 z (MHz)   loss at f0 (dB)\n');
for zc = [0 2.5 5 10]
  z = zc/100; S = exp(-(fg - f0).^2/(2*sf^2)).*exp(-a0p*fg*z);
  fcen = sum(fg.*S)/sum(S);
  fprintf(' %5.1f      %6.3f          %6.3f                        %5.1f\n', zc, fcen/1e6, (f0 - sf^2*a0p*z)/1e6, a0*f0/1e6*zc);
end
fprintf('Round trip to 5 cm (10 cm path): shift %.2f MHz.\n', sf^2*a0p*0.10/1e6);

%% Pulse after propagation (time domain via FFT), y = 1.1 for a visible dispersion
y = 1.1; a0 = 0.5; zc = 8; z = zc/100;
fs = 100e6; N = 4096; t = ((0:N-1) - N/2)/fs;
st = 1/(2*pi*sf); s = exp(-t.^2/(2*st^2)).*cos(2*pi*f0*t);
F = fft(ifftshift(s)); fk = (0:N-1)*fs/N; fk(fk > fs/2) = fk(fk > fs/2) - fs;
af = abs(fk)/1e6; alpha = a0/dBperNp*100*af.^y;              % Np/m
cf = cKK(a0, y, max(abs(fk), 1e3)); cref = cKK(a0, y, f0);
ph = -2*pi*fk*z.*(1./cf - 1/cref);                             % phase relative to a frame moving at c(f0)
G = F.*exp(-alpha*z).*exp(1i*ph);
r = fftshift(real(ifft(G)));
% centroid and envelope peak
R = abs(fft(ifftshift(r))); half = fk > 0;
fprintf('\nAfter %d cm with y = %.1f: spectral centroid %.2f MHz, peak amplitude %.3f of input\n', zc, y, ...
  sum(fk(half).*R(half).^1)/sum(R(half))/1e6, max(abs(r)));
figure('Position', [60 60 900 600]);
subplot(2,2,1); plot(f/1e6, cKK(0.5,1,f) - c0, 'b', f/1e6, cKK(0.5,1.1,f) - c0, 'r', f/1e6, cKK(1,1,f) - c0, 'k', 'LineWidth', 1.3);
xlabel('f (MHz)'); ylabel('c(f) - c(1 MHz) (m/s)'); grid on; legend('0.5, y=1', '0.5, y=1.1', '1.0, y=1', 'Location', 'southeast');
title('Kramers-Kronig dispersion');
subplot(2,2,2); hold on;
for zc2 = [0 5 10]
  S = exp(-(fg - f0).^2/(2*sf^2)).*exp(-a0p*fg*zc2/100); plot(fg/1e6, S, 'LineWidth', 1.3);
end
xlabel('f (MHz)'); ylabel('amplitude'); legend('0 cm', '5 cm', '10 cm'); grid on; title('Downshift, 0.5 dB/cm/MHz');
subplot(2,1,2); plot(t*1e6, s, 'Color', [0.6 0.6 0.6]); hold on; plot(t*1e6, r/max(abs(r)), 'b', 'LineWidth', 1.3);
xlim([-1.2 1.2]); xlabel('t - z/c(f_0) (\mus)'); ylabel('normalized p');
legend('input', sprintf('after %d cm (normalized)', zc)); grid on; title('Attenuated, dispersed pulse');

%% m15_fubini_harmonics.m
% Purpose : plane-wave nonlinear propagation. Shock distance xbar = rho c^3/(beta omega p0),
%           Fubini harmonics B_n = 2 J_n(n sigma)/(n sigma), and a coupled-harmonic solution of
%           the Burgers equation with power-law absorption (Haran & Cook style):
%             dB_n/dsigma = (n/4) sum_{k=1}^{n-1} B_k B_{n-k} - (n/2) sum_{k>=1} B_k B_{k+n} - alpha_n xbar B_n
%           where p/p0 = sum B_n sin(n omega tau), sigma = x/xbar.
% Change  : medium table, p0, f0, number of harmonics N.
% Expect  : water (B/A 5.0) 3 MHz 1 MPa: xbar = 49.2 mm; liver (B/A 7.0) 49.1 mm; fat (B/A 10) 25.6 mm
%           Fubini at sigma = 0.5: B2 = 0.230 (B2/B1 = -12.5 dB); ODE lossless matches to 1e-4
%           liver with absorption (Gamma = 1.06): B2 peaks near 24 mm (23.6) at -20.9 dB re p0
% Runtime : about 3 s.
clear; close all; rng(1);

%% Media: rho (kg/m^3), c (m/s), alpha0 (dB/cm/MHz^y), y, B/A
names = {'water 20 C', 'blood', 'liver', 'muscle', 'fat'};
M = [ 998 1482 0.0022 2.0  5.0;
     1060 1584 0.15   1.2  6.1;
     1060 1578 0.5    1.1  7.0;
     1050 1580 1.0    1.1  7.4;
      950 1450 0.6    1.0 10.0];
f0 = 3e6; p0 = 1e6; dBperNp = 20*log10(exp(1));
fprintf('%-11s  beta   xbar(mm)  Mach eps   alpha(f0) Np/m   Gamma = 1/(alpha xbar)\n', 'medium');
for i = 1:numel(names)
  rho = M(i,1); c = M(i,2); beta = 1 + M(i,5)/2;
  xbar = rho*c^3/(beta*2*pi*f0*p0);
  alpha = M(i,3)*(f0/1e6)^M(i,4)/dBperNp*100;
  fprintf('%-11s  %4.2f   %6.1f    %.2e      %7.2f          %6.2f\n', names{i}, beta, xbar*1e3, p0/(rho*c)/c, alpha, 1/(alpha*xbar));
end

%% Fubini versus the coupled-harmonic ODE (lossless), RK4
N = 40;                                  % harmonics kept
nn = (1:N)';
fub = @(n, s) 2*besselj(n, n*s)./(n*s);
s1 = @(cc) [0; cc(1:N-1)];               % sum_{k=1}^{n-1} B_k B_{n-k}  from conv(B,B)
s2 = @(r) [r(N-1:-1:1); 0];              % sum_k B_k B_{k+n}           from conv(B,flipud(B))
dBds = @(B, att) (nn/4).*s1(conv(B, B)) - (nn/2).*s2(conv(B, flipud(B))) - att.*B;
B = [1; zeros(N-1,1)]; h = 0.0025; att0 = zeros(N,1); sig = 0;
for k = 1:round(0.9/h)
  k1 = dBds(B, att0); k2 = dBds(B + h/2*k1, att0); k3 = dBds(B + h/2*k2, att0); k4 = dBds(B + h*k3, att0);
  B = B + h/6*(k1 + 2*k2 + 2*k3 + k4); sig = sig + h;
  if abs(sig - 0.5) < h/2, B05 = B; end
end
fprintf('\nLossless check (N = %d):   B1       B2       B3\n', N);
fprintf(' Fubini  sigma = 0.5   %.4f   %.4f   %.4f\n', fub(1,0.5), fub(2,0.5), fub(3,0.5));
fprintf(' ODE     sigma = 0.5   %.4f   %.4f   %.4f\n', B05(1), B05(2), B05(3));
fprintf(' Fubini  sigma = 0.9   %.4f   %.4f   %.4f\n', fub(1,0.9), fub(2,0.9), fub(3,0.9));
fprintf(' ODE     sigma = 0.9   %.4f   %.4f   %.4f\n', B(1), B(2), B(3));
fprintf(' small sigma: B2 ~ sigma/2 ; at sigma = 0.1: Fubini %.4f\n', fub(2, 0.1));

%% Liver with power-law absorption: split step (exact absorption, RK4 nonlinearity)
N = 60; nn = (1:N)';
s1 = @(cc) [0; cc(1:N-1)]; s2 = @(r) [r(N-1:-1:1); 0];
dBnl = @(B) (nn/4).*s1(conv(B, B)) - (nn/2).*s2(conv(B, flipud(B)));
rho = 1060; c = 1578; a0 = 0.5; y = 1.1; beta = 1 + 7.0/2;
xbar = rho*c^3/(beta*2*pi*f0*p0);
alphan = a0*(nn*f0/1e6).^y/dBperNp*100;          % Np/m for each harmonic
xmax = 0.12; h = 0.005; ns = round(xmax/xbar/h); hx = xmax/ns; hs = hx/xbar;
E = exp(-alphan*hx/2);                            % half-step absorption factor
B = [1; zeros(N-1,1)]; X = (0:ns)*hx; H = zeros(4, ns+1); H(:,1) = B(1:4);
for k = 1:ns
  B = E.*B;
  k1 = dBnl(B); k2 = dBnl(B + hs/2*k1); k3 = dBnl(B + hs/2*k2); k4 = dBnl(B + hs*k3);
  B = B + hs/6*(k1 + 2*k2 + 2*k3 + k4);
  B = E.*B; H(:,k+1) = B(1:4);
end
[b2max, im] = max(abs(H(2,:)));
fprintf('\nLiver, 3 MHz, 1 MPa: xbar = %.1f mm, Gamma = %.2f\n', xbar*1e3, 1/(a0*3^y/dBperNp*100*xbar));
fprintf(' x (mm)   B1 (dB re p0)   B2 (dB re p0)   B2/B1 (dB)\n');
for xm = [10 25 50 75 100]
  [~, j] = min(abs(X - xm/1e3));
  fprintf(' %5.0f      %6.1f          %6.1f         %6.1f\n', xm, 20*log10(abs(H(1,j))), 20*log10(abs(H(2,j))), 20*log10(abs(H(2,j)/H(1,j))));
end
fprintf(' B2 maximum %.1f dB re p0 at x = %.1f mm\n', 20*log10(b2max), X(im)*1e3);

%% Figures
figure('Position', [60 60 900 600]);
subplot(2,2,1); ss = linspace(1e-3, 1, 200); hold on;
for n = 1:4, plot(ss, fub(n, ss), 'LineWidth', 1.3); end
xlabel('\sigma = x / xbar'); ylabel('B_n'); legend('n=1','n=2','n=3','n=4'); grid on; title('Fubini (lossless)');
subplot(2,2,2); plot(X*1e3, 20*log10(abs(H')), 'LineWidth', 1.3); ylim([-60 0]);
xlabel('x (mm)'); ylabel('dB re p_0'); legend('f_0','2f_0','3f_0','4f_0'); grid on; title('Liver, 3 MHz, 1 MPa, with absorption');
subplot(2,1,2); th = linspace(0, 4*pi, 800); hold on;
for sg = [0 0.5 0.95]
  % implicit Poisson solution P = sin(theta + sigma P), by fixed-point iteration (converges for sigma < 1)
  P = sin(th); for it = 1:200, P = sin(th + sg*P); end
  plot(th/(2*pi), P, 'LineWidth', 1.3);
end
xlabel('\omega\tau / 2\pi (periods)'); ylabel('p / p_0'); legend('\sigma = 0', '0.5', '0.95 (near shock)'); grid on;
title('Waveform steepening (lossless, implicit solution)');

%% m15_harmonic_beam.m
% Purpose : why tissue-harmonic beams are narrower and cleaner. Quasi-linear focal-plane model:
%           the second-harmonic source is proportional to the square of the fundamental
%           transmit pattern, and it is received at 2 f0 by the same array.
%             fundamental two-way : [E_f0 AF_f0]_tx * [E_f0 AF_f0]_rx
%             harmonic two-way    : [E_f0 AF_f0]_tx^2 * [E_2f0 AF_2f0]_rx
%             linear at 2 f0      : [E_2f0 AF_2f0]^2   (same probe driven at 2 f0)
% Change  : f0, number of elements, pitch, steering angle, apodization.
% Expect  : continuous apertures (u = x D / (lambda0 F)): -6 dB full widths
%           sinc 1.207, sinc^2 0.886 (x 0.734), sinc^2(u) sinc(2u) 0.511, sinc^2(2u) 0.443;
%           first sidelobes -13.3, -26.5, -27.9, -26.5 dB.
%           Array, 3 MHz, 64 elements, pitch 0.3 mm, steered 30 deg: the 2 f0 receive pattern has a
%           grating lobe near -20.8 deg; in the harmonic two-way pattern it is suppressed by the
%           low transmit level there (see printed levels).
% Runtime : under 1 s.
clear; close all; rng(1);

%% Continuous aperture: widths and sidelobes
u = linspace(0, 6, 600001);
sincf = @(x) (abs(x) < 1e-12) + (abs(x) >= 1e-12).*sin(pi*x)./(pi*x + (abs(x) < 1e-12));
pats = {sincf(u), sincf(u).^2, sincf(u).^2.*sincf(2*u), sincf(2*u).^2};
lab  = {'fundamental one-way sinc', 'fundamental two-way / harmonic tx', 'harmonic two-way', 'linear at 2f0 two-way'};
fprintf('%-36s  -6 dB FW (lambda0 F/D)   peak sidelobe (dB)\n', 'pattern');
for i = 1:4
  a = abs(pats{i}); i6 = find(a < 0.5, 1); d = diff(a); i0 = find(d(1:end-1) < 0 & d(2:end) >= 0, 1) + 1;
  fprintf('%-36s       %6.3f                %7.2f\n', lab{i}, 2*u(i6), 20*log10(max(a(i0:end))));
end

%% Linear array in the focal plane (far-field array factor in sin(theta))
c = 1540; f0 = 3e6; Nel = 64; pitch = 0.3e-3; kerf = 0.03e-3; th0 = 30; apod = 'rect';
w = pitch - kerf; xn = ((0:Nel-1) - (Nel-1)/2)*pitch;
if strcmp(apod, 'hann'), wn = 0.5 - 0.5*cos(2*pi*((0:Nel-1) + 1)/(Nel + 1)); else, wn = ones(1, Nel); end
th = linspace(-90, 90, 3601); s = sind(th); s0 = sind(th0);
AF = @(f) abs(exp(1i*2*pi*f/c*(s(:) - s0)*xn)*wn(:)).'/sum(wn);
EF = @(f) abs(sincf(w*f/c*s));                          % rigid-baffle element factor
P1 = EF(f0).*AF(f0); P2 = EF(2*f0).*AF(2*f0);
fund = P1.^2; harm = P1.^2.*P2; lin2 = P2.^2;
ipk = find(abs(th - th0) < 2);                          % normalize to the steered main lobe
dBn = @(p) 20*log10(p/max(p(ipk)) + 1e-12);
thg1 = asind(max(-1, min(1, s0 - c/f0/pitch))); thg2 = asind(s0 - c/(2*f0)/pitch);
fprintf('\nArray: %d el, pitch %.2f mm, steer %d deg, f0 %.1f MHz\n', Nel, pitch*1e3, th0, f0/1e6);
if abs(s0 - c/f0/pitch) <= 1, g1 = sprintf('%.1f deg', thg1); else, g1 = 'none (|sin| > 1)'; end
fprintf(' predicted grating lobe: tx/rx at f0 %s ; rx at 2f0 %.1f deg\n', g1, thg2);
win = abs(th - thg2) < 3;
names = {'fundamental two-way', 'harmonic two-way', 'linear 2f0 two-way', 'rx only at 2f0'}; P = {fund, harm, lin2, P2};
fprintf(' %-22s  -6 dB width (deg)   level at the 2f0 grating lobe (dB)\n', 'pattern');
for i = 1:4
  d = dBn(P{i}); [~, k0] = min(abs(th - th0)); a = k0; b = k0;   % contiguous -6 dB region around th0
  while a > 1 && d(a-1) >= -6, a = a - 1; end
  while b < numel(th) && d(b+1) >= -6, b = b + 1; end
  fprintf(' %-22s      %6.2f             %7.1f\n', names{i}, th(b) - th(a), max(d(win)));
end

figure('Position', [60 60 900 420]);
plot(th, dBn(fund), 'b', th, dBn(harm), 'r', th, dBn(lin2), 'k--', 'LineWidth', 1.3);
ylim([-80 0]); xlim([-90 90]); grid on; xlabel('angle (deg)'); ylabel('two-way level (dB)');
legend('fundamental', 'harmonic (tx^2 x rx at 2f_0)', 'linear at 2f_0', 'Location', 'southwest');
title(sprintf('Two-way beam, %d elements, pitch %.2f mm, steered %d deg', Nel, pitch*1e3, th0));

%% m15_pulse_inversion.m
% Purpose : pulse inversion (PI) and amplitude modulation (AM) as Volterra-order filters.
%           Tissue: the transmitted pulse distorts by lossless plane-wave propagation (simple-wave
%           solution, each sample of the waveform travels at c0 + beta u), then scatters linearly.
%           "Bubble-like": a memoryless polynomial echo y = x + a2 x^2 + a3 x^3 (illustrative only).
%           PI  : e(+s) + e(-s)          -> odd orders cancel, even orders double
%           AM  : e(s) - 2 e(s/2)        -> linear cancels, quadratic x 1/2, cubic x 3/4
%           Motion between pulses leaves 2|sin(pi f0 delta)| of the fundamental, delta = 2 v / (c PRF).
% Change  : f0, FBW, sigma (distortion = x / xbar), a2, a3, velocity and PRF.
% Expect  : tissue, sigma = 0.3: single echo 2f0/f0 band ratio about -18 dB; PI removes the fundamental
%           band to the level of the even-order spill (< -30 dB re single), 2f0 band +5.8 dB (ideal +6.0);
%           motion residual -26.2 dB at 10 mm/s, PRF 5 kHz, 3 MHz (formula and simulation agree).
%           Band separation: Gaussian f0 and 2f0 spectra meet at -20 dB when FBW = 45 %.
% Runtime : under 1 s.
clear; close all; rng(1);

%% Pulse and helpers
c = 1540; f0 = 3e6; FBW = 0.5; fs = 200e6; N = 4096;
t = ((0:N-1) - N/2)/fs; sf = FBW*f0/(2*sqrt(2*log(2))); st = 1/(2*pi*sf);
s = exp(-t.^2/(2*st^2)).*sin(2*pi*f0*t);                % unit-amplitude transmit pulse
fk = (0:N-1)*fs/N;
bandE = @(x, f1, f2) sum(abs(fft(x)).^2 .* (fk >= f1 & fk < f2));
Ef = @(x) bandE(x, 0.5*f0, 1.5*f0); Eh = @(x) bandE(x, 1.5*f0, 2.5*f0);
dB10 = @(r) 10*log10(r);
% lossless simple-wave distortion: a sample of value A*s(t) arrives at t - sigma*A*s(t)/(2 pi f0)
% (sigma is defined for amplitude 1; amplitude A scales the distortion by A). Valid before the shock.
prop = @(A, sig) interp1(t - sig*A*s/(2*pi*f0), A*s, t, 'linear', 0);

%% Tissue: PI and AM
sig = 0.3;
eP = prop(1, sig); eN = prop(-1, sig); eH = prop(0.5, sig);
pi_out = eP + eN; am_out = eP - 2*eH;
fprintf('Tissue, f0 = %.1f MHz, FBW = %.0f %%, sigma = %.2f\n', f0/1e6, FBW*100, sig);
fprintf(' single echo: 2f0 band re f0 band          %6.1f dB\n', dB10(Eh(eP)/Ef(eP)));
fprintf(' PI : f0 band re single f0 band             %6.1f dB\n', dB10(Ef(pi_out)/Ef(eP)));
fprintf(' PI : 2f0 band re single 2f0 band           %6.1f dB  (ideal +6.0)\n', dB10(Eh(pi_out)/Eh(eP)));
fprintf(' AM : f0 band re single f0 band             %6.1f dB\n', dB10(Ef(am_out)/Ef(eP)));
fprintf(' AM : 2f0 band re single 2f0 band           %6.1f dB  (quadratic retained x 1/2: -6.0)\n', dB10(Eh(am_out)/Eh(eP)));

%% Bubble-like polynomial echo
a2 = 0.3; a3 = 0.2; ech = @(x) x + a2*x.^2 + a3*x.^3;
bP = ech(s); bN = ech(-s); bH = ech(0.5*s);
fprintf('\nPolynomial echo (a2 = %.1f, a3 = %.1f): PI keeps even orders, AM also keeps cubic terms\n', a2, a3);
fprintf(' PI : f0 band re single f0 band             %6.1f dB\n', dB10(Ef(bP + bN)/Ef(bP)));
fprintf(' AM : f0 band re single f0 band             %6.1f dB  (nonlinear fundamental from x^3)\n', dB10(Ef(bP - 2*bH)/Ef(bP)));
fprintf(' AM : expected cubic-fundamental level      %6.1f dB\n', dB10(Ef(0.75*a3*s.^3)/Ef(bP)));

%% Motion between pulses
fprintf('\nMotion residual of the fundamental after PI (linear medium), f0 = %.1f MHz\n', f0/1e6);
fprintf('  v (mm/s)  PRF (kHz)   delta (ns)   formula (dB)   simulated (dB)\n');
for vv = [1 10 50]
  for prf = [5e3]
    d = 2*vv*1e-3/c/prf; s2 = -interp1(t, s, t - d, 'linear', 0);
    fprintf('  %6.0f     %5.1f      %7.2f       %7.1f        %7.1f\n', vv, prf/1e3, d*1e9, ...
      20*log10(2*abs(sin(pi*f0*d))), dB10(Ef(s + s2)/Ef(s)));
  end
end

%% Band separation for filter-based harmonic imaging (Gaussian spectra)
% fundamental half-width at L dB: (B/2) f0 sqrt(L/6); the 2f0 spectrum of s^2 is sqrt(2) wider in Hz
for L = [6 20 30]
  fprintf('f0 and 2f0 spectra meet at -%2d dB when FBW = %.0f %%\n', L, 100*2/((1 + sqrt(2))*sqrt(L/6)));
end

%% Figure
figure('Position', [60 60 900 600]);
win = abs(t) < 1.2e-6;
subplot(2,2,1); plot(t(win)*1e6, eP(win), 'b', t(win)*1e6, eN(win), 'r'); grid on;
xlabel('t (\mus)'); ylabel('echo'); legend('+s', '-s'); title(sprintf('Tissue echoes, \\sigma = %.1f', sig));
subplot(2,2,2); plot(t(win)*1e6, pi_out(win), 'k', t(win)*1e6, am_out(win), 'm'); grid on;
xlabel('t (\mus)'); legend('PI sum', 'AM: e(s) - 2e(s/2)'); title('Combined outputs');
subplot(2,1,2); half = fk < 12e6; ref = max(abs(fft(eP)));
S1 = 20*log10(abs(fft(eP))/ref + 1e-9); S2 = 20*log10(abs(fft(pi_out))/ref + 1e-9); S3 = 20*log10(abs(fft(am_out))/ref + 1e-9);
plot(fk(half)/1e6, S1(half), 'b', fk(half)/1e6, S2(half), 'k', fk(half)/1e6, S3(half), 'm', 'LineWidth', 1.2);
ylim([-80 5]); grid on; xlabel('f (MHz)'); ylabel('dB re single-echo peak');
legend('single echo', 'PI sum', 'AM'); title('Spectra: PI removes f_0 and doubles 2f_0');