{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Accelerated factorized backprojection (AFBP)\n",
    "\n",
    "FFBP and CFBP (see their examples) reduce the cost of backprojection by imaging\n",
    "subapertures on coarse grids and merging them recursively, paying for the speedup with\n",
    "interpolation error at every merge level. AFBP removes the interpolation completely by\n",
    "merging in the wavenumber domain:\n",
    "\n",
    "1. Split the aperture into subapertures and backproject each one onto the same\n",
    "   polar grid, decimated in azimuth. Unlike FFBP, the grid origin is the same same on all subapertures,\n",
    "   a target sits at the same $(r, \\theta)$ cell in every subaperture image.\n",
    "3. Because the grids share an origin, each subaperture image holds a contiguous\n",
    "   patch of the full image's azimuth wavenumber spectrum, aliased onto the\n",
    "   decimated grid. The full-resolution image is assembled by FFT of the subaperture\n",
    "   images, placing each patch at its true position in the full spectrum, and inverse\n",
    "   FFT. No interpolation needed.\n",
    "\n",
    "The fusion is a single level (no recursion), so the backprojection cost drops by the\n",
    "number of subapertures while the fusion costs a few FFTs. The output matches direct\n",
    "backprojection including pixel phase far more accurately than the recursive\n",
    "interpolating merges, which makes AFBP attractive on short apertures, either\n",
    "standalone or as the base layer of FFBP (`ffbp(..., afbp_nsub=...)`).\n",
    "\n",
    "The library exposes this as a single call, `torchbp.ops.afbp`. In this notebook the\n",
    "algorithm is built by hand from `backprojection_polar_2d` and a few lines of torch to\n",
    "demonstrate how it works.\n",
    "\n",
    "Reference: L. Zhang, H.-l. Li, Z.-j. Qiao and Z.-w. Xu, \"A Fast BP Algorithm With\n",
    "Wavenumber Spectrum Fusion for High-Resolution Spotlight SAR Imaging,\" IEEE Geoscience\n",
    "and Remote Sensing Letters, vol. 11, no. 9, pp. 1460-1464, Sept. 2014."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import hamming\n",
    "\n",
    "import torchbp\n",
    "from torchbp.ops import backprojection_polar_2d\n",
    "\n",
    "device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n",
    "print(\"Device:\", device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Scene and radar parameters\n",
    "\n",
    "A single point target imaged on a heavily oversampled polar grid zoomed around the\n",
    "target, same setup as the FFBP example. The platform track is along y (cross-range),\n",
    "boresight along x, `theta` is the sine of the azimuth angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 6e9          # RF center frequency (Hz)\n",
    "bw = 200e6        # RF bandwidth (Hz), positive for rising sweep\n",
    "tsweep = 100e-6   # Sweep length (s)\n",
    "fs = 2e6          # Sampling frequency (Hz)\n",
    "nsweeps = 512     # Number of pulses in the aperture\n",
    "oversample = 2    # Range FFT oversampling (reduces interpolation error)\n",
    "nsamples = int(fs * tsweep)\n",
    "\n",
    "# Heavily oversampled polar grid, zoomed on the target. theta is sin(angle).\n",
    "nr, ntheta = 256, 512\n",
    "grid = {\"r\": (96.0, 104.0), \"theta\": (-0.06, 0.06), \"nr\": nr, \"ntheta\": ntheta}\n",
    "\n",
    "# Single point target at 100 m, broadside.\n",
    "target_pos = torch.tensor([[100.0, 0.0, 0.0]], dtype=torch.float32, device=device)\n",
    "target_rcs = torch.ones((1, 1), dtype=torch.float32, device=device)\n",
    "\n",
    "# Straight platform track along y at 50 m height. Element spacing = lambda/4.\n",
    "pos = torch.zeros((nsweeps, 3), dtype=torch.float32, device=device)\n",
    "pos[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps, device=device) * 0.25 * 3e8 / fc\n",
    "pos[:, 2] = 50.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Simulate and range-compress the data\n",
    "\n",
    "Identical to the backprojection example. Generate raw data, apply a range window, and\n",
    "FFT to range-compress. `data_fmod` shifts the data spectrum to DC and `alias_fmod`\n",
    "centers the image spectrum, both reduce interpolation error inside backprojection and\n",
    "must be passed consistently everywhere. The wavenumber fusion accounts for both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw = torchbp.util.generate_fmcw_data(target_pos, target_rcs, pos, fc, bw, tsweep, fs)\n",
    "w = torch.tensor(hamming(raw.shape[-1])[None, :], dtype=torch.float32, device=device)\n",
    "data = torch.fft.ifft(raw * w, dim=-1, n=nsamples * oversample)\n",
    "\n",
    "data_fmod = -math.pi * (1 - (oversample - 1) / oversample)\n",
    "data = data * torch.exp(1j * data_fmod * torch.arange(data.shape[-1], device=device))[None, :]\n",
    "\n",
    "r_res = 3e8 / (2 * abs(bw) * oversample)           # range bin size in the data\n",
    "dr = (grid[\"r\"][1] - grid[\"r\"][0]) / grid[\"nr\"]    # range bin size in the image\n",
    "im_margin = oversample * r_res / dr - 1\n",
    "alias_fmod = -math.pi * (1 - im_margin / (1 + im_margin))\n",
    "print(f\"r_res={r_res:.3f} m, image dr={dr:.3f} m, data_fmod={data_fmod:.3f}, alias_fmod={alias_fmod:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "A small helper to display a polar image in dB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "extent = [*grid[\"r\"], *grid[\"theta\"]]\n",
    "\n",
    "def show_polar(img, title, ax=None, dyn=30, ref_max=None):\n",
    "    img = img.detach().cpu()\n",
    "    db = 20 * torch.log10(torch.abs(img) + 1e-12)\n",
    "    m = db.max() if ref_max is None else ref_max\n",
    "    if ax is None:\n",
    "        _, ax = plt.subplots()\n",
    "    im = ax.imshow(db.numpy().T, origin=\"lower\", vmin=m - dyn, vmax=m,\n",
    "                   extent=extent, aspect=\"auto\")\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"Range (m)\")\n",
    "    ax.set_ylabel(\"Angle (sin)\")\n",
    "    return m\n",
    "\n",
    "def relerr(a, b):\n",
    "    return (torch.linalg.norm(a - b) / torch.linalg.norm(a)).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Reference: Direct backprojection\n",
    "\n",
    "`dealias=True` removes the range carrier $\\exp(j \\pi k_\\text{eff} d(r))$, $d(r) =\n",
    "\\sqrt{r^2 + z_0^2}$, from the image. This is required for AFBP subapertures to get a compact spectrum (see SAR image interpolation example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_bp = backprojection_polar_2d(\n",
    "    data, grid, fc, r_res, pos, dealias=True, data_fmod=data_fmod, alias_fmod=alias_fmod\n",
    ").squeeze()\n",
    "\n",
    "m = show_polar(img_bp, \"Direct backprojection (reference)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Step 1: Backproject the subapertures on the shared coarse grid\n",
    "\n",
    "Split the aperture into `NSUB` parts and backproject each onto the global grid with\n",
    "`ntheta / NSUB` azimuth samples. This is the key difference to FFBP: the positions are\n",
    "not re-centered per subaperture and every subaperture uses the\n",
    "same grid extents, so all images are in the same coordinate frame. Each subaperture has `1 / NSUB` of the pulses\n",
    "and `1 / NSUB` of the pixels, and there are `NSUB` of them. The total backprojection\n",
    "cost is `1 / NSUB` of the direct one.\n",
    "\n",
    "The subaperture images look identical to FFBP's at first glance, coarse azimuth\n",
    "resolution, full range resolution, but the target sits at the same cell in each and,\n",
    "each image keeps the azimuth carrier of its off-center aperture position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "NSUB = 4\n",
    "n_c = ntheta // NSUB\n",
    "grid_c = dict(grid, ntheta=n_c)\n",
    "n = nsweeps // NSUB\n",
    "\n",
    "imgs, x_u = [], []\n",
    "for u in range(NSUB):\n",
    "    sl = slice(u * n, (u + 1) * n)\n",
    "    img_u = backprojection_polar_2d(\n",
    "        data[sl], grid_c, fc, r_res, pos[sl],\n",
    "        dealias=True, data_fmod=data_fmod, alias_fmod=alias_fmod,\n",
    "    ).squeeze()\n",
    "    imgs.append(img_u)\n",
    "    x_u.append(float(pos[sl, 1].mean()))  # along-track subaperture centers\n",
    "imgs = torch.stack(imgs)\n",
    "\n",
    "print(\"subaperture image shape:\", tuple(imgs.shape[1:]), \"(full grid is\", (nr, ntheta), \")\")\n",
    "print(\"subaperture centers x_u:\", [f\"{x:.2f}\" for x in x_u])\n",
    "fig, axs = plt.subplots(1, 2, figsize=(11, 4))\n",
    "show_polar(imgs[0], \"Subaperture 0\", ax=axs[0])\n",
    "show_polar(imgs[1], \"Subaperture 1\", ax=axs[1])\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Why the subapertures are spectrum patches\n",
    "\n",
    "For a pixel at $(r, \\alpha)$ ($\\alpha = \\sin \\theta$) and a pulse at along-track\n",
    "position $x$, the two-way distance is $d = \\sqrt{r^2 - 2 r \\alpha x \\cdot r/d_s + \\dots}$\n",
    "and to first order in $x$ the dealiased image of subaperture $u$ carries the azimuth\n",
    "phase\n",
    "\n",
    "$$\\exp\\!\\big({-j k \\, x_u \\, \\tfrac{r}{d_s(r)} \\, (\\alpha - \\alpha_p)}\\big),\n",
    "\\qquad d_s(r) = \\sqrt{r^2 + z_0^2},$$\n",
    "\n",
    "where $k = 4\\pi f_c / c - f_\\text{mod,data} / \\Delta_r$ is the effective range\n",
    "wavenumber of the image and $\\alpha_p$ the target angle. A linear phase in $\\alpha$ is\n",
    "a shifted spectrum: subaperture $u$ occupies the azimuth frequency band centered at\n",
    "$-k\\, x_u\\, r/d_s / (2\\pi)$ cycles per unit $\\alpha$, with width set by the subaperture\n",
    "length. The factor $r / d_s$ is the ground-to-slant range ratio: the classical\n",
    "algorithm assumes a slant-plane image where it is 1, on the torchbp ground-range grid\n",
    "it scales the patch positions (here evaluated at the swath center; the library op\n",
    "follows its range variation with range blocks).\n",
    "\n",
    "The plot shows the azimuth spectrum of each subaperture image formed on the full\n",
    "grid: four flat patches that tile the full aperture band (dashed lines are\n",
    "the predicted patch edges). The coarse `ntheta / NSUB` grid of step 1 samples each\n",
    "patch critically, that is why one subaperture needs only `1 / NSUB` of the azimuth\n",
    "samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "C0 = 299792458.0\n",
    "krc = 4 * math.pi * fc / C0 - data_fmod / r_res  # effective range wavenumber (rad/m)\n",
    "z0 = float(pos[:, 2].mean())\n",
    "r_mid = 0.5 * (grid[\"r\"][0] + grid[\"r\"][1])\n",
    "fac = r_mid / math.sqrt(r_mid**2 + z0**2)        # ground-to-slant factor\n",
    "dalpha = (grid[\"theta\"][1] - grid[\"theta\"][0]) / ntheta\n",
    "L_sub = float(pos[n, 1] - pos[0, 1])             # subaperture length\n",
    "\n",
    "nu = np.fft.fftshift(np.fft.fftfreq(ntheta, d=dalpha))\n",
    "plt.figure(figsize=(9, 4))\n",
    "for u in range(NSUB):\n",
    "    sl = slice(u * n, (u + 1) * n)\n",
    "    img_f = backprojection_polar_2d(\n",
    "        data[sl], grid, fc, r_res, pos[sl],\n",
    "        dealias=True, data_fmod=data_fmod, alias_fmod=alias_fmod,\n",
    "    ).squeeze()\n",
    "    s = np.abs(np.fft.fftshift(np.fft.fft(img_f.cpu().numpy(), axis=-1), axes=-1)).mean(0)\n",
    "    line, = plt.plot(nu, 20 * np.log10(s + 1e-9), label=f\"subaperture {u}\")\n",
    "    center = -krc * x_u[u] * fac / (2 * math.pi)\n",
    "    hw = krc * L_sub * fac / (4 * math.pi)\n",
    "    plt.axvline(center - hw, color=line.get_color(), ls=\"--\", lw=0.8)\n",
    "    plt.axvline(center + hw, color=line.get_color(), ls=\"--\", lw=0.8)\n",
    "s = np.abs(np.fft.fftshift(np.fft.fft(img_bp.cpu().numpy(), axis=-1), axes=-1)).mean(0)\n",
    "plt.plot(nu, 20 * np.log10(s + 1e-9), \"k\", lw=1, label=\"direct BP\")\n",
    "plt.xlim(-300, 300); plt.ylim(-50, 0)\n",
    "plt.xlabel(\"Azimuth frequency (cycles / unit sin theta)\")\n",
    "plt.ylabel(\"dB\")\n",
    "plt.title(\"Subaperture azimuth spectra tile the full aperture band\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Step 2: Wavenumber fusion\n",
    "\n",
    "On the coarse grid each patch is aliased, but nothing is lost: the coarse grid samples\n",
    "every `NSUB`-th point of the fine grid, so a coarse spectrum bin is the sum of the fine\n",
    "spectrum aliases, $S_c[k] = \\frac{1}{N_\\text{sub}} \\sum_j S_f[k + j\\, n_c]$. Inside the\n",
    "patch the other aliases are just the patch's own leakage tails, so the fine bin $m$ is\n",
    "recovered with plain indexing:\n",
    "\n",
    "$$S_f[m] \\approx N_\\text{sub} \\; S_c[m \\bmod n_c].$$\n",
    "\n",
    "There is no explicit circular shift and no interpolation, the patch placement is\n",
    "the index arithmetic. Two refinements make it accurate:\n",
    "\n",
    "- Trapezoidal placement: the patch center is proportional to the range wavenumber\n",
    "  $k_r$, not just its center value, so the assembly runs per range-wavenumber row\n",
    "  (range FFT first, `x_eq` below is the along-track position owning each row/bin).\n",
    "- Overlapping regions summed: a discrete patch has leakage tails past its rect\n",
    "  edges. Assigning each fine bin to exactly one subaperture (hard tiling) truncates\n",
    "  the tails at every seam, which costs surprisingly much (~-20 dB). Instead every\n",
    "  subaperture reconstructs a full coarse band centered on its own patch and the\n",
    "  overlapping regions of adjacent subapertures are summed, near a seam the true\n",
    "  spectrum really is the sum of both patches edge transitions.\n",
    "\n",
    "Note that linear track and $(r, \\sin \\theta)$ grid assumptions are used here. If either\n",
    "assumption doesn't hold then the subaperture spectrum relationship is more complicated.\n",
    "\n",
    "Everything is exact bookkeeping of the FFT bin frequencies: `nu_c` is where the\n",
    "`data_fmod`/`alias_fmod` modulations put the range spectrum center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = ntheta\n",
    "Sa = torch.fft.fft(imgs, dim=-1)   # azimuth FFT -> aliased patches\n",
    "S = torch.fft.fft(Sa, dim=-2)      # range FFT -> range wavenumber rows\n",
    "\n",
    "# Range wavenumber of every row, unwrapped around the modulation-shifted center.\n",
    "nu_c = alias_fmod / (2 * dr) - data_fmod / (2 * r_res)\n",
    "fr = torch.fft.fftfreq(nr, d=dr, dtype=torch.float64, device=device)\n",
    "half = 1 / (2 * dr)\n",
    "kr = krc + 2 * math.pi * (torch.remainder(fr - nu_c + half, 2 * half) - half)\n",
    "\n",
    "# Along-track position owning each (range wavenumber row, fine azimuth bin).\n",
    "nua = torch.fft.fftfreq(N, d=dalpha, dtype=torch.float64, device=device)\n",
    "x_eq = -2 * math.pi * nua[None, :] / (kr[:, None] * fac)\n",
    "# Half width of one coarse band in the same units.\n",
    "x_half = math.pi * n_c / (N * dalpha) / (kr * fac)\n",
    "\n",
    "cols = torch.arange(N, device=device) % n_c\n",
    "fine = S.new_zeros((nr, N))\n",
    "for u in range(NSUB):\n",
    "    mask = (x_eq - x_u[u]).abs() <= x_half[:, None].abs()\n",
    "    fine += S[u, :, cols] * mask\n",
    "fine *= NSUB\n",
    "\n",
    "img_afbp = torch.fft.ifft(torch.fft.ifft(fine, dim=-2), dim=-1)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(11, 4))\n",
    "show_polar(img_bp, \"Direct backprojection\", ax=axs[0], ref_max=m)\n",
    "show_polar(img_afbp, f\"AFBP, {NSUB} subapertures (manual)\", ax=axs[1], ref_max=m)\n",
    "plt.tight_layout();\n",
    "\n",
    "print(f\"relative error vs direct BP: {relerr(img_bp, img_afbp):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "The fused image matches the reference. The azimuth cut through the peak overlays the\n",
    "direct result and the error floor sits far below the FFBP interpolation error on the\n",
    "same scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "pk = np.unravel_index(torch.abs(img_bp).cpu().argmax(), img_bp.shape)\n",
    "theta_axis = np.linspace(grid[\"theta\"][0], grid[\"theta\"][1], ntheta)\n",
    "cut_bp = 20 * torch.log10(torch.abs(img_bp[pk[0]]).cpu() + 1e-12)\n",
    "cut_af = 20 * torch.log10(torch.abs(img_afbp[pk[0]]).cpu() + 1e-12)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
    "axs[0].plot(theta_axis, cut_bp - cut_bp.max(), label=\"direct BP\")\n",
    "axs[0].plot(theta_axis, cut_af - cut_bp.max(), \"--\", label=\"AFBP (manual)\")\n",
    "axs[0].set_xlim(-0.02, 0.02); axs[0].set_ylim(-40, 1)\n",
    "axs[0].set_xlabel(\"Angle (sin)\"); axs[0].set_ylabel(\"dB\")\n",
    "axs[0].set_title(\"Azimuth cut through peak\"); axs[0].legend()\n",
    "show_polar(img_bp - img_afbp, \"Difference (AFBP - direct)\", ax=axs[1], ref_max=m, dyn=80)\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Library op `torchbp.ops.afbp`\n",
    "\n",
    "Everything above is what `torchbp.ops.afbp` does internally, plus:\n",
    "\n",
    "- a theta guard band (`guard_theta`): the subaperture images treat theta as\n",
    "  circular, so a target near the grid edge wraps around with a broadband spectral\n",
    "  skirt. The guard moves the wrap point outside the scene and is cropped from the\n",
    "  output\n",
    "- range blocks that follow the ground-to-slant factor $r / d_s(r)$ over a wide\n",
    "  swath instead of using one swath-center value\n",
    "- uneven pulse/theta splits, batched subaperture backprojection, and validity\n",
    "  warnings (patch fit on the decimated grid, quadratic phase error bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_lib4 = torchbp.ops.afbp(data, grid, fc, r_res, pos, nsub=4, dealias=True,\n",
    "                            data_fmod=data_fmod, alias_fmod=alias_fmod)\n",
    "img_lib8 = torchbp.ops.afbp(data, grid, fc, r_res, pos, nsub=8, dealias=True,\n",
    "                            data_fmod=data_fmod, alias_fmod=alias_fmod)\n",
    "\n",
    "print(f\"afbp nsub=4 vs direct BP: {relerr(img_bp, img_lib4):.4f}\")\n",
    "print(f\"afbp nsub=8 vs direct BP: {relerr(img_bp, img_lib8):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## AFBP as the FFBP base layer\n",
    "\n",
    "The single-level fusion assumes a straight track and a quadratic-phase bound\n",
    "$\\alpha_\\text{max} \\, l \\le r_0 / 4$ per subaperture, so it cannot replace FFBP for\n",
    "very long, curved apertures. There it can still speed up the base level, where the\n",
    "subapertures are short and the direct backprojections dominate the cost:\n",
    "`ffbp(..., afbp_nsub=N)` computes each base-level image with AFBP, cutting the base\n",
    "cost roughly `N` times.\n",
    "\n",
    "Note what this does and does not buy. As a drop-in addition it does not improve\n",
    "accuracy: the FFBP merge interpolation error dominates, and the AFBP base only adds a\n",
    "much smaller error on top (compare the first two lines below). The benefit is the\n",
    "cost trade, because the base level gets `N` times cheaper, an interpolating merge\n",
    "stage can be dropped at roughly constant total cost. `stages=1, afbp_nsub=4` has about\n",
    "the same base-level cost as `stages=3` but only one interpolating merge instead of\n",
    "three, so it runs at the deep-recursion speed with the shallow-recursion accuracy\n",
    "(last line)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "for stages, ans in [(2, 1), (2, 4), (1, 1), (1, 4)]:\n",
    "    img = torchbp.ops.ffbp(data, grid, fc, r_res, pos, stages=stages, dealias=True,\n",
    "                           data_fmod=data_fmod, alias_fmod=alias_fmod, afbp_nsub=ans)\n",
    "    print(f\"ffbp stages={stages} afbp_nsub={ans} vs direct BP: {relerr(img_bp, img):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## Benefit and downside\n",
    "\n",
    "AFBP gives the subaperture speedup of one factorization level with no interpolation\n",
    "error: the fusion is exact up to the discrete patch-edge effects handled by the\n",
    "overlap-sum and the guard band, and in practice sits one to two orders of magnitude\n",
    "below the FFBP merge error while being faster on short apertures. It also preserves\n",
    "pixel phase, so interferometric or autofocus processing downstream is unaffected, and\n",
    "the whole pipeline is differentiable with respect to the data.\n",
    "\n",
    "The downsides are the assumptions the spectrum placement is built on: an\n",
    "approximately straight track along y, patch positions linear in $\\alpha$ (the\n",
    "quadratic phase bound $\\alpha_\\text{max} l \\le r_0/4$), and an azimuth grid fine\n",
    "enough that one subaperture patch fits the decimated band\n",
    "($\\Delta\\alpha \\le \\lambda_\\text{min} / (2 L / N_\\text{sub})$, the op warns when\n",
    "violated). Wide total angles, long curved tracks, or strong altitude variation are handled by FFBP; there `ffbp(..., afbp_nsub=N)` trades interpolating merge stages for\n",
    "the interpolation-free fusion, keeping the deep-recursion speed with fewer merge stages worth of error."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}