diff --git a/doc/api/next_api_changes/development/32211_CDH.rst b/doc/api/next_api_changes/development/32211_CDH.rst new file mode 100644 index 000000000000..c51e2e76c068 --- /dev/null +++ b/doc/api/next_api_changes/development/32211_CDH.rst @@ -0,0 +1,6 @@ +Power Spectral Density +~~~~~~~~~~~~~~~~~~~~~~ +An optional keyword *Funits* is added to `.Axes.psd` so +units that differ from the default 'Hz' can be specified. +This change is backward compatible (no change if the optional +keyword is omitted.) diff --git a/doc/release/next_whats_new/new_psd_feature.rst b/doc/release/next_whats_new/new_psd_feature.rst new file mode 100644 index 000000000000..4338354436c7 --- /dev/null +++ b/doc/release/next_whats_new/new_psd_feature.rst @@ -0,0 +1,30 @@ +Sampling frequency units can be specified for `.Axes.psd` +--------------------------------------------------------------- + +When creating a power spectral density (psd) plot, the units of the +sampling frequency can be specified. (Units were previously always +assumed to be Hz.) + +:: + + import matplotlib.pyplot as plt + import numpy as np + + # Sampling period in units of days + dt = 1/24 + + # Create example signal: sinusoid with red noise + np.random.seed(19680801) # Fixing random state for reproducibility + t = np.arange(0, 20, dt) + nse = np.random.randn(len(t)) + r = np.exp(-t / 0.05) + cnse = np.convolve(nse, r) * dt + cnse = cnse[:len(t)] + s = 0.1 * np.sin(2 * np.pi * t) + cnse + + # Show signal and power spectral density + fig, (ax0, ax1) = plt.subplots(2, 1, layout='constrained') + ax0.plot(t,s) + ax0.set(xlabel='Time (d)', ylabel='Signal') + ax1.psd(s, NFFT=256, Fs=1 / dt, Funits='cpd') + plt.show() diff --git a/galleries/examples/statistics/psd_demo.py b/galleries/examples/statistics/psd_demo.py index bf564df7542c..fd5c4646cb29 100644 --- a/galleries/examples/statistics/psd_demo.py +++ b/galleries/examples/statistics/psd_demo.py @@ -33,6 +33,9 @@ ax0.set_xlabel('Time (s)') ax0.set_ylabel('Signal') ax1.psd(s, NFFT=512, Fs=1 / dt) +# If dt had other units (e.g. days instead of seconds), +# then the units of Fs (e.g. cycles per day or cps) can be specified +# by the keyword Funits (e.g. Funits='cpd') in psd. plt.show() diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py index 09aef856dc48..2fdc21108094 100644 --- a/lib/matplotlib/axes/_axes.py +++ b/lib/matplotlib/axes/_axes.py @@ -8113,7 +8113,8 @@ def ecdf(self, x, weights=None, *, complementary=False, @_docstring.interpd def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None, window=None, noverlap=None, pad_to=None, - sides=None, scale_by_freq=None, return_line=None, **kwargs): + sides=None, scale_by_freq=None, return_line=None, Funits=None, + **kwargs): r""" Plot the power spectral density. @@ -8147,6 +8148,10 @@ def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None, return_line : bool, default: False Whether to include the line object plotted in the returned values. + Funits : str, default: 'Hz' + Units for the sampling frequency *Fs*. It is used to label the + xaxis and yaxis. + Returns ------- Pxx : 1-D array @@ -8194,6 +8199,8 @@ def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None, """ if Fc is None: Fc = 0 + if Funits is None: + Funits = 'Hz' pxx, freqs = mlab.psd(x=x, NFFT=NFFT, Fs=Fs, detrend=detrend, window=window, noverlap=noverlap, pad_to=pad_to, @@ -8201,12 +8208,12 @@ def psd(self, x, NFFT=None, Fs=None, Fc=None, detrend=None, freqs += Fc if scale_by_freq in (None, True): - psd_units = 'dB/Hz' + psd_units = 'dB/%s' % Funits else: psd_units = 'dB' line = self.plot(freqs, 10 * np.log10(pxx), **kwargs) - self.set_xlabel('Frequency') + self.set_xlabel('Frequency (%s)' % Funits) self.set_ylabel('Power Spectral Density (%s)' % psd_units) self.grid(True) diff --git a/lib/matplotlib/axes/_axes.pyi b/lib/matplotlib/axes/_axes.pyi index 227eba9be6d4..855ba30b07f8 100644 --- a/lib/matplotlib/axes/_axes.pyi +++ b/lib/matplotlib/axes/_axes.pyi @@ -654,6 +654,7 @@ class Axes(_AxesBase): scale_by_freq: bool | None = ..., return_line: bool | None = ..., data: DataParamType = ..., + Funits: str | None = ..., **kwargs ) -> tuple[np.ndarray, np.ndarray] | tuple[np.ndarray, np.ndarray, Line2D]: ... def csd( diff --git a/lib/matplotlib/mlab.py b/lib/matplotlib/mlab.py index a694308384c1..2063c9744c31 100644 --- a/lib/matplotlib/mlab.py +++ b/lib/matplotlib/mlab.py @@ -363,8 +363,8 @@ def _spectral_helper(x, y=None, NFFT=None, Fs=None, detrend_func=None, result[slc] *= scaling_factor # MATLAB divides by the sampling frequency so that density function - # has units of dB/Hz and can be integrated by the plotted frequency - # values. Perform the same scaling here. + # has units of V**2/Hz, if x is measured in units of V and the sampling + # frequency is measured in Hz. Perform the same scaling here. if scale_by_freq: result /= Fs # Scale the spectrum by the norm of the window to compensate for @@ -470,10 +470,10 @@ def _single_spectrum_helper( `.detrend_mean`. 'linear' calls `.detrend_linear`. scale_by_freq : bool, default: True - Whether the resulting density values should be scaled by the scaling - frequency, which gives density in units of 1/Hz. This allows for - integration over the returned frequency values. The default is True for - MATLAB compatibility.""") + Whether the resulting density values should be divided by the sampling + frequency, which gives density in units of 1/Hz, if the sampling rate + is measured in Hz. This allows for integration over the returned + frequency values. The default is True for MATLAB compatibility.""") @_docstring.interpd diff --git a/lib/matplotlib/pyplot.py b/lib/matplotlib/pyplot.py index 8315056c81a2..6671e07af64c 100644 --- a/lib/matplotlib/pyplot.py +++ b/lib/matplotlib/pyplot.py @@ -4070,6 +4070,7 @@ def psd( sides: Literal["default", "onesided", "twosided"] | None = None, scale_by_freq: bool | None = None, return_line: bool | None = None, + Funits: str | None = None, *, data: DataParamType = None, **kwargs, @@ -4086,6 +4087,7 @@ def psd( sides=sides, scale_by_freq=scale_by_freq, return_line=return_line, + Funits=Funits, **({"data": data} if data is not None else {}), **kwargs, )