Eigendecomposition & SVD

%% m06_svd_cyst.m  --  Module 06: low-rank approximation of a B-mode cyst phantom
% Purpose : Build a synthetic B-mode image (random scatterers * 2-D RF PSF,
%           anechoic cyst, point targets, layers), detect the envelope with an
%           FFT-based analytic signal, log-compress, then study rank-k
%           approximations: images at ranks 1, 2, 5, 20, 60, the singular-value
%           spectrum, the Eckart-Young identities ||B - B_k||_F^2 = sum_{i>k} sigma_i^2
%           and ||B - B_k||_2 = sigma_{k+1}, a rival rank-k approximation that
%           must lose, and cyst contrast / CNR vs k.
% Change  : f0, PSF widths, cyst radius, dynamic range DR, the list of ranks ks.
% Expect  : both Eckart-Young identities hold (Frobenius to ~1e-14, spectral norm to
%           the printed digits); the rival rank-k approximation loses by 1.5-2.2x;
%           rank 1 already holds 97.3 % of the energy (the -30 dB mean level), so
%           judge the images, not the percentage. Cyst contrast about -21 dB at full
%           rank, capped by the 60 dB display floor (the cyst interior is clipped),
%           and only -6 dB at k = 1; CNR peaks at k = 2 (9.4 vs 3.8 at full rank)
%           because low rank averages speckle, at the price of separable streaks.
% Runtime : ~4 s in Octave.
clear; close all; rng(1);

%% 1. Geometry and RF point-spread function
c0 = 1540; f0 = 5e6; fs = 40e6;              % m/s, Hz, Hz
dz = c0/(2*fs);                              % axial sample spacing (m), 19.25 um
dx = 0.15e-3;                                % lateral line spacing (m)
z = (10e-3:dz:40e-3)'; x = (-19.2e-3:dx:19.2e-3-dx);  % 1559 x 256 grid
Nz = numel(z); Nx = numel(x);
sz = 0.12e-3;                                % axial envelope std (m), FWHM 0.28 mm
sx = 0.30e-3;                                % lateral beam std (m), FWHM 0.71 mm
kz = (-ceil(4*sz/dz):ceil(4*sz/dz))'*dz;  kx = (-ceil(4*sx/dx):ceil(4*sx/dx))*dx;
psf = exp(-kz.^2/(2*sz^2)) .* cos(2*pi*f0*2*kz/c0) * exp(-kx.^2/(2*sx^2));  % RF PSF (two-way phase)

%% 2. Scatterer field: echogenicity map x Gaussian random amplitudes
[X, Z] = meshgrid(x, z);
echo = ones(Nz, Nx);
echo(Z > 14e-3 & Z < 15.5e-3) = 6;           % bright layer (+7.8 dB)
echo(Z > 31e-3 & Z < 33e-3)   = 0.3;         % dark layer (-5.2 dB)
cyst = (X - 2e-3).^2 + (Z - 24e-3).^2 < (4e-3)^2;
echo(cyst) = 0;                              % anechoic cyst, radius 4 mm
incl = ((X + 9e-3)/3e-3).^2 + ((Z - 24e-3)/2.5e-3).^2 < 1;
echo(incl) = 3;                              % hyperechoic inclusion (+4.8 dB)
scat = sqrt(echo) .* randn(Nz, Nx);          % one random scatterer per sample
pts = [-12e-3 20e-3; -4e-3 20e-3; 4e-3 20e-3; 12e-3 20e-3; 10e-3 36e-3];
for p = 1:size(pts,1)
    [~, iz] = min(abs(z - pts(p,2))); [~, ix] = min(abs(x - pts(p,1)));
    scat(iz, ix) = scat(iz, ix) + 400;        % point targets (wires)
end
% 2-D convolution via zero-padded FFTs (base MATLAB conv2 would also work)
Pz = Nz + numel(kz) - 1; Px = Nx + numel(kx) - 1;
rf = real(ifft2(fft2(scat, Pz, Px) .* fft2(psf, Pz, Px)));
oz = (numel(kz)-1)/2; ox = (numel(kx)-1)/2;
rf = rf(oz+(1:Nz), ox+(1:Nx));
rf = rf + 0.002*std(rf(:))*randn(Nz, Nx);    % electronic noise, -54 dB re mean speckle

%% 3. Envelope (analytic signal via FFT, no toolbox hilbert) and log compression
N = Nz; h = zeros(N,1); h(1) = 1;
if mod(N,2) == 0, h(N/2+1) = 1; h(2:N/2) = 2; else, h(2:(N+1)/2) = 2; end
env = abs(ifft(fft(rf) .* repmat(h, 1, Nx)));
env = env(1:4:end, :); zd = z(1:4:end);      % decimate axially to 77 um pixels
cystd = cyst(1:4:end, :);
bgd = ~cystd & abs(Z(1:4:end,:) - 24e-3) < 4e-3 & abs(X(1:4:end,:) - 2e-3) > 6e-3 & abs(X(1:4:end,:) + 9e-3) > 5e-3;
DR = 60;                                     % dynamic range (dB)
B = 20*log10(env/max(env(:)));  B = max(B, -DR);   % log-compressed image (dB), A in the text
[m, n] = size(B);
fprintf('B-mode image %d x %d pixels (%.0f um axial x %.0f um lateral), DR = %d dB\n', m, n, 4*dz*1e6, dx*1e6, DR);

%% 4. SVD, spectrum, Eckart-Young check, cyst metrics vs rank
[U, S, V] = svd(B, 'econ'); s = diag(S);
inC = cystd & ((X(1:4:end,:) - 2e-3).^2 + (Z(1:4:end,:) - 24e-3).^2 < (3e-3)^2);  % cyst core
ks = [1 2 5 10 20 40 80 n];
fprintf('\nEckart-Young: ||B-Bk||_F^2 = sum_{i>k} sigma_i^2 and ||B-Bk||_2 = sigma_{k+1}\n');
fprintf('   k   energy%%   ||B-Bk||_F^2    sum_{i>k} s^2   rel.diff   ||B-Bk||_2  sigma_k+1  contrast(dB)  CNR\n');
con = zeros(size(ks)); cnr = con;
for j = 1:numel(ks)
    k = ks(j);
    Bk = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';
    e2 = norm(B - Bk, 'fro')^2; tail = sum(s(k+1:end).^2);
    if k < n, e2s = norm(B - Bk, 2); sk1 = s(k+1); else, e2s = 0; sk1 = 0; end
    mc = mean(Bk(inC)); mb = mean(Bk(bgd));
    con(j) = mc - mb;                        % contrast on the dB image
    cnr(j) = abs(mc - mb)/sqrt(var(Bk(inC)) + var(Bk(bgd)));
    if tail > 0, rd = abs(e2 - tail)/tail; else, rd = 0; end
    fprintf('%4d   %7.3f   %13.6e   %13.6e   %8.1e   %9.3f  %9.3f   %7.1f     %5.2f\n', k, 100*sum(s(1:k).^2)/sum(s.^2), e2, tail, rd, e2s, sk1, con(j), cnr(j));
end
fprintf('sigma_1 / sigma_n = %.3g;  ranks for 99%% energy: %d;  for 99.9%%: %d\n', s(1)/s(end), ...
    find(cumsum(s.^2)/sum(s.^2) >= 0.99, 1), find(cumsum(s.^2)/sum(s.^2) >= 0.999, 1));
fprintf('(Energy is dominated by the mean level: judge the residual image, not the percentage.)\n');
% Optimality: any other rank-k matrix does worse. Rival: project B onto the span of
% its k highest-energy columns (column-subset selection, a common "cheap" low rank).
[~, ord] = sort(sum(B.^2, 1), 'descend');
fprintf('\nRival rank-k approximation (span of the k strongest columns) vs truncated SVD, Frobenius error ratio:\n');
for k = [1 5 20 60]
    Q = orth(B(:, ord(1:k))); Bc = Q*(Q'*B);
    fprintf('  k = %2d: ||B-Bc||_F / ||B-Bk||_F = %.3f  (>= 1 by Eckart-Young)\n', k, norm(B - Bc, 'fro')/sqrt(sum(s(k+1:end).^2)));
end

%% 5. Figures
figure('Position', [50 50 1000 640]);
show = [n 1 2 5 20 60]; ttl = {'original (full rank)', 'rank 1', 'rank 2', 'rank 5', 'rank 20', 'rank 60'};
for j = 1:numel(show)
    k = show(j); Bk = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';
    subplot(2,3,j); imagesc(x*1e3, zd*1e3, Bk, [-DR 0]); colormap(gray); axis image;
    title(ttl{j}); xlabel('lateral (mm)'); if mod(j,3) == 1, ylabel('depth (mm)'); end
end
figure('Position', [50 50 1000 360]);
subplot(1,2,1); semilogy(s/s(1), 'k.-'); grid on; hold on;
semilogy(ks, s(ks)/s(1), 'ro');
xlabel('index i'); ylabel('\sigma_i / \sigma_1'); title('Singular values of the log-compressed image');
subplot(1,2,2); [ax, h1, h2] = plotyy(ks, con, ks, cnr); set(ax, 'XScale', 'log');
set(h1, 'Marker', 'o'); set(h2, 'Marker', 's'); grid on;
xlabel('rank k'); ylabel(ax(1), 'cyst contrast (dB)'); ylabel(ax(2), 'CNR');
title('Cyst metrics of the rank-k image');

%% m06_svd_clutter_filter.m  --  Module 06: SVD clutter filter on a simulated ultrafast ensemble
% Purpose : Simulate an ultrafast IQ ensemble (pixels x frames): moving tissue
%           speckle, a slow-flow vessel, a faster vessel, a few microbubble tracks,
%           white noise. Build the Casorati matrix, take its SVD, look at the
%           singular-value curve (with the Marchenko-Pastur noise edge), the
%           Doppler spectra of the temporal singular vectors, and power Doppler
%           for: raw data, an SVD band, and three FIXED temporal filters (first
%           difference, 2nd-order Butterworth IIR high-pass run forward-backward,
%           polynomial-regression filter). Then sweep tissue velocity and watch
%           the best SVD band move while every fixed temporal filter collapses.
%           Last, check the fast "Gram route" (eig of the n x n matrix S'S) in
%           double and single precision.
% Change  : vt_main (tissue peak axial velocity), the sweep list vts, PRF, n (frames),
%           blood_dB, noise_dB, the cutoffs fcs and polynomial degrees pdeg.
% Expect  : (rng(1)) at vt = 0 every filter works (SVD 24.7 dB, IIR 23.1 dB contrast).
%           At vt = 8 mm/s (52 Hz peak tissue Doppler) the SVD band 9-46 gives
%           15.2 dB vessel/background contrast; IIR, regression and first difference
%           all stay near 0 dB. The oracle k_lo rises 2 -> 5 -> 9 -> 11 -> 15 for
%           vt = 0, 4, 8, 16, 24 mm/s; k_hi from the MP edge is 45-47. The Gram
%           route matches svd() to ~1e-11 in double; in single it misplaces
%           sigma_k by > 1 % about 66-69 dB below sigma_1.
%           These are oracle (ground-truth-mask) choices: in vivo you have no masks.
% Runtime : ~20 s in Octave (five simulated ensembles; the simulation dominates).
clear; close all; rng(1);

%% 1. Parameters
c0 = 1540; f0 = 5e6; lam = c0/f0*1e3;        % wavelength (mm), 0.308 mm
k2 = 4*pi/lam;                               % two-way wavenumber (rad/mm)
PRF = 1000; n = 128;                         % frame rate (Hz), frames -> 128 ms
NX = 48; NZ = 48; dx = 0.1;                  % pixels, pixel pitch (mm)
xg = (0:NX-1)*dx; zg = (0:NZ-1)*dx; m = NX*NZ;
sx = 0.17; sz = 0.10;                        % PSF std lateral / axial (mm)
blood_dB = -30;                              % blood scatterer amplitude re tissue scatterer
noise_dB = -50;                              % noise power re mean tissue pixel power
fm = 3;                                      % tissue motion frequency (Hz), e.g. breathing/pulsation
vt_main = 8;                                 % tissue peak axial velocity for the main figure (mm/s)
vts = [0 4 8 16 24];                         % velocity sweep (mm/s)
fcs = [25 50 100 150];                       % high-pass cutoffs to try (Hz)
pdeg = [1 2 4 8];                            % polynomial-regression degrees to try

% Scatterers are drawn once and reused for every velocity, so only motion changes.
L = NX*dx; Ntis = 4000;
tis = [rand(Ntis,1)*(L+2)-1, rand(Ntis,1)*(L+2)-1, (randn(Ntis,1)+1i*randn(Ntis,1))/sqrt(2)];
ves = [1.2 2.2 10 0.30  5;                    % [x0 z0 angle(deg) radius(mm) vmax(mm/s)]: slow vessel
       3.4 2.6 70 0.22 15];                   % faster, steep vessel
blood = [];                                   % [vessel, s0, r, amplitude]
for v = 1:size(ves,1)
    nb = round(500*ves(v,4)/0.3);
    blood = [blood; v*ones(nb,1), rand(nb,1)*8-4, (2*rand(nb,1)-1)*ves(v,4), ...
             10^(blood_dB/20)*(randn(nb,1)+1i*randn(nb,1))/sqrt(2)]; %#ok<AGROW>
end
% microbubbles: [xc zc radius(mm) speed(mm/s) start angle(rad)] on circular arcs
mb = [2.4 1.2 1.5 12 0.3; 1.6 3.6 1.2 20 2.0; 3.8 3.8 1.8 9 4.0; 0.8 1.0 1.0 16 5.5];
mbAmp = 1.0;                                  % 0 dB re a tissue scatterer

PSFx = @(x) exp(-(xg(:) - x(:).').^2/(2*sx^2));        % NX x Ns
PSFz = @(z) exp(-(zg(:) - z(:).').^2/(2*sz^2));        % NZ x Ns

%% 2. Simulate one ensemble per tissue velocity; each gives a Casorati matrix S (m x n)
S_all = cell(numel(vts),1); masks = [];
for iv = 1:numel(vts)
    vt = vts(iv); w = 2*pi*fm;
    S = zeros(m, n);
    for f = 1:n
        t = (f-1)/PRF;
        dz0 = vt/w*sin(w*t);                  % bulk axial displacement (mm), peak velocity vt
        % mild non-rigid motion: displacement grows 30 % across the field (depth-dependent)
        dzs = @(z) dz0*(1 + 0.3*(z - L/2)/L);
        dxs = 0.35*dz0;
        zt = tis(:,2) + dzs(tis(:,2)); xt = tis(:,1) + dxs;
        amp = tis(:,3).*exp(-1i*k2*zt);
        img = PSFx(xt)*(amp.*PSFz(zt).');       % NX x NZ
        % blood: parabolic profile, riding on the tissue motion
        a = ves(blood(:,1),3)*pi/180; ux = cos(a); uz = sin(a); R = ves(blood(:,1),4);
        vel = ves(blood(:,1),5).*(1 - (blood(:,3)./R).^2);
        sp = mod(blood(:,2) + vel*t + 4, 8) - 4;
        xb = ves(blood(:,1),1) + ux.*sp - uz.*blood(:,3);
        zb = ves(blood(:,1),2) + uz.*sp + ux.*blood(:,3);
        zb = zb + dzs(zb); xb = xb + dxs;
        img = img + PSFx(xb)*((blood(:,4).*exp(-1i*k2*zb)).*PSFz(zb).');
        % microbubbles on arcs
        ang = mb(:,5) + mb(:,4)./mb(:,3)*t;
        xm = mb(:,1) + mb(:,3).*cos(ang) + dxs; zm = mb(:,2) + mb(:,3).*sin(ang);
        zm = zm + dzs(zm);
        img = img + PSFx(xm)*((mbAmp*exp(-1i*k2*zm)).*PSFz(zm).');
        S(:, f) = img(:);                     % column = one frame (x-major, NX*NZ)
    end
    if iv == 1, p0 = mean(abs(S(:,1)).^2); end % mean tissue pixel power (static case)
    sn = sqrt(p0*10^(noise_dB/10));           % noise std per entry
    S = S + sn*(randn(m,n) + 1i*randn(m,n))/sqrt(2);
    S_all{iv} = S;
end

%% 3. Masks (ground truth, at the ensemble-mean position; only used for scoring)
[XX, ZZ] = ndgrid(xg, zg);                    % NX x NZ, matches img(:)
dn = inf(NX, NZ); dv = cell(1, size(ves,1));
for v = 1:size(ves,1)
    a = ves(v,3)*pi/180;
    dv{v} = abs(-(XX - ves(v,1))*sin(a) + (ZZ - ves(v,2))*cos(a))/ves(v,4);   % distance in radii
    dn = min(dn, dv{v});
end
dmb = inf(NX, NZ); tt = linspace(0, (n-1)/PRF, 60);
for b = 1:size(mb,1)
    ang = mb(b,5) + mb(b,4)/mb(b,3)*tt;
    for q = 1:numel(tt)
        dmb = min(dmb, hypot(XX - mb(b,1) - mb(b,3)*cos(ang(q)), ZZ - mb(b,2) - mb(b,3)*sin(ang(q))));
    end
end
vesM = dn(:) < 0.6;                           % both vessel cores
slowM = dv{1}(:) < 0.6; fastM = dv{2}(:) < 0.6;
bgM = dn(:) > 2.5 & dmb(:) > 0.5;             % background, away from vessels and bubble tracks
contrast = @(pd) 10*log10(mean(pd(vesM))/mean(pd(bgM)));
cSlow = @(pd) 10*log10(mean(pd(slowM))/mean(pd(bgM)));
cFast = @(pd) 10*log10(mean(pd(fastM))/mean(pd(bgM)));
c3 = @(pd) [contrast(pd), cSlow(pd), cFast(pd)];

%% 4. Filters
pdOf = @(Y) mean(abs(Y).^2, 2);               % power Doppler per pixel
fd = @(Y) Y(:,2:end) - Y(:,1:end-1);          % first difference: |H| = 2|sin(pi f/PRF)|
% 2nd-order Butterworth high-pass via the bilinear transform (no toolbox)
bwb = @(fc) [1, -2, 1] / (1 + sqrt(2)*tan(pi*fc/PRF) + tan(pi*fc/PRF)^2);
bwa = @(fc) [1, 2*(tan(pi*fc/PRF)^2 - 1), 1 - sqrt(2)*tan(pi*fc/PRF) + tan(pi*fc/PRF)^2] ...
            / (1 + sqrt(2)*tan(pi*fc/PRF) + tan(pi*fc/PRF)^2);
hp1 = @(Y, fc) filter(bwb(fc), bwa(fc), Y - Y(:,1), [], 2);   % one pass; subtract first sample to limit the step
rev = @(Y) Y(:, end:-1:1);
iir = @(Y, fc) rev(hp1(rev(hp1(Y, fc)), fc));                 % forward-backward: zero phase, |H|^2 (like filtfilt)
cols = @(Y, idx) Y(:, idx);
iirPD = @(Y, fc) mean(abs(cols(iir(Y, fc), 21:n-20)).^2, 2);  % drop 20-frame transients at both ends
% Polynomial regression filter: remove each pixel's projection onto polynomials of
% degree 0..p in slow time. Like the SVD filter it is a projection, S - S*Q*Q',
% but onto a FIXED temporal subspace chosen before seeing the data.
tn = linspace(-1, 1, n)';                     % normalized slow time
regPD = @(Y, p) pdOf(Y - (Y*orth(tn.^(0:p)))*orth(tn.^(0:p))');

%% 5. Main case: singular values, Doppler of temporal singular vectors, PD images
iv = find(vts == vt_main); S = S_all{iv};
[U, Sg, V] = svd(S, 'econ'); s = diag(Sg);
sn = sqrt(p0*10^(noise_dB/10));
edgeMP = sn*(sqrt(m) + sqrt(n));              % Marchenko-Pastur upper edge of the noise bulk
k_hi = find(s > 1.05*edgeMP, 1, 'last');      % keep everything above the noise bulk
fprintf('Casorati %d x %d, vt = %g mm/s, tissue peak Doppler 2vt/lambda = %.1f Hz\n', m, n, vt_main, 2*vt_main/lam);
fprintf('sigma_k (dB re sigma_1), k = 1..10: %s\n', sprintf('%.0f ', 20*log10(s(1:10)/s(1))));
fprintf('MP noise edge = %.1f dB re sigma_1  ->  k_hi = %d (measured last sigma = %.1f dB)\n', ...
    20*log10(edgeMP/s(1)), k_hi, 20*log10(s(end)/s(1)));
NF = 256; fax = ((0:NF-1) - NF/2)*PRF/NF;
Pk = abs(fftshift(fft(conj(V), NF, 1), 1)).^2;           % NF x n, spectrum of each temporal vector
fbar = (abs(fax)*Pk)./sum(Pk, 1);                         % mean |Doppler| per rank (Hz)
fprintf('mean |Doppler| of temporal singular vectors, k = 1..12 (Hz): %s\n', sprintf('%.0f ', fbar(1:12)));
% PD for every band start k_lo (with k_hi fixed): uses PD = sum sigma^2 |u|^2 / n
E = (abs(U).^2) .* (s.'.^2) / n;              % m x n, energy of each rank per pixel
conK = zeros(1, k_hi);
for klo = 1:k_hi, conK(klo) = contrast(sum(E(:, klo:k_hi), 2)); end
[cbest, klo_best] = max(conK);
pdRaw = pdOf(S); pdSVD = sum(E(:, klo_best:k_hi), 2);
pdFD = pdOf(fd(S));
cI3 = zeros(numel(fcs), 3); pdI = cell(1, numel(fcs));
for j = 1:numel(fcs), pdI{j} = iirPD(S, fcs(j)); cI3(j,:) = c3(pdI{j}); end
[cIb, jI] = max(cI3(:,1)); pdIIR = pdI{jI};
fprintf('\nPower Doppler contrast re background (dB):   both   slow   fast vessel\n');
fprintf('  raw (no filter)                        %6.1f %6.1f %6.1f\n', c3(pdRaw));
fprintf('  SVD band %2d-%2d (oracle k_lo)           %6.1f %6.1f %6.1f\n', klo_best, k_hi, c3(pdSVD));
fprintf('  first difference                       %6.1f %6.1f %6.1f\n', c3(pdFD));
for j = 1:numel(fcs)
    fprintf('  IIR Butterworth HP, zero-phase %3d Hz  %6.1f %6.1f %6.1f\n', fcs(j), cI3(j,:));
end
for j = 1:numel(pdeg)
    fprintf('  polynomial regression, degree %d        %6.1f %6.1f %6.1f\n', pdeg(j), c3(regPD(S, pdeg(j))));
end
fprintf('Doppler resolution of the ensemble PRF/n = %.1f Hz (%.2f mm/s axial)\n', PRF/n, PRF/n*lam/2);

%% 6. Velocity sweep: the best band moves, the high-pass fails
fprintf('\nContrast both / slow / fast vessel (dB) versus tissue velocity\n');
fprintf(' vt(mm/s) 2vt/lam(Hz)  oracle band   SVD band             best IIR (fc)              best regr. (p)   1st diff\n');
conCurves = zeros(numel(vts), 40);
for iv = 1:numel(vts)
    Sv = S_all{iv}; [Uv, Sv2, ~] = svd(Sv, 'econ'); sv = diag(Sv2);
    kh = find(sv > 1.05*edgeMP, 1, 'last');
    Ev = (abs(Uv).^2) .* (sv.'.^2) / n;
    cc = -inf(1, kh); for klo = 1:kh, cc(klo) = contrast(sum(Ev(:, klo:kh), 2)); end
    [cb, kb] = max(cc); conCurves(iv, 1:min(40,kh)) = cc(1:min(40,kh));
    ci = -inf; for j = 1:numel(fcs), pj = iirPD(Sv, fcs(j)); cj = contrast(pj); if cj > ci, ci = cj; jb = j; pib = pj; end, end
    cr = -inf; for j = 1:numel(pdeg), cj = contrast(regPD(Sv, pdeg(j))); if cj > cr, cr = cj; jr = j; end, end
    fprintf('  %5.0f     %6.1f       %2d-%2d   %5.1f %5.1f %5.1f    %5.1f %5.1f %5.1f (%3d Hz)    %5.1f (%d)       %5.1f\n', vts(iv), 2*vts(iv)/lam, kb, kh, ...
        c3(sum(Ev(:, kb:kh), 2)), c3(pib), fcs(jb), cr, pdeg(jr), contrast(pdOf(fd(Sv))));
end

%% 7. The fast route: eigendecomposition of the n x n Gram matrix S'S
% sigma_k^2 = eig(S'S), V = its eigenvectors; the filter S*Vb*Vb' never needs U.
% Forming S'S squares the condition number: a sigma_k below sqrt(eps)*sigma_1 is lost.
S = S_all{find(vts == vt_main)};
G = S'*S; G = (G + G')/2;                     % n x n Hermitian, cost O(m n^2)
[W, L] = eig(G); [lg, o] = sort(real(diag(L)), 'descend'); W = W(:, o);
sgD = sqrt(max(lg, 0));
Gs = single(S)'*single(S); Gs = (Gs + Gs')/2;
sgS = sqrt(max(double(sort(real(eig(Gs)), 'descend')), 0));
reS = abs(sgS - s)./s; kbad = find(reS > 0.01, 1);
Vb = W(:, klo_best:k_hi); pdGram = sum(abs(S*Vb).^2, 2)/n;   % = (1/n) sum_t |S Vb Vb'|^2, since Vb'Vb = I
fprintf('\nGram route (eig of %dx%d S''S): max rel. error of sigma_1..k_hi = %.1e (double)\n', n, n, max(abs(sgD(1:k_hi) - s(1:k_hi))./s(1:k_hi)));
fprintf('  PD from Gram route vs svd: max rel. difference %.1e\n', max(abs(pdGram - pdSVD)./pdSVD));
fprintf('  single precision: first sigma with > 1 %% error is k = %d at %.1f dB re sigma_1 (sqrt(eps_single) = %.1f dB)\n', ...
    kbad, 20*log10(s(kbad)/s(1)), 20*log10(sqrt(eps('single'))));

%% 8. Figures
figure('Position', [40 40 1100 640]);
subplot(2,3,1);
plot(20*log10(s/s(1)), 'k.-'); hold on; grid on;
plot([1 n], 20*log10(edgeMP/s(1))*[1 1], 'b--');
plot([klo_best klo_best], [-90 0], 'r'); plot([k_hi k_hi], [-90 0], 'r');
xlabel('singular value index k'); ylabel('\sigma_k (dB re \sigma_1)'); ylim([-90 5]);
title(sprintf('Singular values, vt = %g mm/s', vt_main)); legend('\sigma_k', 'MP noise edge', 'band', 'Location', 'northeast');
subplot(2,3,2);
Pn = 10*log10(Pk./max(Pk, [], 1));
imagesc(1:n, fax, Pn, [-30 0]); axis xy; colormap(gray); hold on;
plot([1 n], 2*vt_main/lam*[1 1], 'r--'); plot([1 n], -2*vt_main/lam*[1 1], 'r--');
xlabel('rank k'); ylabel('Doppler frequency (Hz)'); title('Spectra of temporal singular vectors (dB)');
subplot(2,3,3);
for iv = 1:numel(vts)
    kk = find(isfinite(conCurves(iv,:)) & conCurves(iv,:) ~= 0);
    plot(kk, conCurves(iv,kk), 'LineWidth', 1.5); hold on;
end
grid on; xlabel('k_{lo} (k_{hi} from MP edge)'); ylabel('contrast (dB)'); xlim([1 30]);
title('Best k_{lo} moves with tissue velocity');
legend(arrayfun(@(v) sprintf('%g mm/s', v), vts, 'UniformOutput', false), 'Location', 'southeast');
imgs = {pdRaw, pdSVD, pdIIR};
ttl = {'raw PD', sprintf('SVD band %d-%d', klo_best, k_hi), sprintf('IIR high-pass %d Hz', fcs(jI))};
for j = 1:3
    subplot(2,3,3+j);
    P = reshape(imgs{j}, NX, NZ).';           % rows = depth
    imagesc(xg, zg, 10*log10(P/max(P(:))), [-35 0]); axis image;
    xlabel('lateral (mm)'); ylabel('depth (mm)');
    title(sprintf('%s: %.1f dB', ttl{j}, contrast(imgs{j})));
end

%% m06_noise_floor_mp.m  --  Module 06: the noise floor of the SVD
% Purpose : Singular values of a pure-noise matrix fill the Marchenko-Pastur
%           bulk  sigma*(sqrt(m)-sqrt(n)) ... sigma*(sqrt(m)+sqrt(n)).
%           A rank-1 coherent component of strength theta separates from the
%           bulk only if theta > sigma*(m*n)^(1/4) (BBP transition).
% Change  : m, n (Casorati shape: pixels x frames), sig (noise std),
%           ratios (spike strength / threshold).
% Expect  : measured noise edges within ~2 % of theory for m = 1600, n = 80;
%           sigma_1 and singular-vector overlaps match the BBP formulas.
% Runtime : ~3 s in Octave.
clear; close all; rng(1);

%% 1. Pure complex Gaussian noise: bulk edges
sig = 1.0;                                  % noise std per entry, E|z|^2 = sig^2
shapes = [1600 80; 1600 400; 4096 128; 96 96];
fprintf('--- pure noise, complex Gaussian, sigma = %.1f\n', sig);
fprintf('   m      n   max(s) theory  min(s) theory  spread(dB)\n');
for i = 1:size(shapes,1)
    m = shapes(i,1); n = shapes(i,2);
    X = sig*(randn(m,n) + 1i*randn(m,n))/sqrt(2);
    s = svd(X);
    fprintf('%5d %5d  %6.2f %6.2f  %6.2f %6.2f   %5.1f\n', m, n, s(1), sig*(sqrt(m)+sqrt(n)), ...
        s(end), sig*(sqrt(m)-sqrt(n)), 20*log10(s(1)/s(end)));
end
fprintf('Gavish-Donoho optimal hard threshold (square, known sigma): %.3f*sqrt(n)*sigma (edge is 2*sqrt(n)*sigma)\n', 4/sqrt(3));

%% 2. Histogram of singular values against the Marchenko-Pastur density
m = 1600; n = 400; c = n/m;
X = sig*(randn(m,n) + 1i*randn(m,n))/sqrt(2);
s = svd(X);
a = (1-sqrt(c))^2; b = (1+sqrt(c))^2;      % eigenvalue support of X'X/(sig^2 m)
sg = linspace(sig*sqrt(m*a), sig*sqrt(m*b), 400);
x = sg.^2/(sig^2*m);
p = sqrt(max((b-x).*(x-a),0))./(2*pi*c*x) .* 2.*sg/(sig^2*m);   % density of s
figure('Position',[100 100 900 340]);
subplot(1,2,1);
[cnt, ctr] = hist(s, 30); bw = ctr(2)-ctr(1);
bar(ctr, cnt/(sum(cnt)*bw), 1, 'FaceColor', [0.75 0.8 0.9]); hold on;
plot(sg, p, 'r', 'LineWidth', 2);
xlabel('singular value'); ylabel('density');
title(sprintf('Pure noise %dx%d: MP density', m, n));
legend('svd(X)', 'Marchenko-Pastur', 'Location', 'north');

%% 3. Rank-1 coherent component in noise: BBP detectability threshold
m = 1600; n = 80; c = n/m;
thc = sig*(m*n)^(1/4);                     % detection threshold on theta
edge = sig*(sqrt(m)+sqrt(n));
fprintf('\n--- rank-1 spike in %dx%d noise: theta_c = %.2f, bulk edge = %.2f\n', m, n, thc, edge);
fprintf('theta/theta_c  s1(meas) s1(pred)  |<v1,v>|^2 meas/pred  |<u1,u>|^2 meas/pred\n');
ratios = [0.5 0.8 1.0 1.2 1.5 2 3 4];
ntr = 10; s1m = zeros(size(ratios));
for i = 1:numel(ratios)
    th = ratios(i)*thc; acc = zeros(1,3);
    for r = 1:ntr
        u = randn(m,1) + 1i*randn(m,1); u = u/norm(u);   % spatial pattern
        v = randn(n,1) + 1i*randn(n,1); v = v/norm(v);   % time course
        Y = th*(u*v') + sig*(randn(m,n) + 1i*randn(m,n))/sqrt(2);
        [U,S,V] = svd(Y, 'econ');
        acc = acc + [S(1,1), abs(V(:,1)'*v)^2, abs(U(:,1)'*u)^2]/ntr;
    end
    tp = th/(sig*sqrt(m));                 % normalized strength
    if ratios(i) > 1
        s1p = sqrt((th^2 + sig^2*m)*(th^2 + sig^2*n))/th;
        pv = 1 - c*(1+tp^2)/(tp^2*(tp^2+c)); pu = 1 - (c+tp^2)/(tp^2*(tp^2+1));
    else
        s1p = edge; pv = 0; pu = 0;
    end
    s1m(i) = acc(1);
    fprintf('   %4.1f        %6.2f   %6.2f      %5.2f / %5.2f          %5.2f / %5.2f\n', ...
        ratios(i), acc(1), s1p, acc(2), pv, acc(3), pu);
end
subplot(1,2,2);
tt = linspace(0.05, 4.2, 200)*thc;
s1th = edge*ones(size(tt)); k = tt > thc;
s1th(k) = sqrt((tt(k).^2 + sig^2*m).*(tt(k).^2 + sig^2*n))./tt(k);
plot(tt/thc, s1th, 'r', 'LineWidth', 2); hold on;
plot(ratios, s1m, 'ko', 'MarkerFaceColor', 'k');
plot([0 4.2], [edge edge], 'b--');
xlabel('\theta / \theta_c'); ylabel('\sigma_1');
title('Spike emerges only above \theta_c = \sigma (mn)^{1/4}');
legend('BBP prediction', 'measured', 'bulk edge', 'Location', 'northwest');
fprintf('Note: theta_c = %.1f sigma is far below the bulk edge %.1f sigma: a coherent\n', thc/sig, edge/sig);
fprintf('component is detectable even when sigma_1 barely clears the noise bulk.\n');

%% m06_svd_geometry.m  --  Module 06: SVD geometry, eig vs svd, power iteration, DORT
% Purpose : (1) the unit circle maps to an ellipse with semi-axes sigma_i*u_i;
%           (2) eigenvalues and singular values agree only for normal matrices;
%           (3) power iteration converges as |lambda2/lambda1|^k;
%           (4) DORT: the SVD of an array transfer matrix gives one singular
%               value per point scatterer and a focusing law for each.
% Change  : A (2x2), the gap ratio r, scatterer positions/reflectivities.
% Expect  : sigma = 1.618, 0.618 for the shear; power-iteration slope = log10(r);
%           two significant singular values for two point scatterers.
% Runtime : ~2 s in Octave.
clear; close all; rng(1);

%% 1. Circle -> ellipse
A = [1.2 0.9; -0.3 0.6];
[U,S,V] = svd(A);
th = linspace(0, 2*pi, 361); C = [cos(th); sin(th)]; E = A*C;
fprintf('A = [%.2f %.2f; %.2f %.2f]\n', A');
fprintf('sigma = %.4f, %.4f   |det A| = %.4f = sigma1*sigma2 = %.4f\n', S(1,1), S(2,2), abs(det(A)), S(1,1)*S(2,2));
fprintf('||A V - U S|| = %.1e,  ||A||_F^2 = %.4f = sum sigma^2 = %.4f\n', norm(A*V - U*S), norm(A,'fro')^2, sum(diag(S).^2));
fprintf('v1 at %.1f deg, u1 at %.1f deg (signs of (u_i, v_i) pairs are arbitrary)\n', ...
    atan2(V(2,1),V(1,1))*180/pi, atan2(U(2,1),U(1,1))*180/pi);
figure('Position',[100 100 900 420]);
subplot(1,2,1);
plot(C(1,:), C(2,:), 'Color', [0.6 0.6 0.6]); hold on; plot(E(1,:), E(2,:), 'b', 'LineWidth', 2);
for i = 1:2
    plot([0 V(1,i)], [0 V(2,i)], '--', 'Color', [0.6 0.6 0.6], 'LineWidth', 1.5);
    plot([0 S(i,i)*U(1,i)], [0 S(i,i)*U(2,i)], 'r', 'LineWidth', 2);
end
axis equal; axis([-1.8 1.8 -1.8 1.8]); grid on; xlabel('x_1'); ylabel('x_2');
title('unit circle (grey) -> A x (blue); v_i (dashed) -> \sigma_i u_i (red)');

%% 2. Eigenvalues versus singular values
names = {'symmetric', 'non-symmetric', 'shear', 'rotation 35 deg x 1.1'};
t = 35*pi/180;
mats = {[1.5 0.5; 0.5 0.8], [1.2 0.9; 0.15 0.5], [1 1; 0 1], 1.1*[cos(t) -sin(t); sin(t) cos(t)]};
fprintf('\n%-24s %-26s %-16s %s\n', 'matrix', 'eig', 'svd', 'angle(v1,v2)');
for i = 1:numel(mats)
    M = mats{i}; [X,L] = eig(M); l = diag(L); s = svd(M);
    if isreal(l)
        ang = acos(min(1, abs(X(:,1)'*X(:,2))/(norm(X(:,1))*norm(X(:,2)))))*180/pi;
        fprintf('%-24s %8.4f %8.4f          %6.4f %6.4f   %5.1f deg\n', names{i}, l(1), l(2), s(1), s(2), ang);
    else
        fprintf('%-24s %6.3f%+6.3fj (pair)       %6.4f %6.4f   complex\n', names{i}, real(l(1)), abs(imag(l(1))), s(1), s(2));
    end
end

%% 3. Power iteration on an 8x8 symmetric matrix
n = 8; [Q,~] = qr(randn(n));
subplot(1,2,2);
rs = [0.5 0.7 0.95]; cols = {'b-', 'r--', 'k-.'}; hp = zeros(1, numel(rs));
for j = 1:numel(rs)
    r = rs(j);
    lam = [1, r, r*0.75.^(1:n-2)];
    A8 = Q*diag(lam)*Q';
    x = randn(n,1); K = 80; err = zeros(1,K+1);
    for k = 0:K
        x = x/norm(x); c1 = Q(:,1)'*x;
        err(k+1) = norm(x - c1*Q(:,1))/abs(c1);          % tan of angle to q1
        x = A8*x;
    end
    hp(j) = semilogy(0:K, err, cols{j}, 'LineWidth', 1.6); hold on;
    kk = 5:min(K, find(err > 1e-13, 1, 'last'));
    p = polyfit(kk, log10(err(kk+1)), 1);
    fprintf('power iteration r = %.2f: slope %.4f (theory %.4f), iterations to 1e-15 = %d\n', ...
        r, p(1), log10(r), ceil(-15/log10(r)));
end
ylim([1e-16 1e2]); grid on; xlabel('iteration k'); ylabel('tan \angle(x_k, q_1)');
title('Power iteration: error ~ r^k'); legend(hp, 'r = 0.5', 'r = 0.7', 'r = 0.95', 'Location', 'southwest');

%% 4. DORT: SVD of the inter-element transfer matrix
c0 = 1540; f0 = 5e6; lam0 = c0/f0; kw = 2*pi/lam0;   % m, rad/m
Ne = 64; pitch = 0.3e-3; xe = ((0:Ne-1) - (Ne-1)/2)*pitch;  % element x (m)
sc = [-2e-3 30e-3 1.0; 3e-3 32e-3 0.5];               % [x z reflectivity] of two scatterers
G = zeros(Ne, size(sc,1));
for q = 1:size(sc,1)
    d = sqrt((xe - sc(q,1)).^2 + sc(q,2)^2);
    G(:,q) = exp(-1i*kw*d(:))./sqrt(d(:));             % 2-D Green's function (far-field form)
end
K = G*diag(sc(:,3))*G.';                               % K = K.' by reciprocity
K = K + 1e-3*norm(K)*(randn(Ne) + 1i*randn(Ne))/sqrt(2*Ne);
[Uk,Sk,Vk] = svd(K); sk = diag(Sk);
fprintf('\nDORT, 64 elements, two scatterers (refl. 1.0, 0.5): sigma_1..4 / sigma_1 = %.3f %.3f %.4f %.4f\n', sk(1:4)/sk(1));
fprintf('lateral resolution lambda*z/D = %.2f mm; lateral separation of scatterers = %.1f mm\n', lam0*31e-3/(Ne*pitch)*1e3, (sc(2,1)-sc(1,1))*1e3);
figure('Position',[100 100 900 340]);
subplot(1,2,1); semilogy(sk/sk(1), 'ko-'); grid on; xlabel('index i'); ylabel('\sigma_i / \sigma_1');
title('DORT: one singular value per scatterer');
xg = linspace(-8e-3, 8e-3, 401); zf = 31e-3;          % back-propagate onto a line at 31 mm
subplot(1,2,2); hq = zeros(1,2);
for q = 1:2
    w = conj(Uk(:,q));                                  % phase-conjugate (time-reversed) law
    P = zeros(size(xg));
    for ix = 1:numel(xg)
        d = sqrt((xe - xg(ix)).^2 + zf^2);
        P(ix) = abs(sum(w(:).*exp(-1i*kw*d(:))./sqrt(d(:))));
    end
    hq(q) = plot(xg*1e3, 20*log10(P/max(P)), 'LineWidth', 1.6); hold on;
    [~, im] = max(P); fprintf('singular vector %d focuses at x = %.2f mm\n', q, xg(im)*1e3);
end
hq(3) = plot(sc(:,1)*1e3, [0 0], 'kv', 'MarkerFaceColor', 'k');
ylim([-40 3]); grid on; xlabel('lateral position (mm)'); ylabel('field (dB)');
title('Back-propagated u_1, u_2 at z = 31 mm'); legend(hq, 'conj(u_1)', 'conj(u_2)', 'scatterers', 'Location', 'northwest');

%% m06_tsvd_picard.m  --  Module 06: truncated SVD, Tikhonov and the Picard plot
% Purpose : Deconvolve a blurred, noisy 1-D reflectivity profile b = A x + noise,
%           where A is a 64 x 64 Gaussian blur (the same problem as Demo 06.5).
%           Compare the plain pseudo-inverse, truncated SVD (TSVD) and Tikhonov,
%           read the Picard plot (sigma_i, |u_i' b|, |u_i' b|/sigma_i) and see
%           why the truncation index sits where |u_i' b| meets the noise floor.
% Change  : blur width sg (samples), noise levels nl (dB re signal RMS), n.
% Expect  : kappa(A) = 5.7e12 (sigma_64/sigma_1 = 1.75e-13); MATLAB's default
%           pinv tolerance (7e-15) keeps all 64 terms, so pinv(A)*b has relative
%           error ~1e10 at -40 dB noise. Best TSVD k = 27 (error 0.40), and the
%           Picard rule "last i with |u_i' b0| above the noise std" predicts 26;
%           Tikhonov at its best alpha gives 0.40 too. At -60 dB the best k moves
%           only to 31: sigma_i of a Gaussian blur falls super-exponentially.
% Runtime : < 1 s.
clear; close all; rng(1);

%% 1. Forward model
n = 64; sg = 2.5;                              % samples, Gaussian blur std (samples)
[I, J] = ndgrid(1:n, 1:n);
A = exp(-(I - J).^2/(2*sg^2))/(sqrt(2*pi)*sg);  % symmetric Toeplitz blur, rows sum to ~1
x = zeros(n,1);
x(11:18) = 0.6;                                % a uniform layer
x(39:50) = 0.35 + 0.02*((39:50) - 39);         % a ramped layer
x(28) = 1; x(31) = -0.7;                       % two point reflectors 3 samples apart
b0 = A*x;

[U, S, V] = svd(A); s = diag(S);
fprintf('kappa(A) = sigma_1/sigma_n = %.2e  (sigma_n/sigma_1 = %.2e)\n', s(1)/s(end), s(end)/s(1));
fprintf('default rank/pinv tolerance max(m,n)*eps*sigma_1 = %.1e -> numerical rank %d of %d\n', ...
    max(size(A))*eps(s(1)), sum(s > max(size(A))*eps(s(1))), n);

%% 2. Spectral filters at two noise levels
nl = [-40 -60];                                % noise re signal RMS (dB)
rmsb = sqrt(mean(b0.^2));
fprintf('\n noise   pinv err   best TSVD k  err     best Tikhonov alpha  err\n');
for q = 1:numel(nl)
    sn = rmsb*10^(nl(q)/20);                   % noise std per sample
    b = b0 + sn*randn(n,1);
    beta = U'*b;                               % data coefficients u_i' b
    errK = zeros(n,1);
    for k = 1:n
        xk = V(:,1:k)*(beta(1:k)./s(1:k));
        errK(k) = norm(xk - x)/norm(x);
    end
    [eK, kB] = min(errK);
    alphas = logspace(-6, 0, 121); errA = zeros(size(alphas));
    for j = 1:numel(alphas)
        f = s.^2./(s.^2 + alphas(j)^2);        % Tikhonov filter factors
        errA(j) = norm(V*(f.*beta./s) - x)/norm(x);
    end
    [eA, jA] = min(errA);
    ep = norm(pinv(A)*b - x)/norm(x);
    fprintf(' %3d dB   %8.1e      %2d     %.3f        %.1e          %.3f\n', nl(q), ep, kB, eK, alphas(jA), eA);
    % Picard prediction: keep terms while |u_i' b0| > noise std
    kP = find(abs(U'*b0) > sn, 1, 'last');
    fprintf('          last i with |u_i''b0| > noise std: %d\n', kP);
    if q == 1, b1 = b; beta1 = beta; errK1 = errK; kB1 = kB; sn1 = sn; end
end

%% 3. Figures (first noise level)
figure('Position', [60 60 1000 380]);
subplot(1,2,1);
semilogy(1:n, s, 'k.-', 1:n, abs(beta1), 'mo', 1:n, abs(beta1)./s, 'b.'); hold on;
semilogy([1 n], sn1*[1 1], 'm--'); semilogy([kB1 kB1]+0.5, [1e-16 1e8], 'r-');
ylim([1e-16 1e8]); grid on; xlabel('index i'); ylabel('magnitude');
legend('\sigma_i', '|u_i^T b|', '|u_i^T b| / \sigma_i', 'noise std', 'best k', 'Location', 'southwest');
title(sprintf('Picard plot, noise %d dB', nl(1)));
subplot(1,2,2);
xk = V(:,1:kB1)*(beta1(1:kB1)./s(1:kB1));
plot(1:n, x, 'k--', 1:n, b1, 'm', 1:n, xk, 'b', 'LineWidth', 1.4); grid on;
xlabel('depth sample'); ylabel('reflectivity');
legend('true x', 'data b', sprintf('TSVD, k = %d', kB1), 'Location', 'northeast');
title(sprintf('TSVD estimate, relative error %.2f', errK1(kB1)));