Fix polar harmonic annotations

python/src/adctoolbox/examples/02_spectrum/exp_s11_polar_memory_effect.py

@@ -60,21 +60,21 @@
 signal_1 = sig_ideal + k3 * sig_ideal**3 + DC + np.random.randn(N) * base_noise
 plt.sca(axes[0, 0])
 result_1 = analyze_spectrum_polar(signal_1, fs=Fs, fixed_radial_range=120)
-axes[0, 0].set_title(f'HD3={hd3_dB}dB, k3>0\n(Thermal Noise: 500 uVrms)', pad=20, fontsize=12, fontweight='bold')
+axes[0, 0].set_title(f'HD3={hd3_dB}dB, k3>0\n(Thermal Noise: {base_noise*1e6:.0f} uVrms)', pad=20, fontsize=12, fontweight='bold')
 print(f"[HD3={hd3_dB}dB, k3>0] SNDR={result_1['sndr_dbc']:.2f}dB, THD={result_1['thd_dbc']:.2f}dB, HD3={result_1['harmonics_dbc'][1]:.2f}dB")
 
 # Case 2: HD3 only, k3 negative
 signal_2 = sig_ideal - k3 * sig_ideal**3 + DC + np.random.randn(N) * base_noise
 plt.sca(axes[0, 1])
 result_2 = analyze_spectrum_polar(signal_2, fs=Fs, fixed_radial_range=120)
-axes[0, 1].set_title(f'HD3={hd3_dB}dB, k3<0\n(Thermal Noise: 500 uVrms)', pad=20, fontsize=12, fontweight='bold')
+axes[0, 1].set_title(f'HD3={hd3_dB}dB, k3<0\n(Thermal Noise: {base_noise*1e6:.0f} uVrms)', pad=20, fontsize=12, fontweight='bold')
 print(f"[HD3={hd3_dB}dB, k3<0] SNDR={result_2['sndr_dbc']:.2f}dB, THD={result_2['thd_dbc']:.2f}dB, HD3={result_2['harmonics_dbc'][1]:.2f}dB")
 
 # Case 3: HD2 + HD3 combined
 signal_3 = sig_ideal + k2 * sig_ideal**2 - k3 * sig_ideal**3 + DC + np.random.randn(N) * base_noise
 plt.sca(axes[0, 2])
 result_3 = analyze_spectrum_polar(signal_3, fs=Fs, fixed_radial_range=120)
-axes[0, 2].set_title(f'HD2={hd2_dB}dB + HD3={hd3_dB}dB\n(Thermal Noise: 500 uVrms)', pad=20, fontsize=12, fontweight='bold')
+axes[0, 2].set_title(f'HD2={hd2_dB}dB + HD3={hd3_dB}dB\n(Thermal Noise: {base_noise*1e6:.0f} uVrms)', pad=20, fontsize=12, fontweight='bold')
 print(f"[HD2+HD3] SNDR={result_3['sndr_dbc']:.2f}dB, THD={result_3['thd_dbc']:.2f}dB, HD2={result_3['harmonics_dbc'][0]:.2f}dB, HD3={result_3['harmonics_dbc'][1]:.2f}dB")
 
 print()

python/src/adctoolbox/spectrum/_harmonics.py

@@ -259,16 +259,16 @@ def _calculate_harmonic_phases(
     hd2_phase_deg = 0.0
     hd3_phase_deg = 0.0
 
-    # Extract HD2 phase (harmonic 2, index 1)
-    if len(harmonic_bins) > 1:
-        hd2_bin = harmonic_bins[1]
+    # Extract HD2 phase (harmonic 2, index 0)
+    if len(harmonic_bins) > 0:
+        hd2_bin = harmonic_bins[0]
         if hd2_bin < len(spec_coherent_normalized):
             hd2_phase_rad = np.angle(spec_coherent_normalized[hd2_bin])
             hd2_phase_deg = np.degrees(hd2_phase_rad)
 
-    # Extract HD3 phase (harmonic 3, index 2)
-    if len(harmonic_bins) > 2:
-        hd3_bin = harmonic_bins[2]
+    # Extract HD3 phase (harmonic 3, index 1)
+    if len(harmonic_bins) > 1:
+        hd3_bin = harmonic_bins[1]
         if hd3_bin < len(spec_coherent_normalized):
             hd3_phase_rad = np.angle(spec_coherent_normalized[hd3_bin])
             hd3_phase_deg = np.degrees(hd3_phase_rad)

python/src/adctoolbox/spectrum/analyze_spectrum_polar.py

@@ -54,6 +54,7 @@ def analyze_spectrum_polar(
         fs=fs,
         win_type=win_type,
         coherent_averaging=True,
+        max_harmonic=max(3, harmonic),
     )
 
     # 2. --- Optional Plotting ---

python/src/adctoolbox/spectrum/plot_spectrum_polar.py

@@ -7,8 +7,6 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from adctoolbox.spectrum._harmonics import _locate_harmonic_bins
-from adctoolbox.spectrum._harmonics import _calculate_harmonic_phases
 
 def plot_spectrum_polar(analysis_results, show_metrics=True, harmonic=5, fixed_radial_range=None, ax=None):
     """
@@ -27,19 +25,14 @@ def plot_spectrum_polar(analysis_results, show_metrics=True, harmonic=5, fixed_r
 
     spec = plot_data['complex_spectrum']
     bin_idx = plot_data['fundamental_bin']
-    fundamental_bin_fractional = plot_data['fundamental_bin_fractional']
     harmonic_bins = plot_data['harmonic_bins']
     collided_harmonics = plot_data.get('collided_harmonics', [])
-    N_fft = analysis_results['N']
 
     # Calculate noise floor for polar plot (1st percentile)
     mag_db = 20 * np.log10(np.abs(spec) + 1e-20)
     minR_dB = np.percentile(mag_db, 1)
     minR_dB = -100 if np.isinf(minR_dB) else minR_dB
 
-    # Calculate harmonic phases
-    hd2_phase_deg, hd3_phase_deg = _calculate_harmonic_phases(spec, harmonic_bins)
-
     # Setup axes
     if ax is None:
         ax = plt.gca()
@@ -95,14 +88,13 @@ def plot_spectrum_polar(analysis_results, show_metrics=True, harmonic=5, fixed_r
         ax.plot([0, phase[bin_idx]], [0, mag[bin_idx]], 'b-', linewidth=2, zorder=10)
 
     # Mark harmonics
-    for h in range(2, harmonic + 1):
+    max_harmonic_to_plot = min(harmonic, len(harmonic_bins) + 1)
+    for h in range(2, max_harmonic_to_plot + 1):
         # Skip if this harmonic collides with fundamental (collision makes plotting meaningless)
         if h in collided_harmonics:
             continue
 
-        harmonic_bin = (bin_idx * h) % N_fft
-        if harmonic_bin > N_fft // 2:
-            harmonic_bin = N_fft - harmonic_bin
+        harmonic_bin = int(harmonic_bins[h - 2])
 
         # Skip if this harmonic aliases to DC (bin 0)
         if harmonic_bin == 0:
@@ -119,13 +111,18 @@ def plot_spectrum_polar(analysis_results, show_metrics=True, harmonic=5, fixed_r
 
     # Add metrics annotation
     if show_metrics and metrics:
-        # Extract HD2 and HD3 from harmonic_powers array
-        harmonic_powers = metrics.get('harmonic_powers', [])
-        hd2_db = harmonic_powers[0] if len(harmonic_powers) > 0 else -150
-        hd3_db = harmonic_powers[1] if len(harmonic_powers) > 1 else -150
-
-        hd2_str = f"HD2 = {hd2_db:7.2f} dB ∠{hd2_phase_deg:6.1f}°"
-        hd3_str = f"HD3 = {hd3_db:7.2f} dB ∠{hd3_phase_deg:6.1f}°"
+        harmonics_dbc = metrics.get('harmonics_dbc', [])
+        harmonic_lines = []
+        for h in (2, 3):
+            harmonic_idx = h - 2
+            harmonic_db = harmonics_dbc[harmonic_idx] if len(harmonics_dbc) > harmonic_idx else -150
+            harmonic_phase_deg = 0.0
+            if harmonic_db > -149 and h not in collided_harmonics and len(harmonic_bins) > harmonic_idx:
+                harmonic_bin = int(harmonic_bins[harmonic_idx])
+                if 0 < harmonic_bin < len(spec):
+                    harmonic_phase_deg = np.degrees(np.angle(spec[harmonic_bin]))
+            harmonic_lines.append(f"HD{h} = {harmonic_db:7.2f} dB ∠{harmonic_phase_deg:6.1f}°")
+        harmonic_text = "\n".join(harmonic_lines)
 
         # Build collision warning if present
         collision_warning = ""
@@ -136,12 +133,11 @@ def plot_spectrum_polar(analysis_results, show_metrics=True, harmonic=5, fixed_r
         metrics_text = (
             f"SNR = {metrics.get('snr_dbc', 0):7.2f} dB\n"
             f"THD = {metrics.get('thd_dbc', 0):7.2f} dB\n"
-            f"{hd2_str}\n"
-            f"{hd3_str}"
+            f"{harmonic_text}"
             f"{collision_warning}"
         )
         ax.text(0.02, 0.02, metrics_text, transform=ax.transAxes, fontsize=10,
                verticalalignment='bottom', family='monospace',
                bbox=dict(boxstyle='round', facecolor='white', edgecolor='gray'))
 
-    ax.set_ylim([0, max_radius])
+    ax.set_ylim([0, max_radius])

python/tests/unit/spectrum/test_plot_spectrum_polar.py

@@ -0,0 +1,72 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pytest
+
+from adctoolbox.spectrum._harmonics import _calculate_harmonic_phases
+from adctoolbox.spectrum.analyze_spectrum_polar import analyze_spectrum_polar
+from adctoolbox.spectrum.compute_spectrum import compute_spectrum
+from adctoolbox.spectrum.plot_spectrum_polar import plot_spectrum_polar
+
+
+def test_calculate_harmonic_phases_uses_hd2_hd3_bins():
+    spec = np.ones(16, dtype=complex)
+    harmonic_bins = np.array([2, 3, 4])
+    spec[2] = np.exp(1j * np.deg2rad(20.0))
+    spec[3] = np.exp(1j * np.deg2rad(-45.0))
+    spec[4] = np.exp(1j * np.deg2rad(120.0))
+
+    hd2_phase, hd3_phase = _calculate_harmonic_phases(spec, harmonic_bins)
+
+    assert hd2_phase == pytest.approx(20.0)
+    assert hd3_phase == pytest.approx(-45.0)
+
+
+def test_plot_spectrum_polar_metrics_text_uses_harmonics_dbc():
+    n_fft = 2**13
+    fs = 100e6
+    amplitude = 0.5
+    fin = 123 / n_fft * fs
+    n = np.arange(n_fft)
+    t = n / fs
+
+    hd2_target = -80.0
+    hd3_target = -66.0
+    hd2_amp = 10 ** (hd2_target / 20)
+    hd3_amp = 10 ** (hd3_target / 20)
+    k2 = hd2_amp / (amplitude / 2)
+    k3 = hd3_amp / (amplitude**2 / 4)
+
+    sine = amplitude * np.sin(2 * np.pi * fin * t)
+    signal = sine + k2 * sine**2 + k3 * sine**3
+    result = compute_spectrum(
+        signal,
+        fs=fs,
+        win_type="rectangular",
+        side_bin=0,
+        max_harmonic=5,
+        coherent_averaging=True,
+    )
+
+    fig, ax = plt.subplots(subplot_kw={"projection": "polar"})
+    plot_spectrum_polar(result, harmonic=5, fixed_radial_range=120, ax=ax)
+    metrics_text = next(text.get_text() for text in ax.texts if "SNR =" in text.get_text())
+    plt.close(fig)
+
+    assert "HD2 =  -80." in metrics_text
+    assert "HD3 =  -66." in metrics_text
+    assert "HD4 =" not in metrics_text
+    assert "HD5 =" not in metrics_text
+    assert "HD2 = -150." not in metrics_text
+    assert "HD3 = -150." not in metrics_text
+
+
+def test_analyze_spectrum_polar_harmonic_controls_compute_depth():
+    n_fft = 2**13
+    fs = 100e6
+    fin = 123 / n_fft * fs
+    n = np.arange(n_fft)
+    signal = 0.5 * np.sin(2 * np.pi * fin * n / fs)
+
+    metrics = analyze_spectrum_polar(signal, fs=fs, harmonic=10, create_plot=False)
+
+    assert len(metrics["harmonics_dbc"]) == 9