From bb95a0eda85674d47831292b282b0cf3aa80017f Mon Sep 17 00:00:00 2001 From: chenzc24 <905906303@qq.com> Date: Sun, 28 Jun 2026 17:03:40 +0800 Subject: [PATCH] Fix polar harmonic annotations --- .../exp_s11_polar_memory_effect.py | 6 +- python/src/adctoolbox/spectrum/_harmonics.py | 12 ++-- .../spectrum/analyze_spectrum_polar.py | 1 + .../spectrum/plot_spectrum_polar.py | 38 +++++----- .../unit/spectrum/test_plot_spectrum_polar.py | 72 +++++++++++++++++++ 5 files changed, 99 insertions(+), 30 deletions(-) create mode 100644 python/tests/unit/spectrum/test_plot_spectrum_polar.py diff --git a/python/src/adctoolbox/examples/02_spectrum/exp_s11_polar_memory_effect.py b/python/src/adctoolbox/examples/02_spectrum/exp_s11_polar_memory_effect.py index 613452b3..09091da8 100644 --- a/python/src/adctoolbox/examples/02_spectrum/exp_s11_polar_memory_effect.py +++ b/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() diff --git a/python/src/adctoolbox/spectrum/_harmonics.py b/python/src/adctoolbox/spectrum/_harmonics.py index a3779a36..22c36e90 100644 --- a/python/src/adctoolbox/spectrum/_harmonics.py +++ b/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) diff --git a/python/src/adctoolbox/spectrum/analyze_spectrum_polar.py b/python/src/adctoolbox/spectrum/analyze_spectrum_polar.py index 397138f8..fa9322ea 100644 --- a/python/src/adctoolbox/spectrum/analyze_spectrum_polar.py +++ b/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 --- diff --git a/python/src/adctoolbox/spectrum/plot_spectrum_polar.py b/python/src/adctoolbox/spectrum/plot_spectrum_polar.py index a5e86bf7..01c57615 100644 --- a/python/src/adctoolbox/spectrum/plot_spectrum_polar.py +++ b/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]) \ No newline at end of file + ax.set_ylim([0, max_radius]) diff --git a/python/tests/unit/spectrum/test_plot_spectrum_polar.py b/python/tests/unit/spectrum/test_plot_spectrum_polar.py new file mode 100644 index 00000000..0dcb856e --- /dev/null +++ b/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