Sampling, Fourier & LTI systems

%% m04_fft_1d_basics.m  -- Module 04: frequency axes, aliasing, zero padding, Parseval, windows
% Purpose : build the FFT frequency axis correctly, see aliasing, show that zero padding
%           interpolates but does not resolve (a longer record does), check Parseval, and
%           measure the figures of merit of the classic windows.
% Change  : fs, N, the tone frequencies, the tone separation sepT (in units of 1/T),
%           the window list.
% Expect  : alias of 7 MHz at fs = 10 MHz -> 3.000 MHz; legal bandpass rates: [3.5, 6.5] MHz ->
%           fs >= 13 or 6.5..7.0 MHz, [3, 7] MHz -> only fs >= 14 MHz (band edges at 6.75 MHz
%           land at 3.25 and 0.25 MHz); padded short record shows 1 peak,
%           4x longer record shows 2; Parseval gap ~1e-12; window table
%           rect -13.3 dB / Hann -31.5 dB / Hamming -42.7 dB / Blackman -58.1 dB,
%           roll-off about -6, -18, -6, -18 dB/octave.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% 1. Frequency axis and fftshift
fs = 40;                 % sampling rate (MHz)
N  = 1024;               % record length (samples)
T  = N/fs;               % record duration (us)
df = fs/N;               % bin spacing (MHz)
f_unshift = (0:N-1)*df;                   % bins as returned by fft (0 ... fs)
f_shift   = (-N/2:N/2-1)*df;              % axis after fftshift (-fs/2 ... fs/2)
n  = 0:N-1;
x  = cos(2*pi*5*n/fs + 0.3);              % 5 MHz tone
X  = fftshift(fft(x));
[~, ipk] = max(abs(X(N/2+1:end)));        % positive half only
fprintf('--- frequency axis ---\n');
fprintf('fs = %g MHz, N = %d: bin spacing fs/N = %.2f kHz, T = %.1f us, 1/T = %.2f kHz\n', fs, N, df*1e3, T, 1e3/T);
fprintf('5 MHz tone: peak at %.3f MHz (positive half of the fftshifted axis)\n', f_shift(N/2+ipk));
fprintf('unshifted bin k = N-1 is %.4f MHz = -df (negative frequency)\n', f_unshift(end) - fs);

%% 2. Aliasing: a 7 MHz tone sampled at 10 MHz
fsa = 10; Na = 4096; na = 0:Na-1;
w_hann = 0.5 - 0.5*cos(2*pi*na/Na);               % periodic Hann, written out (no toolbox)
xa = cos(2*pi*7*na/fsa + 0.6) .* w_hann;
Xa = abs(fft(xa)); [~, ka] = max(Xa(1:Na/2+1));
f_alias_formula = abs(7 - fsa*round(7/fsa));
fprintf('\n--- aliasing ---\n');
fprintf('7 MHz at fs = 10 MHz: formula |f0 - fs*round(f0/fs)| = %.3f MHz, FFT peak = %.3f MHz\n', ...
        f_alias_formula, (ka-1)*fsa/Na);

%% 2b. Bandpass sampling: legal rates 2 fH/n <= fs <= 2 fL/(n-1), n = 1..floor(fH/B)
fprintf('\n--- bandpass sampling ---\n');
bands = [3.5 6.5; 3 7];                            % [fL fH] in MHz
for ib = 1:size(bands, 1)
    fL = bands(ib, 1); fH = bands(ib, 2); Bw = fH - fL; nmax = floor(fH/Bw);
    fprintf('band [%.1f, %.1f] MHz, B = %.1f, n_max = %d:', fL, fH, Bw, nmax);
    for nz = 1:nmax
        if nz == 1, fprintf('  n=1: fs >= %.2f', 2*fH);
        else, fprintf('  n=%d: %.2f..%.2f', nz, 2*fH/nz, 2*fL/(nz-1)); end
    end
    fprintf(' MHz\n');
end
fsb = 6.75; fe = [3.5 6.5];                        % band edges sampled at 6.75 MHz (zone n = 2)
fa_e = abs(fe - fsb*round(fe/fsb));
fprintf('fs = %.2f MHz: band edges %.1f / %.1f MHz alias to %.2f / %.2f MHz (zone 2: inverted)\n', fsb, fe, fa_e);

%% 3. Zero padding vs a longer record (complex tones, as IQ data)
f1 = 5.0; sepT = 0.6;                     % tone separation in units of 1/T
Ns = 64; Ts = Ns/fs; f2 = f1 + sepT/Ts;   % 0.6/T apart for the SHORT record
mk = @(Nr) exp(1j*2*pi*f1*(0:Nr-1)/fs) + exp(1j*2*pi*f2*(0:Nr-1)/fs);
countPeaks = @(m) sum( m(2:end-1) > m(1:end-2) & m(2:end-1) >= m(3:end) & m(2:end-1) > 0.5*max(m) );
Npad = 4096;
Xs_pad  = abs(fft(mk(Ns), Npad));         % short record, 64x zero padded
Xl_pad  = abs(fft(mk(4*Ns), Npad));       % 4x longer record, same display grid
fax = (0:Npad-1)*fs/Npad;
sel = fax > f1 - 1.5 & fax < f2 + 1.5;
fprintf('\n--- zero padding vs longer record ---\n');
fprintf('tones %.3f and %.3f MHz (%.2f/T apart for N = %d)\n', f1, f2, sepT, Ns);
fprintf('N = %3d padded to %d: %d peak(s) (bins of %.4f MHz, resolution 1/T = %.3f MHz)\n', Ns, Npad, countPeaks(Xs_pad(sel)), fs/Npad, 1/Ts);
fprintf('N = %3d padded to %d: %d peak(s) (resolution 1/T = %.3f MHz)\n', 4*Ns, Npad, countPeaks(Xl_pad(sel)), 1/(4*Ts));

%% 4. Parseval
xr = randn(1, N);
gap = sum(abs(xr).^2) - sum(abs(fft(xr)).^2)/N;
fprintf('\n--- Parseval ---\nsum|x|^2 - sum|X|^2/N = %.2e\n', gap);

%% 5. Window figures of merit (DFT-even windows, N = 256, 64x zero padding)
Nw = 256; P = 64; L = Nw*P; m = 0:Nw-1; a = 2*pi*m/Nw;
I0 = @(z) sum( bsxfun(@rdivide, bsxfun(@power, z(:)/2, 2*(0:40)), factorial(0:40).^2), 2 ).';  % Bessel I0 series
names = {'Rectangular', 'Hann', 'Hamming', 'Blackman', 'Kaiser b=6'};
W = [ ones(1, Nw);
      0.5 - 0.5*cos(a);
      0.54 - 0.46*cos(a);
      0.42 - 0.5*cos(a) + 0.08*cos(2*a);
      I0(6*sqrt(max(0, 1 - (2*m/Nw - 1).^2))) / I0(6) ];
theoryPSL = [-13.26 -31.47 -42.66 -58.11 NaN];
fprintf('\n--- window table (computed) ---\n');
fprintf('%-12s %8s %8s %8s %8s %6s %6s %8s %10s\n', 'window', 'PSL dB', '-3dB bin', '-6dB bin', 'null bin', 'ENBW', 'CG', 'scal dB', 'roll dB/oct');
figure('Name', 'windows', 'Position', [100 100 760 520]);
cols = lines(5);
for iw = 1:size(W, 1)
    w = W(iw, :);
    S = abs(fft(w, L)); S = S/S(1);
    d = 20*log10(S(1:L/2) + 1e-300);
    k = 2; while ~(d(k) < d(k-1) && d(k) <= d(k+1)), k = k + 1; end     % first null
    psl = max(d(k:end));
    w3 = 2*interp1(d(1:k), (0:k-1)/P, -3.0103);
    w6 = 2*interp1(d(1:k), (0:k-1)/P, -6.0206);
    enbw = Nw*sum(w.^2)/sum(w)^2;
    cg = mean(w);
    scal = -d(P/2 + 1);                                                % tone half a bin off
    env = @(b) max(d(b*P+1:(b+1)*P));                                  % sidelobe peak near bin b
    roll = (env(40) - env(10))/2;                                      % 2 octaves, per octave
    fprintf('%-12s %8.2f %8.3f %8.3f %8.3f %6.3f %6.3f %8.2f %10.1f', names{iw}, psl, w3, w6, (k-1)/P, enbw, cg, scal, roll);
    if ~isnan(theoryPSL(iw)), fprintf('   (theory PSL %.2f)\n', theoryPSL(iw)); else, fprintf('\n'); end
    plot((0:L/2-1)/P, d, 'Color', cols(iw, :), 'LineWidth', 1.1); hold on;
end
hold off; xlim([0 12]); ylim([-110 3]); grid on;
xlabel('frequency offset (bins of N)'); ylabel('|W| (dB re peak)');
title('Window spectra: mainlobe width vs sidelobe level'); legend(names, 'Location', 'northeast');

%% 6. Figure: aliasing and zero padding
figure('Name', 'aliasing and zero padding', 'Position', [100 100 760 600]);
subplot(2, 1, 1);
tc = linspace(0, 1.6, 2000); ts = 0:1/fsa:1.6;
plot(tc, cos(2*pi*7*tc + 0.6), 'Color', [0.6 0.6 0.6]); hold on;
plot(tc, cos(2*pi*3*tc - 0.6), 'r--', 'LineWidth', 1.4);
plot(ts, cos(2*pi*7*ts + 0.6), 'ko', 'MarkerFaceColor', 'k'); hold off;
xlim([0 1.6]); xlabel('time (us)'); ylabel('amplitude');
title('7 MHz tone (gray) sampled at 10 MHz (dots) = 3 MHz alias, phase -\phi (red)');
subplot(2, 1, 2);
plot(fax(sel), Xs_pad(sel)/max(Xs_pad), 'b', 'LineWidth', 1.3); hold on;
plot(fax(sel), Xl_pad(sel)/max(Xl_pad), 'r', 'LineWidth', 1.3);
plot([f1 f1], [0 1.05], 'k:', [f2 f2], [0 1.05], 'k:'); hold off;
xlabel('frequency (MHz)'); ylabel('|X| (norm.)');
title('Zero padding (blue, N = 64) vs 4x longer record (red, N = 256)');
legend('N=64, padded 64x', 'N=256, padded 16x', 'Location', 'northeast');

%% m04_fft_2d_kspace.m  -- Module 04: 2D FFT and the k-space of a pulse-echo image
% Purpose : simulate a small RF image (speckle + point targets) with a broadband,
%           focal-plane (Fraunhofer) pulse-echo PSF, then look at log|FFT2|: the axial band
%           sits at kz = +-2 f0/c, the lateral extent is set by the transmit-receive
%           aperture autoconvolution, |kx| <= f/(c F#) at each frequency (a bow-tie).
%           Also: 2D FFT of a single PSF, of its envelope (B-mode PSF), and of separable windows.
% Change  : f0, FBW, Fnum, apod ('rect' or 'hann'), dx, dz, the scatterer density.
% Expect  : 2f0/c = 6.49 cycles/mm at 5 MHz; two-way lateral edge 1/(lambda F#) = 1.62
%           cycles/mm at F# 2; >99 % of the image k-space energy inside the bow-tie;
%           PSF -6 dB widths close to 0.886 lambda F# (rect) and c tau/2;
%           separable rect window: -13.3 dB sidelobes on the axes, -26.5 dB on the diagonal;
%           filtering: a k-space mask matched to the PSF passband keeps ~99 % of the signal energy
%           and raises the SNR of white noise by ~10 log10(1/mask area); a one-sided kz mask gives
%           the analytic image and its envelope (B-mode).
% Model   : paraxial, focal plane, Tx = Rx aperture, first-order Born, shift-invariant
%           over the small field (circular convolution). Base MATLAB only (R2019b+, Octave 8).
clear; close all; rng(1);

%% Parameters
c    = 1.540;      % speed of sound (mm/us)
f0   = 5;          % center frequency (MHz)
FBW  = 0.60;       % -6 dB fractional bandwidth of the two-way pulse
Fnum = 2;          % F-number z/D (same aperture on transmit and receive)
apod = 'rect';     % 'rect' or 'hann'
Nz = 128; Nx = 128;          % image size (axial x lateral)
dz = 0.02; dx = 0.06;        % pixel size (mm): axial Nyquist 25 c/mm, lateral 8.33 c/mm

B   = FBW*f0; sig_f = B/(2*sqrt(2*log(2)));          % amplitude-spectrum std (MHz)
lam0 = c/f0;
sincf = @(u) (sin(pi*u) + (u == 0)) ./ (pi*u + (u == 0));
if strcmp(apod, 'hann')
    A1 = @(u) sincf(u) + 0.5*sincf(u - 1) + 0.5*sincf(u + 1);   % FT of 0.5+0.5cos(2 pi xi/D), peak 1
else
    A1 = @(u) sincf(u);                                         % FT of a uniform aperture
end

%% Axes (fft order: origin at index 1) and shifted axes for display
kz_f = [0:Nz/2-1, -Nz/2:-1]/(Nz*dz);       % axial spatial frequency (cycles/mm)
x_f  = [0:Nx/2-1, -Nx/2:-1]*dx;            % lateral position (mm), origin first
kz_s = (-Nz/2:Nz/2-1)/(Nz*dz); kx_s = (-Nx/2:Nx/2-1)/(Nx*dx);
z_s  = (-Nz/2:Nz/2-1)*dz;      x_s  = (-Nx/2:Nx/2-1)*dx;

%% PSF: for every temporal frequency f = c kz/2, lateral two-way pattern A1(x/(lambda F#))^2
Hz = zeros(Nz, Nx);                         % axial-FFT x lateral-space representation
for k = 2:Nz/2                              % positive kz (skip DC)
    f  = kz_f(k)*c/2;                       % MHz
    Pf = exp(-(f - f0)^2/(2*sig_f^2));      % two-way pulse amplitude spectrum (on axis)
    u  = x_f/((c/f)*Fnum);                  % normalized lateral coordinate x/(lambda F#)
    Hz(k, :) = Pf * A1(u).^2;               % transmit pattern x receive pattern
end
Hz(Nz:-1:Nz/2+2, :) = Hz(2:Nz/2, :);        % Hermitian twin (real, zero-phase PSF)
h  = real(ifft(Hz, [], 1));                 % RF PSF h(z, x), centered at the origin
ha = ifft([Hz(1, :); 2*Hz(2:Nz/2, :); zeros(Nz/2, Nx)], [], 1);   % analytic PSF (axial)
H  = fft2(h);                               % 2D spectrum of the PSF

%% Scatterers: speckle + three point targets; image = PSF (*) scatterers (circular)
S = 0.5*randn(Nz, Nx) .* (rand(Nz, Nx) < 0.2);  % ~20 scatterers per resolution cell
pts = [32 30; 64 64; 100 100];              % (iz, ix) of the point targets
for ip = 1:size(pts, 1), S(pts(ip, 1), pts(ip, 2)) = 10; end
K   = H .* fft2(S);                         % k-space of the image (2D convolution theorem)
img = real(ifft2(K));
S1 = zeros(Nz, Nx); S1(Nz/2+1, Nx/2+1) = 1; % single point at the center
psf = real(ifft2(H .* fft2(S1)));

%% Measurements
envc = abs(ha); envc = envc/max(envc(:));
ex = fftshift(envc(1, :)); ez = fftshift(envc(:, 1)).';
cross = @(y, xa, lv, i) xa(i-1) + (lv - y(i-1)).*(xa(i) - xa(i-1))./(y(i) - y(i-1));   % first crossing below lv
xr = x_s(Nx/2+1:end); er = ex(Nx/2+1:end); zr = z_s(Nz/2+1:end); fr = ez(Nz/2+1:end);
lat6 = 2*cross(er, xr, 0.5, find(er < 0.5, 1));
ax6  = 2*cross(fr, zr, 0.5, find(fr < 0.5, 1));
tau = 4*log(2)/(pi*B);
kz0 = 2*f0/c; fL = f0 - B/2; fH = f0 + B/2;
kedge = 1/(lam0*Fnum);
[~, ir] = min(abs(kz_f - kz0));                                  % row nearest 2 f0/c
row = abs(H(ir, :)); row = fftshift(row/max(row)); rdB = 20*log10(row + 1e-300);
pos = Nx/2+1:Nx;
k30 = cross(rdB(pos), kx_s(pos), -30, find(rdB(pos) < -30, 1));  % -30 dB half-width
[KX, KZ] = meshgrid(kx_s, kz_s);
inside = abs(KX) <= abs(KZ)/(2*Fnum) + 1/(Nx*dx);                % bow-tie (+1 bin)
E = abs(fftshift(K)).^2;
fprintf('--- k-space of a pulse-echo image (f0 = %.1f MHz, FBW = %.0f %%, F# = %.1f, %s) ---\n', f0, 100*FBW, Fnum, apod);
fprintf('axial band center 2 f0/c = %.3f cycles/mm (RF period lambda/2 = %.4f mm)\n', kz0, lam0/2);
fprintf('axial -6 dB band 2fL/c .. 2fH/c = %.3f .. %.3f cycles/mm\n', 2*fL/c, 2*fH/c);
fprintf('lateral edge at f0: two-way 1/(lambda F#) = %.3f c/mm, one-way %.3f c/mm\n', kedge, kedge/2);
fprintf('measured |H| -30 dB lateral half-width at kz = %.2f c/mm: %.3f c/mm (rect triangle: %.3f)\n', kz_f(ir), k30, kedge*(kz_f(ir)/kz0)*(1 - 10^(-1.5)));
fprintf('fraction of image k-space energy inside the bow-tie |kx| <= |kz|/(2F#): %.5f\n', sum(E(inside))/sum(E(:)));
fprintf('PSF lateral -6 dB width: %.3f mm (narrowband theory %.3f mm = %.3f lambda F#)\n', lat6, ...
        lam0*Fnum*(0.886*strcmp(apod, 'rect') + 1.441*strcmp(apod, 'hann')), 0.886*strcmp(apod, 'rect') + 1.441*strcmp(apod, 'hann'));
fprintf('PSF axial   -6 dB width: %.3f mm (c tau/2 = %.3f mm)\n', ax6, c*tau/2);
fprintf('lateral pixel %.3f mm vs Nyquist for the two-way edge at fH: %.3f mm\n', dx, 1/(2*fH/(c*Fnum)));

%% Figure 1: image, k-space; PSF, k-space
toDB = @(A) max(-60, 20*log10(abs(A)/max(abs(A(:))) + 1e-300));
kzmax = kz_s(end);
figure('Name', 'k-space', 'Position', [60 60 900 720]);
subplot(2, 2, 1); imagesc(x_s, z_s, img); colormap(gray); axis image;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title('RF image: speckle + 3 points');
subplot(2, 2, 2); imagesc(kx_s, kz_s, toDB(fftshift(K))); set(gca, 'YDir', 'normal'); hold on;
plot([kx_s(1) kx_s(end)], [kz0 kz0], 'r--', [kx_s(1) kx_s(end)], -[kz0 kz0], 'r--');
plot([0 kzmax/(2*Fnum)], [0 kzmax], 'y--', [0 -kzmax/(2*Fnum)], [0 kzmax], 'y--');
plot([0 kzmax/(2*Fnum)], [0 -kzmax], 'y--', [0 -kzmax/(2*Fnum)], [0 -kzmax], 'y--'); hold off;
xlabel('k_x (cycles/mm)'); ylabel('k_z (cycles/mm)'); title('log|FFT2| (60 dB): red 2f_0/c, yellow bow-tie'); colorbar;
subplot(2, 2, 3); imagesc(x_s, z_s, psf); axis image;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title('single RF PSF');
subplot(2, 2, 4); imagesc(kx_s, kz_s, toDB(fftshift(H))); set(gca, 'YDir', 'normal'); hold on;
plot([kx_s(1) kx_s(end)], [kz0 kz0], 'r--', [kx_s(1) kx_s(end)], -[kz0 kz0], 'r--');
plot([-kedge kedge], [kz0 kz0], 'go', 'MarkerFaceColor', 'g'); hold off;
xlabel('k_x (cycles/mm)'); ylabel('k_z (cycles/mm)'); title('PSF k-space; green: \pm1/(\lambda_0 F#)'); colorbar;

%% Separable windows and the B-mode (envelope) PSF
Lw = 16; nw = 0:Lw-1;
wins = {ones(1, Lw), 0.5 - 0.5*cos(2*pi*nw/Lw)};  % rect, periodic Hann
wname = {'rect x rect', 'Hann x Hann'};
fprintf('\n--- separable 2D windows (%d x %d, 8x zero padded) ---\n', Lw, Lw);
figure('Name', 'separable windows and B-mode PSF', 'Position', [60 60 900 720]); colormap(gray);
for iw = 1:2
    w2 = wins{iw}.' * wins{iw};                   % separable: outer product
    W2 = abs(fftshift(fft2(w2, Nz, Nx))); W2 = W2/max(W2(:));
    ax = 20*log10(W2(Nz/2+1, Nx/2+1:end) + 1e-300);              % along kx at kz = 0
    k = 2; while ~(ax(k) < ax(k-1) && ax(k) <= ax(k+1)), k = k + 1; end
    [pslA, iA] = max(ax(k:end)); iA = iA + k - 1;               % first-sidelobe column
    diagv = 20*log10(W2(Nz/2+iA, Nx/2+iA));                      % product sidelobe on the diagonal
    fprintf('%-12s axis sidelobe %.2f dB, diagonal sidelobe %.2f dB (sum of the two 1D levels %.2f dB)\n', ...
            wname{iw}, pslA, diagv, 2*pslA);
    subplot(2, 2, iw); imagesc((-Nx/2:Nx/2-1)/Nx, (-Nz/2:Nz/2-1)/Nz, max(-80, 20*log10(W2 + 1e-300)));
    set(gca, 'YDir', 'normal'); axis image; xlim([-0.25 0.25]); ylim([-0.25 0.25]); colorbar;
    xlabel('k_x (cycles/sample)'); ylabel('k_z (cycles/sample)'); title(['|FFT2| ' wname{iw} ' (dB)']);
end
envp = abs(ifft2(fft2(abs(ha)) .* fft2(S1)));    % B-mode PSF (envelope), centered
subplot(2, 2, 3); imagesc(x_s, z_s, envp); axis image;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title('B-mode PSF (envelope)');
subplot(2, 2, 4); imagesc(kx_s, kz_s, toDB(fftshift(fft2(abs(ha))))); set(gca, 'YDir', 'normal'); colorbar;
xlabel('k_x (cycles/mm)'); ylabel('k_z (cycles/mm)'); title('envelope k-space: band moved to baseband');


%% Filtering in k-space: matched passband mask, then a one-sided kz mask for the envelope
Nimg  = randn(Nz, Nx); Nimg = Nimg*norm(img(:))/norm(Nimg(:));   % scaled to exactly 0 dB SNR
Kn    = fft2(Nimg);
Mpass = abs(H) > 0.01*max(abs(H(:)));        % passband: |H| above -40 dB (fft order)
fill  = mean(Mpass(:));                      % fraction of k-space kept
sigF = real(ifft2(K .* Mpass));             % filtered signal part
noiF = real(ifft2(Kn .* Mpass));            % filtered noise part
snr0 = round(100*10*log10(sum(img(:).^2)/sum(Nimg(:).^2)))/100 + 0;
snr1 = 10*log10(sum(sigF(:).^2)/sum(noiF(:).^2));
fprintf('\n--- k-space filtering (white noise at %.1f dB SNR) ---\n', snr0);
fprintf('passband mask |H| > -40 dB keeps %.1f %% of k-space; signal energy retained %.2f %%\n', 100*fill, 100*sum(sigF(:).^2)/sum(img(:).^2));
fprintf('SNR after mask %.1f dB: gain %.1f dB (prediction 10 log10(1/fill) = %.1f dB)\n', snr1, snr1 - snr0, 10*log10(1/fill));
% one-sided axial mask: keep kz > 0 (x2), drop kz < 0 -> analytic image along z, |.| = envelope
Mhalf = [ones(1, Nx); 2*ones(Nz/2-1, Nx); ones(1, Nx); zeros(Nz/2-1, Nx)];
anaI  = ifft2((K + Kn) .* Mpass .* Mhalf);
bmode = 20*log10(abs(anaI)/max(abs(anaI(:))) + 1e-300);
fprintf('one-sided kz mask: max |real(analytic) - filtered RF| = %.1e (the real part is the RF image)\n', ...
        max(max(abs(real(anaI) - (sigF + noiF)))));

figure('Name', 'k-space filtering', 'Position', [60 60 900 720]); colormap(gray);
subplot(2, 2, 1); imagesc(x_s, z_s, img + Nimg); axis image;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title(sprintf('RF image + white noise (%.0f dB SNR)', snr0));
subplot(2, 2, 2); imagesc(kx_s, kz_s, toDB(fftshift(K + Kn))); set(gca, 'YDir', 'normal'); hold on;
contour(kx_s, kz_s, double(fftshift(Mpass)), [0.5 0.5], 'g'); hold off; colorbar;
xlabel('k_x (cycles/mm)'); ylabel('k_z (cycles/mm)'); title('noisy k-space; green: passband mask');
subplot(2, 2, 3); imagesc(x_s, z_s, sigF + noiF); axis image;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title(sprintf('after the mask: SNR %.1f dB', snr1));
subplot(2, 2, 4); imagesc(x_s, z_s, max(-40, bmode)); axis image; colorbar;
xlabel('lateral x (mm)'); ylabel('axial z (mm)'); title('envelope via one-sided k_z mask (40 dB)');

%% m04_aperture_to_beam.m  -- Module 04: the beam pattern is the Fourier transform of the aperture
% Purpose : compute the far-field (or focal-plane) pattern of a linear array as the DTFT of
%           its apodization, compare rect vs Hann, one-way vs two-way (Tx = Rx), show that a
%           zero-padded FFT of the weights reproduces the direct far-field sum, and move
%           grating lobes with the pitch (d = lambda/2, lambda, 2 lambda).
% Change  : f0, N (elements), the pitch list, the steering angle th0_deg.
% Expect  : first sidelobe  rect -13.3 dB one-way / -26.5 dB two-way,
%                           Hann -31.5 dB one-way / -63 dB two-way;
%           -6 dB width (units of lambda/D in sin(theta)): rect 1.21 / 0.89, Hann 2.00 / 1.44;
%           FFT vs direct sum agree to ~1e-14; grating lobes at sin(th_g) = sin(th0) + m lambda/d
%           (d = lambda, 20 deg steering: -41.1 deg); element factor (w = 0.9 d) lowers that
%           grating lobe by ~4.4 dB one-way; a 60 % band fills the rect two-way first null
%           (-24 dB) and lowers the level at the narrowband sidelobe position to -28.8 dB.
% Model   : sections 1-3 monochromatic, point-like elements (no element directivity), far field
%           or focal plane; section 4 adds a rigid-baffle element factor sinc(w sin(theta)/lambda);
%           section 5 a Gaussian two-way spectrum (time-integrated energy pattern).
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% Parameters
c   = 1.540;            % mm/us
f0  = 5;                % MHz
lam = c/f0;             % wavelength (mm) = 0.308
N   = 64;               % elements
d   = lam/2;            % pitch for the apodization study (mm)
n   = 0:N-1;
xn  = (n - (N-1)/2)*d;  % element positions (mm), centered
D   = N*d;              % aperture length (mm)
wR  = ones(1, N);
wH  = 0.5 - 0.5*cos(2*pi*(n + 0.5)/N);       % symmetric Hann sampled at element centers
AF  = @(w, x, u) abs( exp(1j*2*pi/lam * u(:) * x) * w(:) ).' / sum(w);   % |sum w e^{jk x u}|, u = sin(th) - sin(th0)

%% 1. First sidelobe and -6 dB width: rect vs Hann, one-way vs two-way
uf = linspace(0, 6*lam/D, 6001);              % fine grid from the main-lobe peak outward
fprintf('--- apodization (N = %d, d = lambda/2, D = %.2f mm) ---\n', N, D);
fprintf('%-10s %-8s %12s %16s %14s\n', 'aperture', 'pattern', '1st SL (dB)', '-6 dB width', 'first null');
labels = {'rect', 'Hann'}; Ws = {wR, wH};
theorySL = [-13.26 -26.52; -31.47 -62.93]; theoryW6 = [1.207 0.886; 2.000 1.441];
for ia = 1:2
    a1 = AF(Ws{ia}, xn, uf);
    for way = 1:2
        a = a1.^way; dBv = 20*log10(a + 1e-300);
        k = 2; while ~(dBv(k) < dBv(k-1) && dBv(k) <= dBv(k+1)), k = k + 1; end      % first null
        j = k; while ~(dBv(j) > dBv(j-1) && dBv(j) >= dBv(j+1)), j = j + 1; end      % first sidelobe
        i6 = find(dBv < -6.0206, 1);
        u6 = uf(i6-1) + (-6.0206 - dBv(i6-1))*(uf(i6) - uf(i6-1))/(dBv(i6) - dBv(i6-1));
        fprintf('%-10s %-8s %7.2f (%6.2f) %9.3f (%5.3f) lam/D %8.3f lam/D\n', labels{ia}, char('one-way' * (way == 1) + 'two-way' * (way == 2)), ...
                dBv(j), theorySL(ia, way), 2*u6*D/lam, theoryW6(ia, way), uf(k)*D/lam);
    end
end
fprintf('(lateral -6 dB width in the focal plane = width in lam/D times lambda F#: rect two-way 0.886 lambda F#)\n');
Fn = 2;                                       % example F-number
fprintf('at F# = %d, %.0f MHz (lambda F# = %.3f mm): rect -6 dB width %.3f mm one-way, %.3f mm two-way\n', ...
        Fn, f0, lam*Fn, 1.207*lam*Fn, 0.886*lam*Fn);
fprintf('Hann/rect -6 dB width ratio: %.2f one-way, %.2f two-way\n', 2.000/1.207, 1.441/0.886);

%% 2. Beam = FFT of the aperture: zero-padded FFT vs direct far-field sum
L = 4096;
Wf = abs(fft(wH, L)) / sum(wH);                % DTFT samples at psi = -2 pi k/L
psi = -2*pi*(0:L-1)/L; psi = mod(psi + pi, 2*pi) - pi;     % wrap to [-pi, pi)
uk = psi*lam/(2*pi*d);                         % psi = k d u  ->  u = sin(theta)
vis = abs(uk) <= 1;
direct = AF(wH, xn, uk(vis));
fprintf('\n--- FFT of the weights vs direct far-field sum (Hann, L = %d) ---\n', L);
fprintf('max | |FFT| - |direct sum| | over %d visible-region samples = %.2e\n', nnz(vis), max(abs(Wf(vis) - direct)));
fprintf('sin(theta) per FFT bin = lambda/(L d) = %.5f\n', lam/(L*d));

%% 3. Grating lobes vs pitch
th0_deg = 20;
pitches = [0.5 1 2]*lam;
uu = linspace(-1, 1, 20001); thd = asind(uu);
fprintf('\n--- grating lobes, N = %d, steering 0 and %d deg ---\n', N, th0_deg);
figure('Name', 'grating lobes', 'Position', [80 80 760 640]);
for ip = 1:3
    dp = pitches(ip); xp = (n - (N-1)/2)*dp;
    for st = [0 th0_deg]
        s0 = sind(st);
        a = AF(wR, xp, uu - s0); adB = 20*log10(a + 1e-300);
        pk = find([adB(1) > adB(2), adB(2:end-1) > adB(1:end-2) & adB(2:end-1) >= adB(3:end), adB(end) > adB(end-1)] & adB > -1.5);
        meas = thd(pk); meas = meas(abs(meas - st) > 0.5);                          % drop the main lobe
        m = -4:4; m = m(m ~= 0); sg = s0 + m*lam/dp; pred = asind(sg(abs(sg) <= 1 + 1e-12));
        sp = sprintf(' %.1f', sort(pred)); if isempty(pred), sp = ' none'; end
        sm = sprintf(' %.1f', sort(meas)); if isempty(meas), sm = ' none'; end
        fprintf('d = %.2f lambda (%.3f mm), steer %2d deg: predicted [%s ] deg, measured [%s ] deg\n', dp/lam, dp, st, sp, sm);
        if st == th0_deg
            subplot(3, 1, ip); plot(thd, max(-60, adB), 'b'); hold on;
            for g = pred, plot([g g], [-60 3], 'r--'); end
            plot([st st], [-60 3], 'k:'); hold off; ylim([-60 3]); xlim([-90 90]); grid on;
            xlabel('angle \theta (deg)'); ylabel('|AF| (dB)');
            title(sprintf('d = %.1f \\lambda, steered %d deg: red = predicted grating lobes', dp/lam, st));
        end
    end
end
fprintf('d <= lambda/2: no grating lobe for any steering; d <= lambda/(1+sin(th_max)) for |th0| <= th_max\n');
fprintf('for th_max = %d deg: d <= %.3f mm (%.3f lambda)\n', th0_deg, lam/(1 + sind(th0_deg)), 1/(1 + sind(th0_deg)));

%% 4. Element directivity attenuates grating lobes (d = lambda, steered 20 deg)
dp = lam; xp = (n - (N-1)/2)*dp; s0 = sind(th0_deg);
we = 0.9*dp;                                       % element width (mm): 10 % kerf
sincf = @(v) (sin(pi*v) + (v == 0)) ./ (pi*v + (v == 0));
EF = @(uabs) abs(sincf(we*uabs/lam));
thg = asind(s0 - lam/dp);                          % predicted grating lobe (deg)
aM = AF(wR, xp, 0)        * EF(s0);                % main lobe (u = 0 at theta = th0)
aG = AF(wR, xp, sind(thg) - s0) * EF(sind(thg));   % grating lobe
fprintf('\n--- element directivity (d = lambda, w = 0.9 d, steer %d deg) ---\n', th0_deg);
fprintf('grating lobe at %.1f deg; element factor there %.2f dB, at the main lobe %.2f dB\n', thg, 20*log10(EF(sind(thg))), 20*log10(EF(s0)));
fprintf('grating lobe re main lobe: %.2f dB one-way, %.2f dB two-way (0 dB without the element factor)\n', 20*log10(aG/aM), 40*log10(aG/aM));

%% 5. Broadband: a 60 % two-way band smears sidelobes and fills nulls (energy pattern)
FBW = 0.60; sf = FBW/(2*sqrt(2*log(2)));           % Gaussian amplitude-spectrum std (units of f0)
fr  = linspace(1 - 4*sf, 1 + 4*sf, 41);            % frequencies (units of f0)
Pw  = exp(-(fr - 1).^2/sf^2);                      % |P(f)|^2 weights
ue  = [0 1 1.43 2 2.46]*lam/D;                     % sin(theta): peak, rect null, rect SL, Hann null, Hann SL
AFl = @(w, x, u, lm) abs( exp(1j*2*pi/lm * u(:) * x) * w(:) ).' / sum(w);
fprintf('\n--- broadband two-way energy pattern, FBW = %.0f %% (levels in dB re peak) ---\n', 100*FBW);
fprintf('%-6s %-11s %8s %8s %8s %8s\n', 'apod', 'model', 'u=1', 'u=1.43', 'u=2', 'u=2.46');
ub2 = linspace(-8*lam/D, 8*lam/D, 4001); Ebb = zeros(2, numel(ub2));
for ia = 1:2
    Eu = zeros(1, numel(ue)); Eb = zeros(1, numel(ub2));
    for kf = 1:numel(fr)
        Eu = Eu + Pw(kf)*AFl(Ws{ia}, xn, ue,  lam/fr(kf)).^4;
        Eb = Eb + Pw(kf)*AFl(Ws{ia}, xn, ub2, lam/fr(kf)).^4;
    end
    Ebb(ia, :) = 10*log10(Eb/max(Eb));
    nb = 40*log10(AF(Ws{ia}, xn, ue) + 1e-300);
    fprintf('%-6s %-11s %8.1f %8.1f %8.1f %8.1f\n', labels{ia}, 'narrowband', max(nb(2:end), -99));
    fprintf('%-6s %-11s %8.1f %8.1f %8.1f %8.1f\n', labels{ia}, 'broadband', 10*log10(Eu(2:end)/Eu(1)));
end
fprintf('(u in units of lambda0/D; -99 marks an exact null)\n');

%% Figures: apodization and FFT-vs-direct
figure('Name', 'aperture to beam', 'Position', [80 80 760 640]);
subplot(2, 2, 1); plot(xn, wR, 'bo', xn, wH, 'r.-', 'MarkerSize', 8); ylim([0 1.1]);
xlabel('element position (mm)'); ylabel('weight'); title('aperture: rect (blue), Hann (red)');
ub = linspace(-8*lam/D, 8*lam/D, 4001);
subplot(2, 2, 2);
plot(ub*D/lam, 20*log10(AF(wR, xn, ub)), 'b', ub*D/lam, 20*log10(AF(wH, xn, ub)), 'r'); ylim([-80 3]); grid on;
xlabel('sin\theta in units of \lambda/D'); ylabel('dB'); title('one-way: -13.3 vs -31.5 dB');
subplot(2, 2, 3);
plot(ub*D/lam, 40*log10(AF(wR, xn, ub)), 'b', ub*D/lam, 40*log10(AF(wH, xn, ub)), 'r'); hold on;
plot(ub2*D/lam, Ebb(1, :), 'b--', ub2*D/lam, Ebb(2, :), 'r--', 'LineWidth', 1.2); hold off; ylim([-100 3]); grid on;
xlabel('sin\theta in units of \lambda/D'); ylabel('dB'); title('two-way: narrowband (solid), 60 % band (dashed)');
subplot(2, 2, 4);
ud = linspace(-1, 1, 4001);
plot(ud, 20*log10(AF(wH, xn, ud)), 'k'); hold on;
iv = find(vis); iv = iv(1:4:end);
plot(uk(iv), 20*log10(Wf(iv) + 1e-300), 'ro', 'MarkerSize', 3); hold off; ylim([-100 3]); xlim([-0.3 0.3]); grid on;
xlabel('sin\theta'); ylabel('dB'); title('Hann: direct sum (black) vs FFT bins (red)');

%% m04_rf_convolution.m  -- Module 04: RF line = pulse * reflectivity, spectra multiply, axial resolution
% Purpose : model an RF line as a Gaussian pulse convolved with a reflectivity profile,
%           verify that spectra multiply (with enough zero padding), measure the -6 dB
%           pulse length and axial resolution, resolve two reflectors, and compute the
%           envelope with an FFT-based analytic signal (no hilbert()).
% Change  : f0, fractional bandwidth FBW, the reflector spacings, fs.
% Expect  : tau(-6 dB) = 4 ln2/(pi B) = 0.294 us and dz = c tau/2 = 0.227 mm at 5 MHz, 60 %;
%           FFT and direct convolution agree to ~1e-13 when Nfft >= L+M-1;
%           Gabor product sigma_t' sigma_f' = 1/(4 pi) = 0.0796.
% Base MATLAB only; runs in MATLAB R2019b+ and GNU Octave 8.
clear; close all; rng(1);

%% Parameters
c   = 1.540;            % speed of sound (mm/us)
fs  = 100;              % sampling rate (MHz)
f0  = 5;                % center frequency (MHz)
FBW = 0.60;             % fractional bandwidth at -6 dB (amplitude)
B   = FBW*f0;                         % -6 dB bandwidth (MHz)
sig_t = sqrt(2*log(2))/(pi*B);        % envelope std (us): B = 2 sqrt(2 ln2) sigma_f, sigma_f = 1/(2 pi sigma_t)
tau_th = 4*log(2)/(pi*B);             % -6 dB pulse duration (us)
dz_th  = c*tau_th/2;                  % axial resolution (mm)

%% Pulse (two-way impulse response)
tp = (-ceil(4.5*sig_t*fs):ceil(4.5*sig_t*fs))/fs;     % time axis of the pulse (us)
p  = exp(-tp.^2/(2*sig_t^2)) .* cos(2*pi*f0*tp);
M  = numel(p);

%% Analytic signal by FFT (replaces hilbert): keep DC, double positive bins, zero negative bins
analytic = @(x, Nf) ifft( fft(x, Nf) .* [1, 2*ones(1, Nf/2-1), 1, zeros(1, Nf/2-1)] );

%% Measured pulse length from its envelope
Nf  = 2^nextpow2(8*M);
env_p = abs(analytic(p, Nf)); env_p = env_p(1:M);
above = find(env_p >= max(env_p)/2);
% sub-sample crossing by linear interpolation at both edges
i1 = above(1); i2 = above(end);
t1 = tp(i1-1) + (max(env_p)/2 - env_p(i1-1))/(env_p(i1) - env_p(i1-1))/fs;
t2 = tp(i2) + (env_p(i2) - max(env_p)/2)/(env_p(i2) - env_p(i2+1))/fs;
tau_meas = t2 - t1;
% Gabor product (rms widths of |envelope|^2 and |P(f)|^2, one-sided spectrum)
e2 = env_p.^2; mt = sum(tp.*e2)/sum(e2); st2 = sqrt(sum((tp-mt).^2.*e2)/sum(e2));
Pf = abs(fft(p, Nf)); fax = (0:Nf-1)*fs/Nf; pos = 1:Nf/2;
q2 = Pf(pos).^2; mf = sum(fax(pos).*q2)/sum(q2); sf2 = sqrt(sum((fax(pos)-mf).^2.*q2)/sum(q2));
fprintf('--- pulse: f0 = %.1f MHz, B(-6 dB) = %.2f MHz ---\n', f0, B);
fprintf('tau(-6 dB): theory 4ln2/(pi B) = %.4f us, measured from FFT envelope = %.4f us\n', tau_th, tau_meas);
fprintf('axial resolution c*tau/2: theory %.4f mm, measured %.4f mm\n', dz_th, c*tau_meas/2);
fprintf('Gabor product sigma_t''*sigma_f'' = %.5f (1/(4 pi) = %.5f)\n', st2*sf2, 1/(4*pi));

%% Reflectivity: speckle line (z = 0..6 mm) and the RF line by FFT convolution
zmax = 6; Ns = round(2*zmax/c*fs);              % samples covering 0..zmax (round trip)
rho = zeros(1, Ns); hit = rand(1, Ns) < 0.15;    % ~20 scatterers per mm (sample = 7.7 um)
rho(hit) = randn(1, nnz(hit));
L = Ns;
Nlin = L + M - 1;                                % linear convolution length
Nfft = 2^nextpow2(Nlin);
r_fft = real(ifft(fft(rho, Nfft) .* fft(p, Nfft)));  r_fft = r_fft(1:Nlin);
r_dir = conv(rho, p);
Nshort = L;                                      % too short (< L+M-1): circular wrap-around
r_circ = real(ifft(fft(rho, Nshort) .* fft(p, Nshort)));
fprintf('\n--- convolution ---\n');
fprintf('L = %d, M = %d, L+M-1 = %d, Nfft = %d: max|FFT conv - conv| = %.2e\n', L, M, Nlin, Nfft, max(abs(r_fft - r_dir)));
fprintf('Nfft = %d (< L+M-1): max error on first %d samples = %.2e (wrap-around)\n', Nshort, M, max(abs(r_circ(1:M) - r_dir(1:M))));
% spectra multiply
R = fft(r_dir, Nfft); S = fft(rho, Nfft); P = fft(p, Nfft);
fprintf('max|R - P.S| / max|R| = %.2e  (convolution in time = multiplication in frequency)\n', max(abs(R - P.*S))/max(abs(R)));

%% Two reflectors: axial resolution test
% The envelope of two echoes depends on their relative carrier phase phi = 4 pi f0 s / c:
% in phase (phi = 0) the envelopes add, in quadrature (phi = pi/2) the intensities add
% (= the average over random phase), in antiphase they cancel at the midpoint.
spac = [0.15 0.20 dz_th 0.25 0.30 0.40];           % spacings (mm)
fprintf('\n--- two equal reflectors at 3 mm: envelope dip between them (dB) ---\n');
fprintf('%10s %10s %12s %10s %10s %12s\n', 'spacing mm', 'spacing/dz', 'phi (deg)', 'dip RF', 'dip phi=0', 'dip phi=90');
zc = 3; Nr = 2^nextpow2(Ns + M);
gEnv = @(t) exp(-t.^2/(2*sig_t^2));               % Gaussian envelope
tt = linspace(-1, 1, 4001);                        % fine time axis (us) for the envelope-only cases
dipOf = @(e, a, b) 20*log10(min(e(a:b))/min(e(a), e(b)));
figure('Name', 'two reflectors', 'Position', [100 100 760 560]);
for is = 1:numel(spac)
    s = spac(is);
    rr = zeros(1, Ns);
    k1 = round(2*(zc - s/2)/c*fs) + 1; k2 = round(2*(zc + s/2)/c*fs) + 1;
    rr(k1) = 1; rr(k2) = rr(k2) + 1;
    rf = conv(rr, p);                                % pulse centered: shift by (M-1)/2
    rf = rf((M+1)/2 : (M+1)/2 + Ns - 1);
    en = abs(analytic(rf, Nr)); en = en(1:Ns);
    dip = dipOf(en, k1, k2);
    phi = mod(2*pi*f0*(k2 - k1)/fs, 2*pi);          % actual relative phase after sample rounding
    dt = s/c*2/2;                                    % half the echo separation (us)
    e1 = gEnv(tt + dt); e2 = gEnv(tt - dt);
    [~, a] = min(abs(tt + dt)); [~, b] = min(abs(tt - dt));
    d0 = dipOf(e1 + e2, a, b);                       % in phase
    d90 = dipOf(sqrt(e1.^2 + e2.^2), a, b);          % quadrature = intensities add
    fprintf('%10.3f %10.2f %12.0f %10.1f %10.1f %12.1f\n', s, s/dz_th, phi*180/pi, dip, d0, d90);
    if is == 2 || is == 5
        zz = (0:Ns-1)/fs*c/2;
        subplot(2, 1, 1 + (is == 5));
        plot(zz, rf, 'Color', [0.2 0.4 0.8]); hold on; plot(zz, en, 'r', 'LineWidth', 1.4);
        plot([zc-s/2 zc+s/2], [0 0], 'kv', 'MarkerFaceColor', 'k'); hold off;
        xlim([zc-1 zc+1]); xlabel('depth z = ct/2 (mm)'); ylabel('RF, envelope');
        title(sprintf('s = %.2f mm = %.2f dz, phase %.0f deg: dip %.1f dB', s, s/dz_th, phi*180/pi, dip));
    end
end
fprintf('At spacing = c tau/2 (%.3f mm) the quadrature (random-phase average) dip is -3.0 dB.\n', dz_th);

%% Figure: pulse, spectra, speckle RF line with FFT envelope
figure('Name', 'pulse and RF line', 'Position', [100 100 760 560]);
subplot(2, 2, 1);
plot(tp, p, 'b', tp, env_p, 'r', tp, -env_p, 'r'); grid on;
xlim([-0.6 0.6]); xlabel('time (us)'); ylabel('p(t)'); title(sprintf('pulse, tau = %.3f us', tau_meas));
subplot(2, 2, 2);
fa = (0:Nfft/2)*fs/Nfft;
Sn = abs(S(1:Nfft/2+1)); Pn = abs(P(1:Nfft/2+1)); Rn = abs(R(1:Nfft/2+1));
plot(fa, Sn/max(Sn), 'Color', [0.6 0.6 0.6]); hold on;
plot(fa, Pn/max(Pn), 'b', 'LineWidth', 1.5); plot(fa, Rn/max(Rn), 'r'); hold off;
xlim([0 12]); xlabel('frequency (MHz)'); ylabel('|.| (norm.)'); title('|S| gray, |P| blue, |PS| red');
subplot(2, 1, 2);
rfs = r_dir((M+1)/2 : (M+1)/2 + Ns - 1);
ens = abs(analytic(rfs, 2^nextpow2(2*Ns))); ens = ens(1:Ns);
zz = (0:Ns-1)/fs*c/2;
plot(zz, rfs, 'Color', [0.2 0.4 0.8]); hold on; plot(zz, ens, 'r', 'LineWidth', 1.2); hold off;
xlabel('depth z = ct/2 (mm)'); ylabel('amplitude'); title('speckle RF line and FFT-analytic envelope');