-
Notifications
You must be signed in to change notification settings - Fork 0
Description
Bug summary
Windows correction is not correct in mlab._spectral_helper:
https://github.com/matplotlib/matplotlib/blob/3418bada1c1f44da1f73916c5603e3ae79fe58c1/lib/matplotlib/mlab.py#L423-L430
The np.abs is not needed, and gives wrong results for windows with negative values, such as flattop. For reference, SciPy's implementation:
https://github.com/scipy/scipy/blob/d9f75db82fdffef06187c9d8d2f0f5b36c7a791b/scipy/signal/_spectral_py.py#L1854-L1859
Code for reproduction
import numpy as np
from matplotlib import mlab
Fs = 1024.0
NFFT = 512
f0 = 128.0
t = np.arange(NFFT) / Fs
x = np.sin(2*np.pi*f0*t)
# Window with negative samples but positive sum
# (simple cosine-shift window for demonstration)
def negative_value_window(N):
n = np.arange(N, dtype=float)
return 0.5 + np.cos(2*np.pi * n / (N - 1))
w = negative_value_window(NFFT)
# Compute spectrum and density
spec_density, f = mlab.psd(x=x, NFFT=NFFT, Fs=Fs, sides='onesided',
window=w, scale_by_freq=True)
spec_spectrum, _ = mlab.psd(x=x, NFFT=NFFT, Fs=Fs, sides='onesided',
window=w, scale_by_freq=False)
# Expected relationship (SciPy semantics):
# spec_density * sum(w**2) == spec_spectrum / Fs * (sum(w))**2
lhs = spec_density * (w**2).sum()
rhs = spec_spectrum / Fs * (w.sum()**2)
print("max abs(lhs - rhs):", np.max(np.abs(lhs - rhs)))
# Current (pre-fix): non-zero discrepancy due to using sum(abs(w))**2 instead of (sum(w))**2
# Expected (post-fix): discrepancy ~ numerical noise (~1e-8 to 1e-6)Actual outcome
The relationship above does not hold for windows with negative samples because the code divides by sum(abs(w))**2 instead of (sum(w))**2 under spectrum scaling (scale_by_freq=False).
Expected outcome
Under spectrum scaling (scale_by_freq=False), PSD should be normalized by (sum(window))**2, not sum(abs(window))**2. Under spectral density scaling (scale_by_freq=True), PSD should be normalized by Fs * sum(window**2).
Additional information
Researcher specification:
- File: lib/matplotlib/mlab.py
- Function:
_spectral_helper - Change under mode == 'psd' and scale_by_freq == False from:
result /= np.abs(window).sum()**2
to:result /= (window.sum())**2
- Rationale: For spectrum scaling, normalization by
(sum(w))**2preserves segment power and matches SciPy'swelch(scale='spectrum'). Usingsum(abs(w))**2biases windows with negative coefficients. - Density scaling (scale_by_freq=True) already correctly uses
Fs * sum(w**2)and should remain unchanged. - Magnitude/phase/angle modes remain unchanged.
Tests to add (pytest):
- Spectrum scaling uses
(sum(w))**2with a real window that includes negative samples. - Spectral density scaling uses
Fs * sum(w**2). - Regression:
spec_density * sum(w**2) == spec_spectrum / Fs * (sum(w))**2.
Operating system
N/A
Matplotlib Version
From branch: matplotlib__matplotlib-25122
Installation
N/A
Task ref: 275