{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Choosing fc for backprojection\n",
    "\n",
    "`fc` in backprojection is the carrier of the per-pixel phase compensation\n",
    "$e^{+j 4\\pi f_c R / c}$ applied to the range-compressed data. It must match the\n",
    "phase-reference convention of the data, and that convention is set by the waveform and how the\n",
    "range profile is read out:\n",
    "\n",
    "* Interpolating readout (= backprojection): reference = first ADC sample, so\n",
    "  `fc = flow` for a rising sweep, `fc = fhigh` for a falling sweep.\n",
    "* Interpolating readout + mid-sweep factor $e^{-j\\pi f T}$: reference = mid\n",
    "  sweep, so `fc = fcenter`, both sweep directions.\n",
    "* Rising sweep, ADC samples time-reversed: reference = last ADC sample, so\n",
    "  `fc = flow + K(N-1)/fs` ≈ `fhigh`.\n",
    "* Rising sweep, delayed ADC start: reference = first ADC sample, so `fc` = the\n",
    "  frequency at that first sample.\n",
    "* FFT + fixed range bin (classic range–Doppler analysis): reference = window\n",
    "  centre (via leakage), so `fc = fcenter`.\n",
    "* Matched-filter pulse compression, DC-centred baseband: reference =\n",
    "  demodulation frequency (LO), so `fc = fcenter`.\n",
    "\n",
    "`flow` is the lowest frequency (start frequency for rising sweep) and `fhigh` is the highest frequency. For falling sweep the start frequency is `fhigh`.\n",
    "\n",
    "\"Phase reference at the first ADC sample\" is a property of the readout, not of\n",
    "the signal: the FFT correlates the whole windowed sweep against\n",
    "$e^{-j2\\pi f t}$, whose zero-phase instant is the first sample, so the measured\n",
    "phase is a fit to the entire sweep, quoted at that instant. `fc = flow` is not a sampling of the hardware frequency\n",
    "at one instant: sweep non-linearity matters through its window-weighted effect\n",
    "over the whole sweep and its value at the reference instant specifically carries no special weight.\n",
    "Only a deviation with a nonzero window-weighted mean shifts the effective reference.\n",
    "\n",
    "This example demonstrates:\n",
    "\n",
    "1. RF-sample-level simulation: the dechirped IF signal is\n",
    "   $e^{\\,j 2\\pi (f_{low}\\tau + K\\tau t - K\\tau^2/2)}$ with $t = 0$ at the\n",
    "   first ADC sample. A plain FFT therefore leaves the data referenced to the\n",
    "   sweep start, and backprojection must use `fc = flow` for a rising sweep.\n",
    "2. `fc` is phase bookkeeping: even a hypothetical DC-start sweep\n",
    "   (`flow = 0`, so `fc = 0` and the phase compensation is identically 1) focuses\n",
    "   perfectly. Azimuth resolution is set by the physical spectrum of the signal,\n",
    "   not by the `fc` parameter.\n",
    "3. Mid-sweep re-referencing: multiplying the conjugated range spectrum by\n",
    "   $e^{-j\\pi f T}$ (over the signed beat-frequency axis $f$, sweep length $T$)\n",
    "   moves the reference to mid sweep, making `fc = fcenter` exact and rising/falling\n",
    "   sweeps symmetric, so identical processing handles both sweep directions.\n",
    "4. Time reversal: reversing the ADC samples moves the phase reference to the frequency\n",
    "   at the last sample in the original order. `fc` must be changed to correctly focus the image\n",
    "   for start-referenced signal.\n",
    "5. Why fixed-bin analyses mislead: reading the phase history from one fixed\n",
    "   FFT bin while the target migrates through range bins is `fcenter`-referenced\n",
    "   through window leakage. This is how classical range-Doppler processing works.\n",
    "   That is a property of that readout, not of backprojection, which interpolates at the beat frequency.\n",
    "6. Pulse compression radar: correlation range compression time origin moves the range origin, not the phase\n",
    "   reference. `fc` equals the demodulation frequency, `fcenter` for the\n",
    "   usual DC-centred baseband convention, regardless of window placement or sample order."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Dechirp signal model\n",
    "\n",
    "Rising-sweep TX phase with chirp rate $K = B/T$:\n",
    "\n",
    "$$\\varphi_{tx}(t) = 2\\pi\\!\\left(f_{low}\\,t + \\tfrac{1}{2}K t^2\\right)$$\n",
    "\n",
    "The mixer forms IF $=$ TX $\\cdot$ conj(RX) with RX $=$ TX$(t-\\tau)$, $\\tau = 2R/c$:\n",
    "\n",
    "$$\\varphi_{IF}(t) = \\varphi_{tx}(t) - \\varphi_{tx}(t-\\tau)\n",
    " = 2\\pi\\!\\left(f_{low}\\,\\tau + K\\tau\\,t - \\tfrac{1}{2}K\\tau^2\\right),\n",
    " \\qquad t = 0 \\text{ at the first ADC sample.}$$\n",
    "\n",
    "Re-parameterizing the same signal about mid sweep, $t' = t - T/2$, gives\n",
    "$2\\pi(f_c\\tau + K\\tau t' - \\tfrac12 K\\tau^2)$ with $f_c = f_{low} + B/2$, identical signal but different `fc`. The choice of where $t=0$ affects the required `fc` value and this depends on the processing method.\n",
    "\n",
    "Range compression with the Fourier transform basis $e^{-j2\\pi f t}$ ($t=0$ at the first\n",
    "sample), evaluated at the beat frequency $f_b = K\\tau$, which is what\n",
    "backprojection's interpolation does, leaves\n",
    "\n",
    "$$\\angle X(f_b) = 2\\pi f_{low}\\tau - \\pi K \\tau^2 .$$\n",
    "\n",
    "In the conjugated-spectrum convention (assuming rising sweep) the target phase is\n",
    "$-4\\pi f_{low} R/c + \\pi K\\tau^2$, so the backprojection compensation\n",
    "$e^{+j4\\pi f_c R/c}$ cancels it only if `fc` = $f_{low}$. The last term ($ -\\pi K \\tau^2$) is the\n",
    "residual video phase (RVP): small, common to every `fc` choice. It can be compensated separately if necessary.\n",
    "\n",
    "Mid-sweep factor: on the conjugated spectrum the target sits at $f = K\\tau$,\n",
    "so multiplying by $e^{-j\\pi f T}$ contributes $e^{-j\\pi K\\tau T} = e^{-j\\pi\n",
    "B\\tau}$, moving the reference phase to $-4\\pi(f_{low}+B/2)R/c$: exactly\n",
    "$f_c = f_{center}$. A falling sweep puts the target at $f = -K\\tau$ and the same\n",
    "factor yields $-4\\pi(f_{stop}-B/2)R/c$ — again $f_{center}$. One factor, both\n",
    "sweep directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import torch\n",
    "import torchbp\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import hamming, hanning\n",
    "\n",
    "c0 = 3e8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## RF-sample-level FMCW simulation\n",
    "\n",
    "TX chirp sampled at 8 GHz, delayed RX, calculate mixer product, zero-phase\n",
    "low-pass filter (no delay to avoid modifying phase), analytic signal, decimation to the 50 MHz\n",
    "ADC rate. No complex-exponential shortcuts used.\n",
    "\n",
    "A realistic sweep (1.0–1.5 GHz) is used here because a real mixer cannot dechirp a\n",
    "DC-start sweep: near $t \\approx \\tau$ the sum and difference products overlap at DC\n",
    "and no filter separates them. That is also why the DC-start sweep used later is\n",
    "hypothetical, there it is generated directly from the analytic IF model that\n",
    "this part validates.\n",
    "\n",
    "Lowpass filtering causes slight error at the edges because of edge-effects, but the match is good against the ideal IF signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow = 1.0e9 # sweep 1.0 -> 1.5 GHz\n",
    "bw1 = 0.5e9\n",
    "tsweep1 = 20e-6\n",
    "K1 = bw1 / tsweep1\n",
    "fc1 = flow + bw1 / 2\n",
    "\n",
    "fs_rf = 8e9 # RF simulation rate\n",
    "n_rf = int(round(tsweep1 * fs_rf))\n",
    "t_rf = np.arange(n_rf) / fs_rf\n",
    "dec = 160\n",
    "fs_adc = fs_rf / dec # 50 MHz ADC\n",
    "n_adc = n_rf // dec\n",
    "t_adc = np.arange(n_adc) / fs_adc\n",
    "phi0 = 0.1234 # Arbitrary TX phase offset\n",
    "\n",
    "def phi_tx(t):\n",
    "    return 2 * np.pi * (flow * t + 0.5 * K1 * t**2) + phi0\n",
    "\n",
    "def rf_dechirp(tau):\n",
    "    '''TX, delayed RX, mixer, LPF, Hilbert, decimate.'''\n",
    "    tx = np.cos(phi_tx(t_rf))\n",
    "    rx = np.where(t_rf >= tau, np.cos(phi_tx(t_rf - tau)), 0.0)\n",
    "    b, a = signal.butter(6, 20e6 / (fs_rf / 2))\n",
    "    beat = signal.filtfilt(b, a, tx * rx) # zero-phase LPF\n",
    "    return signal.hilbert(beat)[::dec]\n",
    "\n",
    "R_test = 45.0\n",
    "tau_test = 2 * R_test / c0\n",
    "sig = rf_dechirp(tau_test)\n",
    "\n",
    "# analytic IF model, t = 0 at the first ADC sample\n",
    "phi_model = 2 * np.pi * (flow * tau_test + K1 * tau_test * t_adc\n",
    "                         - 0.5 * K1 * tau_test**2)\n",
    "valid = (t_adc > tau_test + 200e-9) & (t_adc < tsweep1 - 200e-9)\n",
    "perr = np.angle(sig * np.exp(-1j * phi_model))\n",
    "print(f\"R = {R_test} m, tau = {tau_test*1e9:.1f} ns, \"\n",
    "      f\"beat = {K1*tau_test/1e6:.2f} MHz\")\n",
    "print(f\"IF phase error vs first-sample model: \"\n",
    "      f\"RMS {np.sqrt(np.mean(perr[valid]**2))*1e3:.2f} mrad, \"\n",
    "      f\"max {np.max(np.abs(perr[valid]))*1e3:.2f} mrad\")\n",
    "\n",
    "# sanity check: the mid-sweep parameterization is the SAME signal\n",
    "phi_mid = 2 * np.pi * (fc1 * tau_test + K1 * tau_test * (t_adc - tsweep1 / 2)\n",
    "                       - 0.5 * K1 * tau_test**2)\n",
    "\n",
    "print(f\"first-sample vs mid-sweep parameterization: max difference \"\n",
    "      f\"{np.max(np.abs(phi_model - phi_mid)):.2e} rad\")\n",
    "print(\"   (identical signal to float rounding)\")\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(t_adc[valid] * 1e6, perr[valid] * 180/np.pi)\n",
    "plt.xlabel(\"Fast time from first ADC sample (µs)\")\n",
    "plt.ylabel(\"Phase error (deg)\")\n",
    "plt.title(\"Phase error of RF-level dechirped IF vs model with t=0\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "### Which reference does the interpolating readout leave?\n",
    "\n",
    "Repeat the RF-level simulation at different target distances and read the (windowed,\n",
    "conjugated) DTFT out at the beat frequency, exactly what backprojection's\n",
    "interpolation does. The measured phase is compared against the candidate models\n",
    "$-2\\pi f_{ref}\\tau + \\pi K\\tau^2$. The candidates are far apart ($f_c$ vs\n",
    "$f_{low}$ differ by $2\\pi (B/2)\\tau \\approx$ hundreds of radians here), so the\n",
    "verdict is unambiguous: flow is the correct choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs = np.linspace(20, 90, 15)\n",
    "taus = 2 * Rs / c0\n",
    "win = hanning(n_adc)\n",
    "n = np.arange(n_adc)\n",
    "meas = np.empty(len(Rs))\n",
    "for i, tau in enumerate(taus):\n",
    "    s = rf_dechirp(tau)\n",
    "    fb = K1 * tau\n",
    "    X = (win * s * np.exp(-2j * np.pi * fb * n / fs_adc)).sum().conj()\n",
    "    meas[i] = np.angle(X)\n",
    "\n",
    "plt.figure(figsize=(8, 3.2))\n",
    "for name, f_ref in {\"flow\": flow, \"fc\": fc1, \"fhigh\": flow + bw1}.items():\n",
    "    model = -2 * np.pi * f_ref * taus + np.pi * K1 * taus**2\n",
    "    d = np.angle(np.exp(1j * (meas - model)))\n",
    "    print(f\"readout phase vs {name:6s} model: wrapped-difference \"\n",
    "          f\"RMS {np.sqrt(np.mean(d**2)):7.4f} rad\")\n",
    "    plt.plot(Rs, d, \"o-\", ms=4, label=name)\n",
    "plt.xlabel(\"Target range (m)\")\n",
    "plt.ylabel(\"Measured − model, wrapped (rad)\")\n",
    "plt.title(\"Tracked readout of RF-level data: phase reference is flow\")\n",
    "plt.legend(title=\"reference model\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## The readout functional decides the reference\n",
    "\n",
    "A hypothetical DC-start sweep (`flow = 0`, B = 1.5 GHz) makes the contrast\n",
    "maximal: an `flow`-referenced phase history is flat, while an\n",
    "`fcenter`-referenced one spans ~11 rad over this aperture. One point target\n",
    "migrates ~3.6 range bins across the aperture, and the per-pulse range profile is\n",
    "read out several ways:\n",
    "\n",
    "* fixed-bin: one integer FFT bin for the whole aperture. Once the target\n",
    "  walks off the bin, the window-leakage phase re-references the readout to the\n",
    "  window centre = mid sweep, so it behaves as `fcenter`. Its phase is only\n",
    "  meaningful near the main lobe, so the comparison is amplitude²-weighted. A\n",
    "  classic range–Doppler analysis on a fixed range gate measures this\n",
    "  functional, a common source of the incorrect conclusion that dechirped data\n",
    "  are inherently `fcenter`-referenced.\n",
    "* tracking: DTFT evaluated at the migrating beat frequency, i.e. what\n",
    "  backprojection's interpolation does. Referenced to the first-sample instant:\n",
    "  `flow` for a rising sweep, `fhigh` for a falling one.\n",
    "* tracking + mid-sweep re-referencing: the $e^{-j\\pi f T}$ factor on the\n",
    "  conjugated spectrum: `fcenter`, exactly, for both sweep directions.\n",
    "* tracking, time-reversed samples: flipping the ADC samples moves the\n",
    "  reference to the last-sample instant. The exact reference is the frequency at\n",
    "  the last sample, $f_{last} = f_{low} + K(N{-}1)/f_s$, one sample short of\n",
    "  $f_{stop}$. Visible in the table as the `f_last` column winning over `fhigh`\n",
    "  by a small margin. With the mid-sweep factor the reversed data land on\n",
    "  $f_c - K/f_s$ for the same one-sample reason.\n",
    "\n",
    "Models use the case-appropriate residual-video-phase sign (rising, conjugated:\n",
    "$+\\pi K\\tau^2$; falling, conjugated: $-\\pi K\\tau^2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow2, bw2, tsweep2 = 0.0, 1.5e9, 100e-6 # hypothetical DC-start sweep\n",
    "K2 = bw2 / tsweep2\n",
    "fc2, fhigh2 = bw2 / 2, flow2 + bw2\n",
    "\n",
    "v, R0, L_ap = 50.0, 50.0, 12.0\n",
    "nfast, nslow = 1024, 256\n",
    "fs2 = nfast / tsweep2\n",
    "t2 = np.arange(nfast) / fs2\n",
    "eta = np.linspace(-L_ap / (2 * v), L_ap / (2 * v), nslow)\n",
    "R_eta = np.sqrt((v * eta)**2 + R0**2)\n",
    "tau2 = 2 * R_eta / c0\n",
    "fb2 = K2 * tau2\n",
    "print(f\"target migrates {(R_eta.max()-R_eta.min())/(c0/(2*bw2)):.1f} \"\n",
    "      f\"range bins across the aperture\")\n",
    "\n",
    "w2 = hanning(nfast)\n",
    "def make_pulses(rising=True):\n",
    "    s = np.zeros((nslow, nfast), complex)\n",
    "    for i in range(nslow):\n",
    "        if rising:\n",
    "            ph = 2 * np.pi * (flow2 * tau2[i] + K2 * tau2[i] * t2\n",
    "                              - 0.5 * K2 * tau2[i]**2)\n",
    "        else:\n",
    "            ph = 2 * np.pi * (fhigh2 * tau2[i] - K2 * tau2[i] * t2\n",
    "                              + 0.5 * K2 * tau2[i]**2)\n",
    "        s[i] = w2 * np.exp(1j * ph)\n",
    "    return s\n",
    "\n",
    "s_rise = make_pulses(True)\n",
    "s_fall = make_pulses(False)\n",
    "s_rev = np.flip(s_rise, axis=1)               # time-reversed ADC samples\n",
    "\n",
    "def fixed_bin_conj(sig):\n",
    "    '''One integer FFT bin for the whole aperture (fixed range gate).'''\n",
    "    S = np.conj(np.fft.fft(sig, axis=1))\n",
    "    k = np.argmax(np.abs(S[nslow // 2]))\n",
    "    return np.unwrap(np.angle(S[:, k])), np.abs(S[:, k])\n",
    "\n",
    "def tracking_conj(sig, freqs, midsweep=False):\n",
    "    '''BP readout: conjugated DTFT AT the migrating beat frequency,\n",
    "    optionally with the mid-sweep factor exp(-1j*pi*f*T).'''\n",
    "    nn = np.arange(sig.shape[1])\n",
    "    out = np.empty(nslow, complex)\n",
    "    for i in range(nslow):\n",
    "        val = (sig[i] * np.exp(-2j * np.pi * freqs[i] * nn / fs2)).sum().conj()\n",
    "        if midsweep:\n",
    "            val *= np.exp(-1j * np.pi * freqs[i] * tsweep2)\n",
    "        out[i] = val\n",
    "    return np.unwrap(np.angle(out)), np.abs(out)\n",
    "\n",
    "f_last = flow2 + K2 * (nfast - 1) / fs2     # frequency at the LAST sample\n",
    "\n",
    "cases = {   # name: (readout, RVP sign of the conjugated model)\n",
    "    \"fixed-bin, rising\":             (fixed_bin_conj(s_rise), +1),\n",
    "    \"tracking, rising\":              (tracking_conj(s_rise, fb2), +1),\n",
    "    \"tracking, rising + midsweep\":   (tracking_conj(s_rise, fb2, True), +1),\n",
    "    \"tracking, falling\":             (tracking_conj(s_fall, -fb2), -1),\n",
    "    \"tracking, falling + midsweep\":  (tracking_conj(s_fall, -fb2, True), -1),\n",
    "    \"tracking, time-reversed\":       (tracking_conj(s_rev, -fb2), +1),\n",
    "    \"tracking, reversed + midsweep\": (tracking_conj(s_rev, -fb2, True), +1),\n",
    "}\n",
    "refs = {\"flow\": flow2, \"fc\": fc2, \"fhigh\": fhigh2,\n",
    "        \"f_last\": f_last, \"fc-K/fs\": fc2 - K2 / fs2}\n",
    "\n",
    "def rms_vs(f_ref, rvp_sign, meas_amp):\n",
    "    meas, amp = meas_amp\n",
    "    model = -2 * np.pi * f_ref * tau2 + rvp_sign * np.pi * K2 * tau2**2\n",
    "    d = meas - model\n",
    "    d = d - d[nslow // 2] # constants are unobservable\n",
    "    return np.sqrt(np.average(d**2, weights=amp**2))\n",
    "\n",
    "print(f\"\\nphase-history RMS mismatch (rad) vs reference model \"\n",
    "      f\"(fc span here = {2*np.pi*fc2*(tau2.max()-tau2.min()):.1f} rad):\")\n",
    "print(f\"{'case':32s}\" + \"\".join(f\"{r:>9s}\" for r in refs) + \"   winner\")\n",
    "for name, (meas_amp, rvp) in cases.items():\n",
    "    r = {rn: rms_vs(rf, rvp, meas_amp) for rn, rf in refs.items()}\n",
    "    row = \"\".join(f\"{r[rn]:9.3f}\" for rn in refs)\n",
    "    print(f\"{name:32s}{row}   {min(r, key=r.get)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 4.5))\n",
    "plot_cases = [(\"fixed-bin, rising\", \"-\"),\n",
    "              (\"tracking, rising\", \"-\"),\n",
    "              (\"tracking, rising + midsweep\", \"--\"),\n",
    "              (\"tracking, falling\", \"-\"),\n",
    "              (\"tracking, time-reversed\", \"--\")]\n",
    "for name, ls in plot_cases:\n",
    "    meas, amp = cases[name][0]\n",
    "    y = meas - meas[nslow // 2]\n",
    "    if name.startswith(\"fixed-bin\"):          # phase valid only near main lobe\n",
    "        y = np.where(amp > 0.3 * amp.max(), y, np.nan)\n",
    "    plt.plot(eta * 1e3, y, lw=1.8, ls=ls, label=name)\n",
    "for rn in (\"flow\", \"fc\", \"fhigh\"):\n",
    "    m = -2 * np.pi * refs[rn] * tau2 + np.pi * K2 * tau2**2\n",
    "    m -= m[nslow // 2]\n",
    "    plt.plot(eta * 1e3, m, \"k:\", lw=1)\n",
    "    plt.annotate(rn, (eta[-1] * 1e3, m[-1]), fontsize=9,\n",
    "                 xytext=(4, 0), textcoords=\"offset points\")\n",
    "plt.xlabel(\"Slow time (ms)\")\n",
    "plt.ylabel(\"Relative phase (rad)\")\n",
    "plt.title(\"Measured phase histories vs reference models (dotted)\")\n",
    "plt.legend(fontsize=9)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "Reading the table: the tracking readout of a rising sweep is\n",
    "`flow`-referenced (RMS 0.000), using `fcenter` there would be a 5 rad RMS error\n",
    "over even this small aperture. The mid-sweep factor makes `fcenter` exact for\n",
    "rising and falling sweeps, which is what lets identical processing handle both\n",
    "sweep directions. Time reversal lands on the last-sample frequency `f_last`, not\n",
    "exactly `fhigh`. The fixed-bin readout is `fcenter`-referenced no matter what,\n",
    "because window leakage re-references it, a property of that readout, not of\n",
    "backprojection."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## Backprojection with a DC-start sweep\n",
    "\n",
    "Now with `torchbp.ops.backprojection_polar_2d` on real image formation. The stress\n",
    "test for point 2: `flow = 0`, B = 1.5 GHz. With first-sample-referenced data the\n",
    "correct carrier is `fc = flow = 0`, i.e. no phase compensation at all.\n",
    "\n",
    "`torchbp.util.generate_fmcw_data` generates data referenced to the `fc` argument\n",
    "it is given (phase $2\\pi(-f_c\\tau - K\\tau t + \\tfrac12 K\\tau^2)$). Real dechirped\n",
    "hardware data is referenced to the first ADC sample, so hardware-convention\n",
    "data is generated by passing the first-sample frequency as the carrier:\n",
    "`flow` for a rising sweep, `fhigh` (with `bw < 0`) for a falling one.\n",
    "\n",
    "Three images of the same point target:\n",
    "\n",
    "* A: `fc = 0`, no re-referencing (consistent conventions),\n",
    "* B: `fc = B/2` without re-referencing (mismatched: the image model assumes a\n",
    "  mid-sweep reference the data do not have),\n",
    "* C: `fc = B/2` with the mid-sweep factor (consistent again).\n",
    "\n",
    "A and C should be identical and focused, B should defocus. The azimuth resolution of A\n",
    "is predicted from the physical spectral centroid (750 MHz, the signal energy\n",
    "is what it is), not from the `fc` parameter.\n",
    "\n",
    "Implementing the mid-sweep factor with the backprojection kernel: at range bin\n",
    "index $s_x$ of the oversampled range profile the beat frequency is $f = s_x\n",
    "f_s/M$, so $e^{-j\\pi f T} = e^{-j (\\pi/\\text{oversample})\\, s_x}$: a pure linear\n",
    "phase across range bins. Multiplying the data with it would shift the data's\n",
    "bin-domain spectrum away from DC and degrade the kernel's linear interpolation.\n",
    "Instead it is folded into `data_fmod`, which the kernel removes exactly at the\n",
    "interpolated (fractional) bin position, equivalent, and interpolation accuracy is\n",
    "preserved. In this processing chain the falling-sweep branch conjugates the data,\n",
    "which flips the sign of the beat axis, hence the `sign(bw)` factor; with a signed\n",
    "beat-frequency axis the factor is always $e^{-j\\pi f T}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "def range_compress(raw, bw, fs, tsweep, oversample, midsweep):\n",
    "    '''Window + oversampled range FFT + spectrum centering, torchbp\n",
    "    conventions. Optionally re-references the phase to mid sweep.'''\n",
    "    nsamples = raw.shape[-1]\n",
    "    M = nsamples * oversample\n",
    "    w = torch.tensor(hamming(nsamples)[None, :], dtype=torch.float32)\n",
    "    data_fmod = -np.pi * (1 - (oversample - 1) / oversample)\n",
    "    if bw > 0:  # with a rising sweep the IF frequencies are negative\n",
    "        data = torch.fft.ifft(raw * w, dim=-1, n=M)\n",
    "    else:\n",
    "        data = torch.fft.ifft(raw.conj() * w, dim=-1, n=M).conj()\n",
    "        data_fmod = -data_fmod\n",
    "    data = data * torch.exp(1j * data_fmod * torch.arange(M))[None, :]\n",
    "    if midsweep:\n",
    "        # exp(-1j*pi*f*tsweep) over the beat axis is a pure linear phase of\n",
    "        # sign(bw)*pi/oversample rad/bin: fold it into data_fmod, which the\n",
    "        # kernel removes exactly at the interpolated bin position.\n",
    "        data_fmod = data_fmod + np.sign(bw) * np.pi / oversample\n",
    "    return data, data_fmod\n",
    "\n",
    "def peak_and_widths(img, grid):\n",
    "    '''Peak value, peak (r, theta) and -3 dB widths of a polar image.'''\n",
    "    a = torch.abs(img)\n",
    "    ir, it = np.unravel_index(torch.argmax(a).item(), a.shape)\n",
    "    r_ax = np.linspace(*grid[\"r\"], grid[\"nr\"], endpoint=False)\n",
    "    t_ax = np.linspace(*grid[\"theta\"], grid[\"ntheta\"], endpoint=False)\n",
    "    def w3db(cut, dx):\n",
    "        aa = cut / cut.max()\n",
    "        i0, lev = int(np.argmax(aa)), 10 ** (-3 / 20)\n",
    "        def cross(d):\n",
    "            i = i0\n",
    "            while 0 < i < len(aa) - 1 and aa[i + d] > lev:\n",
    "                i += d\n",
    "            return i + d * (aa[i] - lev) / (aa[i] - aa[i + d])\n",
    "        return (cross(1) - cross(-1)) * dx\n",
    "    return (a[ir, it].item(), r_ax[ir], t_ax[it],\n",
    "            w3db(a[:, it].numpy(), r_ax[1] - r_ax[0]),\n",
    "            w3db(a[ir, :].numpy(), t_ax[1] - t_ax[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow3 = 0.0 # hypothetical DC-start sweep\n",
    "bw3 = 1.5e9\n",
    "tsweep3 = 100e-6\n",
    "fs3 = 12.8e6\n",
    "oversample = 2\n",
    "fcenter3 = flow3 + bw3 / 2\n",
    "nsweeps = 256\n",
    "\n",
    "pos = torch.zeros([nsweeps, 3], dtype=torch.float32)\n",
    "pos[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * 0.1\n",
    "L = float(pos[-1, 1] - pos[0, 1])\n",
    "target_pos = torch.tensor([[100.0, 0, 0]], dtype=torch.float32)\n",
    "target_rcs = torch.tensor([[1.0]], dtype=torch.float32)\n",
    "\n",
    "grid = {\"r\": (95.0, 105.0), \"theta\": (-0.03, 0.03), \"nr\": 256, \"ntheta\": 256}\n",
    "r_res3 = c0 / (2 * bw3 * oversample)\n",
    "\n",
    "# hardware convention: reference = first ADC sample -> generate with flow\n",
    "raw = torchbp.util.generate_fmcw_data(target_pos, target_rcs, pos,\n",
    "                                      flow3, bw3, tsweep3, fs3)\n",
    "\n",
    "imgs = {}\n",
    "for name, fc_bp, ms in [(\"A: fc = 0 (flow), no midsweep\", 0.0, False),\n",
    "                        (\"B: fc = B/2, no midsweep\", fcenter3, False),\n",
    "                        (\"C: fc = B/2, midsweep\", fcenter3, True)]:\n",
    "    data, data_fmod = range_compress(raw, bw3, fs3, tsweep3, oversample, ms)\n",
    "    imgs[name] = torchbp.ops.backprojection_polar_2d(\n",
    "        data, grid, fc_bp, r_res3, pos, data_fmod=data_fmod).squeeze()\n",
    "\n",
    "names = list(imgs)\n",
    "pkA = peak_and_widths(imgs[names[0]], grid)\n",
    "for name, img in imgs.items():\n",
    "    pk, r, th, rw, tw = peak_and_widths(img, grid)\n",
    "    print(f\"{name:30s} peak {20*np.log10(pk/pkA[0]):+6.2f} dB, \"\n",
    "          f\"r = {r:.2f} m, theta = {th:+.4f}\")\n",
    "    print(f\"   -3 dB widths: range {rw:.3f} m, azimuth {tw*r:.3f} m\")\n",
    "\n",
    "lam_c = c0 / fcenter3     # PHYSICAL spectral centroid, not the fc parameter\n",
    "print(f\"\\nazimuth -3 dB width predicted from physical centroid \"\n",
    "      f\"({fcenter3/1e6:.0f} MHz): {0.886 * lam_c * 100.0 / (2 * L):.3f} m\")\n",
    "print(f\"range -3 dB width predicted (Hamming window): \"\n",
    "      f\"{1.3 * c0 / (2 * bw3):.3f} m\")\n",
    "dmag = (imgs[names[0]].abs() - imgs[names[2]].abs()).abs().max().item()\n",
    "print(f\"max ||A| - |C|| / peak = {dmag/pkA[0]:.1e}  (A and C identical)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"axes.grid\": False}):\n",
    "    fig, axs = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)\n",
    "    extent = [*grid[\"theta\"], *grid[\"r\"]]\n",
    "    for ax, name in zip(axs, names):\n",
    "        db = 20 * torch.log10(imgs[name].abs() / pkA[0] + 1e-12)\n",
    "        im = ax.imshow(db.numpy(), origin=\"lower\", vmin=-20, vmax=0,\n",
    "                       extent=extent, aspect=\"auto\")\n",
    "        ax.set_title(name, fontsize=10)\n",
    "        ax.set_xlabel(\"Angle (sin radians)\")\n",
    "    axs[0].set_ylabel(\"Range (m)\")\n",
    "    fig.colorbar(im, ax=axs, shrink=0.9, label=\"dB rel. image A peak\")\n",
    "\n",
    "plt.figure(figsize=(8, 3.2))\n",
    "t_ax = np.linspace(*grid[\"theta\"], grid[\"ntheta\"], endpoint=False)\n",
    "for name, ls in zip(names, [\"-\", \"-\", \"--\"]):\n",
    "    img = imgs[name]\n",
    "    ir = np.unravel_index(img.abs().argmax().item(), img.shape)[0]\n",
    "    cut = img.abs()[ir, :].numpy() / pkA[0]\n",
    "    plt.plot(t_ax * 100.0, 20 * np.log10(cut + 1e-12),\n",
    "             lw=1.8, ls=ls, label=name)\n",
    "plt.ylim(-45, 2)\n",
    "plt.xlabel(\"Cross-range at 100 m (m)\")\n",
    "plt.ylabel(\"dB rel. image A peak\")\n",
    "plt.title(\"Azimuth cuts through the point target (DC-start sweep)\")\n",
    "plt.legend(fontsize=9)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Image A focuses with `fc = 0`, the phase compensation is unity and the\n",
    "azimuth resolution still matches the prediction from the physical 750 MHz\n",
    "spectral centroid, because focusing information lives in the data. Image B, the mismatched convention, is badly defocused and\n",
    "shifted. Image C is identical to A in magnitude: `fc = B/2` becomes correct once\n",
    "the data are actually re-referenced to mid sweep.\n",
    "\n",
    "Note that A and C are only identical in magnitude: a complex SAR image formed\n",
    "with carrier `fc` carries an `fc`-dependent spatial phase ramp (this is why e.g.\n",
    "`polar_to_cart` and interferometric processing need the same `fc` the image was\n",
    "formed with)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "## One processing chain for rising and falling sweeps\n",
    "\n",
    "The practical payoff of mid-sweep re-referencing: with `fc = fcenter` the same\n",
    "processing works for both sweep directions. A realistic 6 GHz / 200 MHz pair,\n",
    "without re-referencing the rising sweep would need `fc = flow` and the falling\n",
    "sweep `fc = fhigh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc_rf = 6e9\n",
    "bw4 = 200e6\n",
    "fs4 = 2e6\n",
    "tsweep4 = 100e-6\n",
    "pos4 = torch.zeros([nsweeps, 3], dtype=torch.float32)\n",
    "pos4[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * 0.25 * c0 / fc_rf\n",
    "grid4 = {\"r\": (95.0, 105.0), \"theta\": (-0.03, 0.03), \"nr\": 256, \"ntheta\": 256}\n",
    "r_res4 = c0 / (2 * bw4 * oversample)\n",
    "\n",
    "# hardware convention: carrier = frequency at the first ADC sample\n",
    "raw_rise = torchbp.util.generate_fmcw_data(\n",
    "    target_pos, target_rcs, pos4, fc_rf - bw4 / 2, bw4, tsweep4, fs4)\n",
    "raw_fall = torchbp.util.generate_fmcw_data(\n",
    "    target_pos, target_rcs, pos4, fc_rf + bw4 / 2, -bw4, tsweep4, fs4)\n",
    "\n",
    "res = {}\n",
    "for name, raw4, bws in [(\"rising\", raw_rise, bw4), (\"falling\", raw_fall, -bw4)]:\n",
    "    data, data_fmod = range_compress(raw4, bws, fs4, tsweep4, oversample, True)\n",
    "    res[name] = torchbp.ops.backprojection_polar_2d(\n",
    "        data, grid4, fc_rf, r_res4, pos4, data_fmod=data_fmod).squeeze()\n",
    "    pk, r, th, rw, tw = peak_and_widths(res[name], grid4)\n",
    "    print(f\"{name:8s} sweep, fc = fcenter + midsweep: \"\n",
    "          f\"r = {r:.2f} m, theta = {th:+.4f}\")\n",
    "    print(f\"   widths: range {rw:.3f} m, azimuth {tw*r:.3f} m\")\n",
    "\n",
    "d = (res[\"rising\"].abs() - res[\"falling\"].abs()).abs().max().item()\n",
    "print(f\"\\nmax ||rising| - |falling|| / peak = \"\n",
    "      f\"{d / res['rising'].abs().max().item():.1e}\")\n",
    "pr = res[\"rising\"].flatten()[res[\"rising\"].abs().argmax()].item()\n",
    "pf = res[\"falling\"].flatten()[res[\"falling\"].abs().argmax()].item()\n",
    "print(f\"peak phase: rising {np.angle(pr):+.3f} rad, \"\n",
    "      f\"falling {np.angle(pf):+.3f} rad\")\n",
    "print(\"   (difference is 2x the residual video phase; opposite signs \"\n",
    "      \"for opposite sweeps)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Pulse compression radar\n",
    "\n",
    "A pulse compression radar transmits a coded pulse, receives the full bandwidth\n",
    "and range compresses by correlating with a replica of the pulse (equivalently\n",
    "FFT, multiply by the conjugated reference spectrum, IFFT). Typically the data\n",
    "entering range compression is at DC-centered baseband: quadrature demodulation\n",
    "with the LO at the band centre, so the pulse spectrum spans $\\pm B/2$ around\n",
    "$\\mathrm{DC}$.\n",
    "\n",
    "The structural difference to dechirp processing is that for correlation delay changes range origin instead of affecting the phase.\n",
    "What sets the phase reference instead is carrier\n",
    "removal: demodulating an echo delayed by $\\tau$ with $e^{-j2\\pi f_d t}$ leaves\n",
    "the compressed peak with phase $-2\\pi f_d \\tau$, so backprojection needs\n",
    "`fc` $= f_d$, the net demodulation frequency, which is the band center for\n",
    "DC-centered data. There is also no residual video phase with matched filtering.\n",
    "\n",
    "The same phase-history experiment as earlier, on a wide-fractional-bandwidth\n",
    "pulse radar (0.5–1.5 GHz chirp pulse), with the manipulations that changed the\n",
    "dechirp reference:\n",
    "\n",
    "* record window delayed: the target moves to a smaller lag (a range-origin\n",
    "  shift), the peak phase is unchanged,\n",
    "* ADC samples time-reversed (with the reference reversed consistently): the\n",
    "  correlation output is exactly the same, compare with FMCW where reversal moved\n",
    "  the dechirp reference to $f_{last}$,\n",
    "* LO not at the band centre (here at the band edge, band $0..B$): the\n",
    "  reference follows the LO, not the band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc5, bw5, tp5 = 1.0e9, 1.0e9, 0.5e-6   # RF band 0.5-1.5 GHz, 0.5 us pulse\n",
    "fs5 = 3e9\n",
    "ntp5 = int(round(tp5 * fs5))\n",
    "nrec5 = int(round(1.0e-6 * fs5))\n",
    "M5 = 8192 # >= nrec5 + ntp5: linear correlation\n",
    "fax5 = np.fft.fftfreq(M5, 1 / fs5)\n",
    "\n",
    "v5, R05, L5 = 50.0, 50.0, 12.0\n",
    "nslow5 = 128\n",
    "eta5 = np.linspace(-L5 / (2 * v5), L5 / (2 * v5), nslow5)\n",
    "tau5 = 2 * np.sqrt((v5 * eta5)**2 + R05**2) / c0\n",
    "\n",
    "def pulse_env(tt, tp, bw):\n",
    "    '''Complex envelope of the TX pulse: DC-centred chirp, -B/2 to +B/2.'''\n",
    "    K = bw / tp\n",
    "    return np.where((tt >= 0) & (tt < tp),\n",
    "                    np.exp(2j * np.pi * (-bw / 2 * tt + 0.5 * K * tt**2)), 0)\n",
    "\n",
    "def compressed_peak_phase(f_demod, t0=0.0, reverse=False):\n",
    "    '''Demodulate the echoes at f_demod, range compress by correlation with\n",
    "    the demodulated replica and evaluate the correlation spectrum at the\n",
    "    exact target lag, which is what backprojection's interpolation does.'''\n",
    "    t_win = t0 + np.arange(nrec5) / fs5\n",
    "    t_ref = np.arange(ntp5) / fs5\n",
    "    foff = fc5 - f_demod              # residual carrier after demodulation\n",
    "    out = []\n",
    "    for tau in tau5:\n",
    "        s = (pulse_env(t_win - tau, tp5, bw5)\n",
    "             * np.exp(2j * np.pi * foff * (t_win - tau))\n",
    "             * np.exp(-2j * np.pi * f_demod * tau))\n",
    "        ref = pulse_env(t_ref, tp5, bw5) * np.exp(2j * np.pi * foff * t_ref)\n",
    "        lag = tau - t0\n",
    "        if reverse:\n",
    "            s, ref = s[::-1], ref[::-1]\n",
    "            lag = (nrec5 - 1) / fs5 - (ntp5 - 1) / fs5 - (tau - t0)\n",
    "        Y = np.fft.fft(s, M5) * np.conj(np.fft.fft(ref, M5))\n",
    "        out.append(np.sum(Y * np.exp(2j * np.pi * fax5 * lag)) / M5)\n",
    "    return np.unwrap(np.angle(np.array(out)))\n",
    "\n",
    "cases5 = {\n",
    "    \"DC-centred baseband (LO = band centre)\":      compressed_peak_phase(fc5),\n",
    "    \"record window delayed 200 ns\":                compressed_peak_phase(fc5, t0=200e-9),\n",
    "    \"ADC samples time-reversed\":                   compressed_peak_phase(fc5, reverse=True),\n",
    "    \"LO at band edge (band 0..B, not DC-centred)\": compressed_peak_phase(fc5 - bw5 / 2),\n",
    "}\n",
    "refs5 = {\"flow\": fc5 - bw5 / 2, \"fc\": fc5, \"fhigh\": fc5 + bw5 / 2}\n",
    "\n",
    "print(\"peak-phase-history RMS mismatch (rad) vs reference model \"\n",
    "      f\"(fc vs flow span = {2*np.pi*(bw5/2)*(tau5.max()-tau5.min()):.1f} rad):\")\n",
    "print(f\"{'case':46s}\" + \"\".join(f\"{r:>9s}\" for r in refs5) + \"   winner\")\n",
    "for name, meas in cases5.items():\n",
    "    r = {}\n",
    "    for rn, fr in refs5.items():\n",
    "        d = meas - (-2 * np.pi * fr * tau5)\n",
    "        d = d - d[nslow5 // 2]        # constants are unobservable\n",
    "        r[rn] = np.sqrt(np.mean(d**2))\n",
    "    row = \"\".join(f\"{r[rn]:9.3f}\" for rn in refs5)\n",
    "    print(f\"{name:46s}{row}   {min(r, key=r.get)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "The winner is the demodulation frequency in every row: the time-origin and\n",
    "sample-order manipulations are invisible to the correlation readout.\n",
    "\n",
    "## Backprojection of pulse-compressed data\n",
    "\n",
    "Same radar, now with `torchbp.ops.backprojection_polar_2d`. The range bin spacing is\n",
    "`r_res` $= c/(2 f_s \\cdot \\text{oversample})$ (oversampling by zero-padding\n",
    "the DC-centred product spectrum), bin 0 is at zero range, and the compressed\n",
    "data is used directly: the peak phase $-4\\pi f_d R/c$ already matches the\n",
    "torchbp convention, no conjugation needed. Four images of the same point target:\n",
    "\n",
    "* A: DC-centred data, `fc = fcenter` (consistent conventions),\n",
    "* B: same data, `fc = flow` (wrong bookkeeping, 500 MHz off),\n",
    "* C: ADC window delayed by 250 ns, `fc = fcenter` unchanged: the delay is a\n",
    "  pure range-origin shift, corrected with `d0` $= -c t_0/2$, and the image\n",
    "  should match A exactly, including phase. A dechirp window shift would\n",
    "  instead have changed the required `fc` by $K t_0$.\n",
    "* D: LO offset 50 MHz below the band centre. The compressed profile then\n",
    "  rides a residual carrier $e^{j2\\pi f_{off}(l - \\tau)}$ across the lag bins\n",
    "  $l$. Folding it into `data_fmod` keeps the kernel's interpolation accurate and moves the\n",
    "  $-2\\pi f_{off}\\tau$ part of the phase into the compensation, so the total\n",
    "  reference is $-2\\pi(f_d + f_{off})\\tau$ and `fc = fcenter` again: the\n",
    "  bookkeeping mirror image of mid-sweep re-referencing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs5b = 1.2e9\n",
    "ntp5b = int(round(tp5 * fs5b))\n",
    "nrec5b = int(round(1.5e-6 * fs5b))\n",
    "M5b = 4096\n",
    "nkeep5 = 2048\n",
    "r_res5 = c0 / (2 * fs5b * oversample)\n",
    "\n",
    "pos5 = torch.zeros([nsweeps, 3], dtype=torch.float32)\n",
    "pos5[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * 0.1\n",
    "L5b = float(pos5[-1, 1] - pos5[0, 1])\n",
    "tau5b = 2 * np.sqrt(100.0**2 + pos5[:, 1].numpy()**2) / c0\n",
    "grid5 = {\"r\": (95.0, 105.0), \"theta\": (-0.06, 0.06), \"nr\": 256, \"ntheta\": 256}\n",
    "\n",
    "def make_pulse_data(t0=0.0, f_demod=fc5):\n",
    "    '''Echoes demodulated at f_demod; the record starts t0 after TX.'''\n",
    "    t_win = t0 + np.arange(nrec5b) / fs5b\n",
    "    foff = fc5 - f_demod\n",
    "    return (pulse_env(t_win[None, :] - tau5b[:, None], tp5, bw5)\n",
    "            * np.exp(2j * np.pi * foff * (t_win[None, :] - tau5b[:, None]))\n",
    "            * np.exp(-2j * np.pi * f_demod * tau5b)[:, None])\n",
    "\n",
    "def range_compress_pc(s, foff=0.0):\n",
    "    '''Correlation range compression with a Hamming band taper and\n",
    "    oversampled (zero-padded product spectrum) output.'''\n",
    "    t_ref = np.arange(ntp5b) / fs5b\n",
    "    ref = pulse_env(t_ref, tp5, bw5) * np.exp(2j * np.pi * foff * t_ref)\n",
    "    fax = np.fft.fftfreq(M5b, 1 / fs5b)\n",
    "    W = ((0.54 + 0.46 * np.cos(2 * np.pi * (fax - foff) / bw5))\n",
    "         * (np.abs(fax - foff) <= bw5 / 2))\n",
    "    Y = (np.fft.fft(s, M5b, axis=-1)\n",
    "         * np.conj(np.fft.fft(ref, M5b))[None, :] * W[None, :])\n",
    "    Yp = np.zeros((s.shape[0], oversample * M5b), complex)\n",
    "    Yp[:, :M5b // 2] = Y[:, :M5b // 2]\n",
    "    Yp[:, -(M5b // 2):] = Y[:, M5b // 2:]\n",
    "    y = np.fft.ifft(Yp, axis=-1) * oversample\n",
    "    return torch.tensor(y[:, :nkeep5], dtype=torch.complex64)\n",
    "\n",
    "t0C = 250e-9                              # delayed ADC window for image C\n",
    "foffD = 50e6                              # LO 50 MHz below band centre for D\n",
    "data_pc = range_compress_pc(make_pulse_data())\n",
    "data_pcC = range_compress_pc(make_pulse_data(t0C))\n",
    "data_pcD = range_compress_pc(make_pulse_data(f_demod=fc5 - foffD), foff=foffD)\n",
    "fmodD = 2 * np.pi * foffD / (fs5b * oversample)  # residual carrier, rad/bin\n",
    "\n",
    "imgs5 = {}\n",
    "for name, dat, fc_bp, d0, fmod in [\n",
    "        (\"A: fc = fcenter\", data_pc, fc5, 0.0, 0.0),\n",
    "        (\"B: fc = flow\", data_pc, fc5 - bw5 / 2, 0.0, 0.0),\n",
    "        (\"C: window +250 ns, fc = fcenter\", data_pcC, fc5, -c0 * t0C / 2, 0.0),\n",
    "        (\"D: LO offset -50 MHz, data_fmod\", data_pcD, fc5, 0.0, fmodD)]:\n",
    "    imgs5[name] = torchbp.ops.backprojection_polar_2d(\n",
    "        dat, grid5, fc_bp, r_res5, pos5, d0=d0, data_fmod=fmod).squeeze()\n",
    "\n",
    "names5 = list(imgs5)\n",
    "pkA5 = peak_and_widths(imgs5[names5[0]], grid5)\n",
    "for name, img in imgs5.items():\n",
    "    pk, r, th, rw, tw = peak_and_widths(img, grid5)\n",
    "    print(f\"{name:32s} peak {20*np.log10(pk/pkA5[0]):+6.2f} dB, \"\n",
    "          f\"r = {r:.2f} m, theta = {th:+.4f}\")\n",
    "    print(f\"   -3 dB widths: range {rw:.3f} m, azimuth {tw*r:.3f} m\")\n",
    "\n",
    "lam5 = c0 / fc5\n",
    "print(f\"\\nazimuth -3 dB width predicted from the 1 GHz band centre: \"\n",
    "      f\"{0.886 * lam5 * 100.0 / (2 * L5b):.3f} m\")\n",
    "print(f\"range -3 dB width predicted (Hamming band taper): \"\n",
    "      f\"{1.3 * c0 / (2 * bw5):.3f} m\")\n",
    "dAC = (imgs5[names5[0]] - imgs5[names5[2]]).abs().max().item()\n",
    "print(f\"max |A - C| / peak = {dAC/pkA5[0]:.1e}\")\n",
    "print(\"   (complex diff: window delay moved range origin, not phase ref)\")\n",
    "dAD = (imgs5[names5[0]] - imgs5[names5[3]]).abs().max().item()\n",
    "print(f\"max |A - D| / peak = {dAD/pkA5[0]:.1e}\")\n",
    "print(\"   (small linear-interpolation residual from carrier-riding bins)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"axes.grid\": False}):\n",
    "    fig, axs = plt.subplots(1, 4, figsize=(13, 3.2), sharey=True)\n",
    "    extent = [*grid5[\"theta\"], *grid5[\"r\"]]\n",
    "    for ax, name in zip(axs, names5):\n",
    "        db = 20 * torch.log10(imgs5[name].abs() / pkA5[0] + 1e-12)\n",
    "        im = ax.imshow(db.numpy(), origin=\"lower\", vmin=-20, vmax=0,\n",
    "                       extent=extent, aspect=\"auto\")\n",
    "        ax.set_title(name, fontsize=9)\n",
    "        ax.set_xlabel(\"Angle (sin radians)\")\n",
    "    axs[0].set_ylabel(\"Range (m)\")\n",
    "    fig.colorbar(im, ax=axs, shrink=0.9, label=\"dB rel. image A peak\")\n",
    "\n",
    "plt.figure(figsize=(8, 3.2))\n",
    "t_ax = np.linspace(*grid5[\"theta\"], grid5[\"ntheta\"], endpoint=False)\n",
    "for name, ls in zip(names5, [\"-\", \"-\", \"--\", \":\"]):\n",
    "    img = imgs5[name]\n",
    "    ir = np.unravel_index(img.abs().argmax().item(), img.shape)[0]\n",
    "    cut = img.abs()[ir, :].numpy() / pkA5[0]\n",
    "    plt.plot(t_ax * 100.0, 20 * np.log10(cut + 1e-12),\n",
    "             lw=1.8, ls=ls, label=name)\n",
    "plt.ylim(-45, 2)\n",
    "plt.xlabel(\"Cross-range at 100 m (m)\")\n",
    "plt.ylabel(\"dB rel. image A peak\")\n",
    "plt.title(\"Azimuth cuts through the point target (pulse compression)\")\n",
    "plt.legend(fontsize=9)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "Image A focuses with the textbook `fc = fcenter` and matches the predicted\n",
    "widths. Image B shows that `fc` with 500 MHz\n",
    "error defocuses azimuth. C matches A\n",
    "to numerical precision in complex value, so the record window placement\n",
    "affects only the range origin, and D shows the offset-LO case restored to\n",
    "`fc = fcenter` by `data_fmod`.\n",
    "\n",
    "Summary of the two waveform families: dechirp readouts inherit their phase\n",
    "reference from time-origin conventions (a single FFT maps delay to frequency,\n",
    "so linear-phase and origin conventions imprint on the phase vs range) and the\n",
    "required `fc` is therefore implementation-sensitive: `flow`, `fcenter`,\n",
    "`fhigh` or none of them, depending on window placement, sample order and\n",
    "re-referencing factors. Matched-filter readouts inherit the reference from the\n",
    "demodulation chain alone: `fc` equals the net frequency removed from the echo\n",
    "and is insensitive to windowing, trigger timing and sample order. With the\n",
    "usual DC-centered baseband convention that is `fcenter`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}