Skip to content

Commit adecc56

Browse files
cdholmescdholmes
authored
Allow setting frequency units for psd plots (#32211)
In psd, allow user to specify the units for the sampling frequency of the analyzed array. psd previously assumed units to be Hz and displayed this on the y-axis label. Users can now specify other units. This change has no effect on the psd calculation. Update psd documentation, example, type hinting, release notes
1 parent af04c71 commit adecc56

6 files changed

Lines changed: 54 additions & 10 deletions

File tree

‎doc/release/next_whats_new/new_psd_feature.rst‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
Sampling frequency units can be specified for `.Axes.psd`
2+
---------------------------------------------------------
3+
4+
When creating a power spectral density (psd) plot, the units of the
5+
sampling frequency can be specified. (Units were previously always
6+
assumed to be Hz.)
7+
8+
.. plot::
9+
:include-source: true
10+
:alt: Time series and its power spectral density (psd), where the psd is correctly labeled with frequency units
11+
12+
# Sampling period in units of days
13+
dt = 1/24
14+
15+
# Create example signal: sinusoid with red noise
16+
np.random.seed(19680801) # Fixing random state for reproducibility.
17+
t = np.arange(0, 20, dt)
18+
nse = np.random.randn(len(t))
19+
r = np.exp(-t / 0.05)
20+
cnse = np.convolve(nse, r) * dt
21+
cnse = cnse[:len(t)]
22+
s = 0.1 * np.sin(2 * np.pi * t) + cnse
23+
24+
# Show signal and power spectral density
25+
fig, (ax0, ax1) = plt.subplots(2, 1, layout='constrained')
26+
ax0.plot(t,s)
27+
ax0.set(xlabel='Time (d)', ylabel='Signal')
28+
ax1.psd(s, NFFT=256, Fs=1 / dt, Funits='cpd')
29+
plt.show()

‎galleries/examples/statistics/psd_demo.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,9 @@
3333
ax0.set_xlabel('Time (s)')
3434
ax0.set_ylabel('Signal')
3535
ax1.psd(s, NFFT=512, Fs=1 / dt)
36+
# If dt had other units (e.g. days instead of seconds),
37+
# then the units of Fs (e.g. cycles per day or cps) can be specified
38+
# by the keyword Funits (e.g. Funits='cpd') in psd.
3639

3740
plt.show()
3841

‎lib/matplotlib/axes/_axes.py‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8113,7 +8113,8 @@ def ecdf(self, x, weights=None, *, complementary=False,
81138113
@_docstring.interpd
81148114
def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None,
81158115
window=None, noverlap=None, pad_to=None,
8116-
sides=None, scale_by_freq=None, return_line=None, **kwargs):
8116+
sides=None, scale_by_freq=None, return_line=None, Funits=None,
8117+
**kwargs):
81178118
r"""
81188119
Plot the power spectral density.
81198120
@@ -8147,6 +8148,12 @@ def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None,
81478148
return_line : bool, default: False
81488149
Whether to include the line object plotted in the returned values.
81498150
8151+
Funits : str, default: 'Hz'
8152+
Units for the sampling frequency *Fs*. It is used to label the
8153+
xaxis and yaxis.
8154+
8155+
.. versionadded:: 3.12
8156+
81508157
Returns
81518158
-------
81528159
Pxx : 1-D array
@@ -8194,19 +8201,21 @@ def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None,
81948201
"""
81958202
if Fc is None:
81968203
Fc = 0
8204+
if Funits is None:
8205+
Funits = 'Hz'
81978206

81988207
pxx, freqs = mlab.psd(x=x, NFFT=NFFT, Fs=Fs, detrend=detrend,
81998208
window=window, noverlap=noverlap, pad_to=pad_to,
82008209
sides=sides, scale_by_freq=scale_by_freq)
82018210
freqs += Fc
82028211

82038212
if scale_by_freq in (None, True):
8204-
psd_units = 'dB/Hz'
8213+
psd_units = 'dB/%s' % Funits
82058214
else:
82068215
psd_units = 'dB'
82078216

82088217
line = self.plot(freqs, 10 * np.log10(pxx), **kwargs)
8209-
self.set_xlabel('Frequency')
8218+
self.set_xlabel('Frequency (%s)' % Funits)
82108219
self.set_ylabel('Power Spectral Density (%s)' % psd_units)
82118220
self.grid(True)
82128221

‎lib/matplotlib/axes/_axes.pyi‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -653,6 +653,7 @@ class Axes(_AxesBase):
653653
sides: Literal["default", "onesided", "twosided"] | None = ...,
654654
scale_by_freq: bool | None = ...,
655655
return_line: bool | None = ...,
656+
Funits: str | None = ...,
656657
data: DataParamType = ...,
657658
**kwargs
658659
) -> tuple[np.ndarray, np.ndarray] | tuple[np.ndarray, np.ndarray, Line2D]: ...

‎lib/matplotlib/mlab.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -362,9 +362,9 @@ def _spectral_helper(x, y=None, NFFT=None, Fs=None, detrend_func=None,
362362

363363
result[slc] *= scaling_factor
364364

365-
# MATLAB divides by the sampling frequency so that density function
366-
# has units of dB/Hz and can be integrated by the plotted frequency
367-
# values. Perform the same scaling here.
365+
# Divide by the sampling frequency so that density function
366+
# has units of V**2/Hz, if x is measured in units of V and the sampling
367+
# frequency is measured in Hz.
368368
if scale_by_freq:
369369
result /= Fs
370370
# Scale the spectrum by the norm of the window to compensate for
@@ -470,10 +470,10 @@ def _single_spectrum_helper(
470470
`.detrend_mean`. 'linear' calls `.detrend_linear`.
471471
472472
scale_by_freq : bool, default: True
473-
Whether the resulting density values should be scaled by the scaling
474-
frequency, which gives density in units of 1/Hz. This allows for
475-
integration over the returned frequency values. The default is True for
476-
MATLAB compatibility.""")
473+
Whether the resulting density values should be divided by the sampling
474+
frequency, which gives density in units of 1/Hz, if the sampling rate
475+
is measured in Hz. This allows for integration over the returned
476+
frequency values. The default is True for MATLAB compatibility.""")
477477

478478

479479
@_docstring.interpd

‎lib/matplotlib/pyplot.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4070,6 +4070,7 @@ def psd(
40704070
sides: Literal["default", "onesided", "twosided"] | None = None,
40714071
scale_by_freq: bool | None = None,
40724072
return_line: bool | None = None,
4073+
Funits: str | None = None,
40734074
*,
40744075
data: DataParamType = None,
40754076
**kwargs,
@@ -4086,6 +4087,7 @@ def psd(
40864087
sides=sides,
40874088
scale_by_freq=scale_by_freq,
40884089
return_line=return_line,
4090+
Funits=Funits,
40894091
**({"data": data} if data is not None else {}),
40904092
**kwargs,
40914093
)

0 commit comments

Comments
 (0)