Aluode/PerceptionLabPortable
0
1"""Tools for spectral analysis.
2"""
3import numpy as np
4import numpy.typing as npt
5from scipy import fft as sp_fft
6from . import _signaltools
7from .windows import get_window
8from ._arraytools import const_ext, even_ext, odd_ext, zero_ext
9import warnings
10from typing import Literal
11
12
13__all__ = ['periodogram', 'welch', 'lombscargle', 'csd', 'coherence',
14 'spectrogram', 'stft', 'istft', 'check_COLA', 'check_NOLA']
15
16
17def lombscargle(
18 x: npt.ArrayLike,
19 y: npt.ArrayLike,
20 freqs: npt.ArrayLike,
21 precenter: bool = False,
22 normalize: bool | Literal["power", "normalize", "amplitude"] = False,
23 *,
24 weights: npt.NDArray | None = None,
25 floating_mean: bool = False,
26) -> npt.NDArray:
27 """
28 Compute the generalized Lomb-Scargle periodogram.
29
30 The Lomb-Scargle periodogram was developed by Lomb [1]_ and further
31 extended by Scargle [2]_ to find, and test the significance of weak
32 periodic signals with uneven temporal sampling. The algorithm used
33 here is based on a weighted least-squares fit of the form
34 ``y(ω) = a*cos(ω*x) + b*sin(ω*x) + c``, where the fit is calculated for
35 each frequency independently. This algorithm was developed by Zechmeister
36 and Kürster which improves the Lomb-Scargle periodogram by enabling
37 the weighting of individual samples and calculating an unknown y offset
38 (also called a "floating-mean" model) [3]_. For more details, and practical
39 considerations, see the excellent reference on the Lomb-Scargle periodogram [4]_.
40
41 When *normalize* is False (or "power") (default) the computed periodogram
42 is unnormalized, it takes the value ``(A**2) * N/4`` for a harmonic
43 signal with amplitude A for sufficiently large N. Where N is the length of x or y.
44
45 When *normalize* is True (or "normalize") the computed periodogram is normalized
46 by the residuals of the data around a constant reference model (at zero).
47
48 When *normalize* is "amplitude" the computed periodogram is the complex
49 representation of the amplitude and phase.
50
51 Input arrays should be 1-D of a real floating data type, which are converted into
52 float64 arrays before processing.
53
54 Parameters
55 ----------
56 x : array_like
57 Sample times.
58 y : array_like
59 Measurement values. Values are assumed to have a baseline of ``y = 0``. If
60 there is a possibility of a y offset, it is recommended to set `floating_mean`
61 to True.
62 freqs : array_like
63 Angular frequencies (e.g., having unit rad/s=2π/s for `x` having unit s) for
64 output periodogram. Frequencies are normally >= 0, as any peak at ``-freq`` will
65 also exist at ``+freq``.
66 precenter : bool, optional
67 Pre-center measurement values by subtracting the mean, if True. This is
68 a legacy parameter and unnecessary if `floating_mean` is True.
69 normalize : bool | str, optional
70 Compute normalized or complex (amplitude + phase) periodogram.
71 Valid options are: ``False``/``"power"``, ``True``/``"normalize"``, or
72 ``"amplitude"``.
73 weights : array_like, optional
74 Weights for each sample. Weights must be nonnegative.
75 floating_mean : bool, optional
76 Determines a y offset for each frequency independently, if True.
77 Else the y offset is assumed to be `0`.
78
79 Returns
80 -------
81 pgram : array_like
82 Lomb-Scargle periodogram.
83
84 Raises
85 ------
86 ValueError
87 If any of the input arrays x, y, freqs, or weights are not 1D, or if any are
88 zero length. Or, if the input arrays x, y, and weights do not have the same
89 shape as each other.
90 ValueError
91 If any weight is < 0, or the sum of the weights is <= 0.
92 ValueError
93 If the normalize parameter is not one of the allowed options.
94
95 See Also
96 --------
97 periodogram: Power spectral density using a periodogram
98 welch: Power spectral density by Welch's method
99 csd: Cross spectral density by Welch's method
100
101 Notes
102 -----
103 The algorithm used will not automatically account for any unknown y offset, unless
104 floating_mean is True. Therefore, for most use cases, if there is a possibility of
105 a y offset, it is recommended to set floating_mean to True. If precenter is True,
106 it performs the operation ``y -= y.mean()``. However, precenter is a legacy
107 parameter, and unnecessary when floating_mean is True. Furthermore, the mean
108 removed by precenter does not account for sample weights, nor will it correct for
109 any bias due to consistently missing observations at peaks and/or troughs. When the
110 normalize parameter is "amplitude", for any frequency in freqs that is below
111 ``(2*pi)/(x.max() - x.min())``, the predicted amplitude will tend towards infinity.
112 The concept of a "Nyquist frequency" limit (see Nyquist-Shannon sampling theorem)
113 is not generally applicable to unevenly sampled data. Therefore, with unevenly
114 sampled data, valid frequencies in freqs can often be much higher than expected.
115
116 References
117 ----------
118 .. [1] N.R. Lomb "Least-squares frequency analysis of unequally spaced
119 data", Astrophysics and Space Science, vol 39, pp. 447-462, 1976
120 :doi:`10.1007/bf00648343`
121
122 .. [2] J.D. Scargle "Studies in astronomical time series analysis. II -
123 Statistical aspects of spectral analysis of unevenly spaced data",
124 The Astrophysical Journal, vol 263, pp. 835-853, 1982
125 :doi:`10.1086/160554`
126
127 .. [3] M. Zechmeister and M. Kürster, "The generalised Lomb-Scargle periodogram.
128 A new formalism for the floating-mean and Keplerian periodograms,"
129 Astronomy and Astrophysics, vol. 496, pp. 577-584, 2009
130 :doi:`10.1051/0004-6361:200811296`
131
132 .. [4] J.T. VanderPlas, "Understanding the Lomb-Scargle Periodogram,"
133 The Astrophysical Journal Supplement Series, vol. 236, no. 1, p. 16,
134 May 2018
135 :doi:`10.3847/1538-4365/aab766`
136
137
138 Examples
139 --------
140 >>> import numpy as np
141 >>> rng = np.random.default_rng()
142
143 First define some input parameters for the signal:
144
145 >>> A = 2. # amplitude
146 >>> c = 2. # offset
147 >>> w0 = 1. # rad/sec
148 >>> nin = 150
149 >>> nout = 1002
150
151 Randomly generate sample times:
152
153 >>> x = rng.uniform(0, 10*np.pi, nin)
154
155 Plot a sine wave for the selected times:
156
157 >>> y = A * np.cos(w0*x) + c
158
159 Define the array of frequencies for which to compute the periodogram:
160
161 >>> w = np.linspace(0.25, 10, nout)
162
163 Calculate Lomb-Scargle periodogram for each of the normalize options:
164
165 >>> from scipy.signal import lombscargle
166 >>> pgram_power = lombscargle(x, y, w, normalize=False)
167 >>> pgram_norm = lombscargle(x, y, w, normalize=True)
168 >>> pgram_amp = lombscargle(x, y, w, normalize='amplitude')
169 ...
170 >>> pgram_power_f = lombscargle(x, y, w, normalize=False, floating_mean=True)
171 >>> pgram_norm_f = lombscargle(x, y, w, normalize=True, floating_mean=True)
172 >>> pgram_amp_f = lombscargle(x, y, w, normalize='amplitude', floating_mean=True)
173
174 Now make a plot of the input data:
175
176 >>> import matplotlib.pyplot as plt
177 >>> fig, (ax_t, ax_p, ax_n, ax_a) = plt.subplots(4, 1, figsize=(5, 6))
178 >>> ax_t.plot(x, y, 'b+')
179 >>> ax_t.set_xlabel('Time [s]')
180 >>> ax_t.set_ylabel('Amplitude')
181
182 Then plot the periodogram for each of the normalize options, as well as with and
183 without floating_mean=True:
184
185 >>> ax_p.plot(w, pgram_power, label='default')
186 >>> ax_p.plot(w, pgram_power_f, label='floating_mean=True')
187 >>> ax_p.set_xlabel('Angular frequency [rad/s]')
188 >>> ax_p.set_ylabel('Power')
189 >>> ax_p.legend(prop={'size': 7})
190 ...
191 >>> ax_n.plot(w, pgram_norm, label='default')
192 >>> ax_n.plot(w, pgram_norm_f, label='floating_mean=True')
193 >>> ax_n.set_xlabel('Angular frequency [rad/s]')
194 >>> ax_n.set_ylabel('Normalized')
195 >>> ax_n.legend(prop={'size': 7})
196 ...
197 >>> ax_a.plot(w, np.abs(pgram_amp), label='default')
198 >>> ax_a.plot(w, np.abs(pgram_amp_f), label='floating_mean=True')
199 >>> ax_a.set_xlabel('Angular frequency [rad/s]')
200 >>> ax_a.set_ylabel('Amplitude')
201 >>> ax_a.legend(prop={'size': 7})
202 ...
203 >>> plt.tight_layout()
204 >>> plt.show()
205
206 """
207
208 # if no weights are provided, assume all data points are equally important
209 if weights is None:
210 weights = np.ones_like(y, dtype=np.float64)
211 else:
212 # if provided, make sure weights is an array and cast to float64
213 weights = np.asarray(weights, dtype=np.float64)
214
215 # make sure other inputs are arrays and cast to float64
216 # done before validation, in case they were not arrays
217 x = np.asarray(x, dtype=np.float64)
218 y = np.asarray(y, dtype=np.float64)
219 freqs = np.asarray(freqs, dtype=np.float64)
220
221 # validate input shapes
222 if not (x.ndim == 1 and x.size > 0 and x.shape == y.shape == weights.shape):
223 raise ValueError("Parameters x, y, weights must be 1-D arrays of "
224 "equal non-zero length!")
225 if not (freqs.ndim == 1 and freqs.size > 0):
226 raise ValueError("Parameter freqs must be a 1-D array of non-zero length!")
227
228 # validate weights
229 if not (np.all(weights >= 0) and np.sum(weights) > 0):
230 raise ValueError("Parameter weights must have only non-negative entries "
231 "which sum to a positive value!")
232
233 # validate normalize parameter
234 if isinstance(normalize, bool):
235 # if bool, convert to str literal
236 normalize = "normalize" if normalize else "power"
237
238 if normalize not in ["power", "normalize", "amplitude"]:
239 raise ValueError(
240 "Normalize must be: False (or 'power'), True (or 'normalize'), "
241 "or 'amplitude'."
242 )
243
244 # weight vector must sum to 1
245 weights *= 1.0 / weights.sum()
246
247 # if requested, perform precenter
248 if precenter:
249 y -= y.mean()
250
251 # transform arrays
252 # row vector
253 freqs = freqs.reshape(1, -1)
254 # column vectors
255 x = x.reshape(-1, 1)
256 y = y.reshape(-1, 1)
257 weights = weights.reshape(-1, 1)
258
259 # store frequent intermediates
260 weights_y = weights * y
261 freqst = freqs * x
262 coswt = np.cos(freqst)
263 sinwt = np.sin(freqst)
264
265 Y = np.dot(weights.T, y) # Eq. 7
266 CC = np.dot(weights.T, coswt * coswt) # Eq. 13
267 SS = 1.0 - CC # trig identity: S^2 = 1 - C^2 Eq.14
268 CS = np.dot(weights.T, coswt * sinwt) # Eq. 15
269
270 if floating_mean:
271 C = np.dot(weights.T, coswt) # Eq. 8
272 S = np.dot(weights.T, sinwt) # Eq. 9
273 CC -= C * C # Eq. 13
274 SS -= S * S # Eq. 14
275 CS -= C * S # Eq. 15
276
277 # calculate tau (phase offset to eliminate CS variable)
278 tau = 0.5 * np.arctan2(2.0 * CS, CC - SS) # Eq. 19
279 freqst_tau = freqst - tau
280
281 # coswt and sinwt are now offset by tau, which eliminates CS
282 coswt_tau = np.cos(freqst_tau)
283 sinwt_tau = np.sin(freqst_tau)
284
285 YC = np.dot(weights_y.T, coswt_tau) # Eq. 11
286 YS = np.dot(weights_y.T, sinwt_tau) # Eq. 12
287 CC = np.dot(weights.T, coswt_tau * coswt_tau) # Eq. 13, CC range is [0.5, 1.0]
288 SS = 1.0 - CC # trig identity: S^2 = 1 - C^2 Eq. 14, SS range is [0.0, 0.5]
289
290 if floating_mean:
291 C = np.dot(weights.T, coswt_tau) # Eq. 8
292 S = np.dot(weights.T, sinwt_tau) # Eq. 9
293 YC -= Y * C # Eq. 11
294 YS -= Y * S # Eq. 12
295 CC -= C * C # Eq. 13, CC range is now [0.0, 1.0]
296 SS -= S * S # Eq. 14, SS range is now [0.0, 0.5]
297
298 # to prevent division by zero errors with a and b, as well as correcting for
299 # numerical precision errors that lead to CC or SS being approximately -0.0,
300 # make sure CC and SS are both > 0
301 epsneg = np.finfo(dtype=y.dtype).epsneg
302 CC[CC < epsneg] = epsneg
303 SS[SS < epsneg] = epsneg
304
305 # calculate a and b
306 # where: y(w) = a*cos(w) + b*sin(w) + c
307 a = YC / CC # Eq. A.4 and 6, eliminating CS
308 b = YS / SS # Eq. A.4 and 6, eliminating CS
309 # c = Y - a * C - b * S
310
311 # store final value as power in A^2 (i.e., (y units)^2)
312 pgram = 2.0 * (a * YC + b * YS)
313
314 # squeeze back to a vector
315 pgram = np.squeeze(pgram)
316
317 if normalize == "power": # (default)
318 # return the legacy power units ((A**2) * N/4)
319
320 pgram *= float(x.shape[0]) / 4.0
321
322 elif normalize == "normalize":
323 # return the normalized power (power at current frequency wrt the entire signal)
324 # range will be [0, 1]
325
326 YY = np.dot(weights_y.T, y) # Eq. 10
327 if floating_mean:
328 YY -= Y * Y # Eq. 10
329
330 pgram *= 0.5 / np.squeeze(YY) # Eq. 20
331
332 else: # normalize == "amplitude":
333 # return the complex representation of the best-fit amplitude and phase
334
335 # squeeze back to vectors
336 a = np.squeeze(a)
337 b = np.squeeze(b)
338 tau = np.squeeze(tau)
339
340 # calculate the complex representation, and correct for tau rotation
341 pgram = (a + 1j * b) * np.exp(1j * tau)
342
343 return pgram
344
345
346def periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant',
347 return_onesided=True, scaling='density', axis=-1):
348 """
349 Estimate power spectral density using a periodogram.
350
351 Parameters
352 ----------
353 x : array_like
354 Time series of measurement values
355 fs : float, optional
356 Sampling frequency of the `x` time series. Defaults to 1.0.
357 window : str or tuple or array_like, optional
358 Desired window to use. If `window` is a string or tuple, it is
359 passed to `get_window` to generate the window values, which are
360 DFT-even by default. See `get_window` for a list of windows and
361 required parameters. If `window` is array_like it will be used
362 directly as the window and its length must be equal to the length
363 of the axis over which the periodogram is computed. Defaults
364 to 'boxcar'.
365 nfft : int, optional
366 Length of the FFT used. If `None` the length of `x` will be
367 used.
368 detrend : str or function or `False`, optional
369 Specifies how to detrend each segment. If `detrend` is a
370 string, it is passed as the `type` argument to the `detrend`
371 function. If it is a function, it takes a segment and returns a
372 detrended segment. If `detrend` is `False`, no detrending is
373 done. Defaults to 'constant'.
374 return_onesided : bool, optional
375 If `True`, return a one-sided spectrum for real data. If
376 `False` return a two-sided spectrum. Defaults to `True`, but for
377 complex data, a two-sided spectrum is always returned.
378 scaling : { 'density', 'spectrum' }, optional
379 Selects between computing the power spectral density ('density')
380 where `Pxx` has units of V**2/Hz and computing the squared magnitude
381 spectrum ('spectrum') where `Pxx` has units of V**2, if `x`
382 is measured in V and `fs` is measured in Hz. Defaults to
383 'density'
384 axis : int, optional
385 Axis along which the periodogram is computed; the default is
386 over the last axis (i.e. ``axis=-1``).
387
388 Returns
389 -------
390 f : ndarray
391 Array of sample frequencies.
392 Pxx : ndarray
393 Power spectral density or power spectrum of `x`.
394
395 See Also
396 --------
397 welch: Estimate power spectral density using Welch's method
398 lombscargle: Lomb-Scargle periodogram for unevenly sampled data
399
400 Notes
401 -----
402 Consult the :ref:`tutorial_SpectralAnalysis` section of the :ref:`user_guide`
403 for a discussion of the scalings of the power spectral density and
404 the magnitude (squared) spectrum.
405
406 .. versionadded:: 0.12.0
407
408 Examples
409 --------
410 >>> import numpy as np
411 >>> from scipy import signal
412 >>> import matplotlib.pyplot as plt
413 >>> rng = np.random.default_rng()
414
415 Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
416 0.001 V**2/Hz of white noise sampled at 10 kHz.
417
418 >>> fs = 10e3
419 >>> N = 1e5
420 >>> amp = 2*np.sqrt(2)
421 >>> freq = 1234.0
422 >>> noise_power = 0.001 * fs / 2
423 >>> time = np.arange(N) / fs
424 >>> x = amp*np.sin(2*np.pi*freq*time)
425 >>> x += rng.normal(scale=np.sqrt(noise_power), size=time.shape)
426
427 Compute and plot the power spectral density.
428
429 >>> f, Pxx_den = signal.periodogram(x, fs)
430 >>> plt.semilogy(f, Pxx_den)
431 >>> plt.ylim([1e-7, 1e2])
432 >>> plt.xlabel('frequency [Hz]')
433 >>> plt.ylabel('PSD [V**2/Hz]')
434 >>> plt.show()
435
436 If we average the last half of the spectral density, to exclude the
437 peak, we can recover the noise power on the signal.
438
439 >>> np.mean(Pxx_den[25000:])
440 0.000985320699252543
441
442 Now compute and plot the power spectrum.
443
444 >>> f, Pxx_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')
445 >>> plt.figure()
446 >>> plt.semilogy(f, np.sqrt(Pxx_spec))
447 >>> plt.ylim([1e-4, 1e1])
448 >>> plt.xlabel('frequency [Hz]')
449 >>> plt.ylabel('Linear spectrum [V RMS]')
450 >>> plt.show()
451
452 The peak height in the power spectrum is an estimate of the RMS
453 amplitude.
454
455 >>> np.sqrt(Pxx_spec.max())
456 2.0077340678640727
457
458 """
459 x = np.asarray(x)
460
461 if x.size == 0:
462 return np.empty(x.shape), np.empty(x.shape)
463
464 if window is None:
465 window = 'boxcar'
466
467 if nfft is None:
468 nperseg = x.shape[axis]
469 elif nfft == x.shape[axis]:
470 nperseg = nfft
471 elif nfft > x.shape[axis]:
472 nperseg = x.shape[axis]
473 elif nfft < x.shape[axis]:
474 s = [np.s_[:]]*len(x.shape)
475 s[axis] = np.s_[:nfft]
476 x = x[tuple(s)]
477 nperseg = nfft
478 nfft = None
479
480 if hasattr(window, 'size'):
481 if window.size != nperseg:
482 raise ValueError('the size of the window must be the same size '
483 'of the input on the specified axis')
484
485 return welch(x, fs=fs, window=window, nperseg=nperseg, noverlap=0,
486 nfft=nfft, detrend=detrend, return_onesided=return_onesided,
487 scaling=scaling, axis=axis)
488
489
490def welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None,
491 detrend='constant', return_onesided=True, scaling='density',
492 axis=-1, average='mean'):
493 r"""
494 Estimate power spectral density using Welch's method.
495
496 Welch's method [1]_ computes an estimate of the power spectral
497 density by dividing the data into overlapping segments, computing a
498 modified periodogram for each segment and averaging the
499 periodograms.
500
501 Parameters
502 ----------
503 x : array_like
504 Time series of measurement values
505 fs : float, optional
506 Sampling frequency of the `x` time series. Defaults to 1.0.
507 window : str or tuple or array_like, optional
508 Desired window to use. If `window` is a string or tuple, it is
509 passed to `get_window` to generate the window values, which are
510 DFT-even by default. See `get_window` for a list of windows and
511 required parameters. If `window` is array_like it will be used
512 directly as the window and its length must be nperseg. Defaults
513 to a Hann window.
514 nperseg : int, optional
515 Length of each segment. Defaults to None, but if window is str or
516 tuple, is set to 256, and if window is array_like, is set to the
517 length of the window.
518 noverlap : int, optional
519 Number of points to overlap between segments. If `None`,
520 ``noverlap = nperseg // 2``. Defaults to `None`.
521 nfft : int, optional
522 Length of the FFT used, if a zero padded FFT is desired. If
523 `None`, the FFT length is `nperseg`. Defaults to `None`.
524 detrend : str or function or `False`, optional
525 Specifies how to detrend each segment. If `detrend` is a
526 string, it is passed as the `type` argument to the `detrend`
527 function. If it is a function, it takes a segment and returns a
528 detrended segment. If `detrend` is `False`, no detrending is
529 done. Defaults to 'constant'.
530 return_onesided : bool, optional
531 If `True`, return a one-sided spectrum for real data. If
532 `False` return a two-sided spectrum. Defaults to `True`, but for
533 complex data, a two-sided spectrum is always returned.
534 scaling : { 'density', 'spectrum' }, optional
535 Selects between computing the power spectral density ('density')
536 where `Pxx` has units of V**2/Hz and computing the squared magnitude
537 spectrum ('spectrum') where `Pxx` has units of V**2, if `x`
538 is measured in V and `fs` is measured in Hz. Defaults to
539 'density'
540 axis : int, optional
541 Axis along which the periodogram is computed; the default is
542 over the last axis (i.e. ``axis=-1``).
543 average : { 'mean', 'median' }, optional
544 Method to use when averaging periodograms. Defaults to 'mean'.
545
546 .. versionadded:: 1.2.0
547
548 Returns
549 -------
550 f : ndarray
551 Array of sample frequencies.
552 Pxx : ndarray
553 Power spectral density or power spectrum of x.
554
555 See Also
556 --------
557 periodogram: Simple, optionally modified periodogram
558 lombscargle: Lomb-Scargle periodogram for unevenly sampled data
559
560 Notes
561 -----
562 An appropriate amount of overlap will depend on the choice of window
563 and on your requirements. For the default Hann window an overlap of
564 50% is a reasonable trade off between accurately estimating the
565 signal power, while not over counting any of the data. Narrower
566 windows may require a larger overlap.
567
568 If `noverlap` is 0, this method is equivalent to Bartlett's method
569 [2]_.
570
571 Consult the :ref:`tutorial_SpectralAnalysis` section of the :ref:`user_guide`
572 for a discussion of the scalings of the power spectral density and
573 the (squared) magnitude spectrum.
574
575 .. versionadded:: 0.12.0
576
577 References
578 ----------
579 .. [1] P. Welch, "The use of the fast Fourier transform for the
580 estimation of power spectra: A method based on time averaging
581 over short, modified periodograms", IEEE Trans. Audio
582 Electroacoust. vol. 15, pp. 70-73, 1967.
583 .. [2] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
584 Biometrika, vol. 37, pp. 1-16, 1950.
585
586 Examples
587 --------
588 >>> import numpy as np
589 >>> from scipy import signal
590 >>> import matplotlib.pyplot as plt
591 >>> rng = np.random.default_rng()
592
593 Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
594 0.001 V**2/Hz of white noise sampled at 10 kHz.
595
596 >>> fs = 10e3
597 >>> N = 1e5
598 >>> amp = 2*np.sqrt(2)
599 >>> freq = 1234.0
600 >>> noise_power = 0.001 * fs / 2
601 >>> time = np.arange(N) / fs
602 >>> x = amp*np.sin(2*np.pi*freq*time)
603 >>> x += rng.normal(scale=np.sqrt(noise_power), size=time.shape)
604
605 Compute and plot the power spectral density.
606
607 >>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
608 >>> plt.semilogy(f, Pxx_den)
609 >>> plt.ylim([0.5e-3, 1])
610 >>> plt.xlabel('frequency [Hz]')
611 >>> plt.ylabel('PSD [V**2/Hz]')
612 >>> plt.show()
613
614 If we average the last half of the spectral density, to exclude the
615 peak, we can recover the noise power on the signal.
616
617 >>> np.mean(Pxx_den[256:])
618 0.0009924865443739191
619
620 Now compute and plot the power spectrum.
621
622 >>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
623 >>> plt.figure()
624 >>> plt.semilogy(f, np.sqrt(Pxx_spec))
625 >>> plt.xlabel('frequency [Hz]')
626 >>> plt.ylabel('Linear spectrum [V RMS]')
627 >>> plt.show()
628
629 The peak height in the power spectrum is an estimate of the RMS
630 amplitude.
631
632 >>> np.sqrt(Pxx_spec.max())
633 2.0077340678640727
634
635 If we now introduce a discontinuity in the signal, by increasing the
636 amplitude of a small portion of the signal by 50, we can see the
637 corruption of the mean average power spectral density, but using a
638 median average better estimates the normal behaviour.
639
640 >>> x[int(N//2):int(N//2)+10] *= 50.
641 >>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
642 >>> f_med, Pxx_den_med = signal.welch(x, fs, nperseg=1024, average='median')
643 >>> plt.semilogy(f, Pxx_den, label='mean')
644 >>> plt.semilogy(f_med, Pxx_den_med, label='median')
645 >>> plt.ylim([0.5e-3, 1])
646 >>> plt.xlabel('frequency [Hz]')
647 >>> plt.ylabel('PSD [V**2/Hz]')
648 >>> plt.legend()
649 >>> plt.show()
650
651 """
652 freqs, Pxx = csd(x, x, fs=fs, window=window, nperseg=nperseg,
653 noverlap=noverlap, nfft=nfft, detrend=detrend,
654 return_onesided=return_onesided, scaling=scaling,
655 axis=axis, average=average)
656
657 return freqs, Pxx.real
658
659
660def csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None,
661 detrend='constant', return_onesided=True, scaling='density',
662 axis=-1, average='mean'):
663 r"""
664 Estimate the cross power spectral density, Pxy, using Welch's method.
665
666 Parameters
667 ----------
668 x : array_like
669 Time series of measurement values
670 y : array_like
671 Time series of measurement values
672 fs : float, optional
673 Sampling frequency of the `x` and `y` time series. Defaults
674 to 1.0.
675 window : str or tuple or array_like, optional
676 Desired window to use. If `window` is a string or tuple, it is
677 passed to `get_window` to generate the window values, which are
678 DFT-even by default. See `get_window` for a list of windows and
679 required parameters. If `window` is array_like it will be used
680 directly as the window and its length must be nperseg. Defaults
681 to a Hann window.
682 nperseg : int, optional
683 Length of each segment. Defaults to None, but if window is str or
684 tuple, is set to 256, and if window is array_like, is set to the
685 length of the window.
686 noverlap: int, optional
687 Number of points to overlap between segments. If `None`,
688 ``noverlap = nperseg // 2``. Defaults to `None`.
689 nfft : int, optional
690 Length of the FFT used, if a zero padded FFT is desired. If
691 `None`, the FFT length is `nperseg`. Defaults to `None`.
692 detrend : str or function or `False`, optional
693 Specifies how to detrend each segment. If `detrend` is a
694 string, it is passed as the `type` argument to the `detrend`
695 function. If it is a function, it takes a segment and returns a
696 detrended segment. If `detrend` is `False`, no detrending is
697 done. Defaults to 'constant'.
698 return_onesided : bool, optional
699 If `True`, return a one-sided spectrum for real data. If
700 `False` return a two-sided spectrum. Defaults to `True`, but for
701 complex data, a two-sided spectrum is always returned.
702 scaling : { 'density', 'spectrum' }, optional
703 Selects between computing the cross spectral density ('density')
704 where `Pxy` has units of V**2/Hz and computing the cross spectrum
705 ('spectrum') where `Pxy` has units of V**2, if `x` and `y` are
706 measured in V and `fs` is measured in Hz. Defaults to 'density'
707 axis : int, optional
708 Axis along which the CSD is computed for both inputs; the
709 default is over the last axis (i.e. ``axis=-1``).
710 average : { 'mean', 'median' }, optional
711 Method to use when averaging periodograms. If the spectrum is
712 complex, the average is computed separately for the real and
713 imaginary parts. Defaults to 'mean'.
714
715 .. versionadded:: 1.2.0
716
717 Returns
718 -------
719 f : ndarray
720 Array of sample frequencies.
721 Pxy : ndarray
722 Cross spectral density or cross power spectrum of x,y.
723
724 See Also
725 --------
726 periodogram: Simple, optionally modified periodogram
727 lombscargle: Lomb-Scargle periodogram for unevenly sampled data
728 welch: Power spectral density by Welch's method. [Equivalent to
729 csd(x,x)]
730 coherence: Magnitude squared coherence by Welch's method.
731
732 Notes
733 -----
734 By convention, Pxy is computed with the conjugate FFT of X
735 multiplied by the FFT of Y.
736
737 If the input series differ in length, the shorter series will be
738 zero-padded to match.
739
740 An appropriate amount of overlap will depend on the choice of window
741 and on your requirements. For the default Hann window an overlap of
742 50% is a reasonable trade off between accurately estimating the
743 signal power, while not over counting any of the data. Narrower
744 windows may require a larger overlap.
745
746 Consult the :ref:`tutorial_SpectralAnalysis` section of the :ref:`user_guide`
747 for a discussion of the scalings of a spectral density and an (amplitude) spectrum.
748
749 .. versionadded:: 0.16.0
750
751 References
752 ----------
753 .. [1] P. Welch, "The use of the fast Fourier transform for the
754 estimation of power spectra: A method based on time averaging
755 over short, modified periodograms", IEEE Trans. Audio
756 Electroacoust. vol. 15, pp. 70-73, 1967.
757 .. [2] Rabiner, Lawrence R., and B. Gold. "Theory and Application of
758 Digital Signal Processing" Prentice-Hall, pp. 414-419, 1975
759
760 Examples
761 --------
762 >>> import numpy as np
763 >>> from scipy import signal
764 >>> import matplotlib.pyplot as plt
765 >>> rng = np.random.default_rng()
766
767 Generate two test signals with some common features.
768
769 >>> fs = 10e3
770 >>> N = 1e5
771 >>> amp = 20
772 >>> freq = 1234.0
773 >>> noise_power = 0.001 * fs / 2
774 >>> time = np.arange(N) / fs
775 >>> b, a = signal.butter(2, 0.25, 'low')
776 >>> x = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
777 >>> y = signal.lfilter(b, a, x)
778 >>> x += amp*np.sin(2*np.pi*freq*time)
779 >>> y += rng.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)
780
781 Compute and plot the magnitude of the cross spectral density.
782
783 >>> f, Pxy = signal.csd(x, y, fs, nperseg=1024)
784 >>> plt.semilogy(f, np.abs(Pxy))
785 >>> plt.xlabel('frequency [Hz]')
786 >>> plt.ylabel('CSD [V**2/Hz]')
787 >>> plt.show()
788
789 """
790 freqs, _, Pxy = _spectral_helper(x, y, fs, window, nperseg, noverlap,
791 nfft, detrend, return_onesided, scaling,
792 axis, mode='psd')
793
794 # Average over windows.
795 if len(Pxy.shape) >= 2 and Pxy.size > 0:
796 if Pxy.shape[-1] > 1:
797 if average == 'median':
798 # np.median must be passed real arrays for the desired result
799 bias = _median_bias(Pxy.shape[-1])
800 if np.iscomplexobj(Pxy):
801 Pxy = (np.median(np.real(Pxy), axis=-1)
802 + 1j * np.median(np.imag(Pxy), axis=-1))
803 else:
804 Pxy = np.median(Pxy, axis=-1)
805 Pxy /= bias
806 elif average == 'mean':
807 Pxy = Pxy.mean(axis=-1)
808 else:
809 raise ValueError(f'average must be "median" or "mean", got {average}')
810 else:
811 Pxy = np.reshape(Pxy, Pxy.shape[:-1])
812
813 return freqs, Pxy
814
815
816def spectrogram(x, fs=1.0, window=('tukey', .25), nperseg=None, noverlap=None,
817 nfft=None, detrend='constant', return_onesided=True,
818 scaling='density', axis=-1, mode='psd'):
819 """Compute a spectrogram with consecutive Fourier transforms (legacy function).
820
821 Spectrograms can be used as a way of visualizing the change of a
822 nonstationary signal's frequency content over time.
823
824 .. legacy:: function
825
826 :class:`ShortTimeFFT` is a newer STFT / ISTFT implementation with more
827 features also including a :meth:`~ShortTimeFFT.spectrogram` method.
828 A :ref:`comparison <tutorial_stft_legacy_stft>` between the
829 implementations can be found in the :ref:`tutorial_stft` section of
830 the :ref:`user_guide`.
831
832 Parameters
833 ----------
834 x : array_like
835 Time series of measurement values
836 fs : float, optional
837 Sampling frequency of the `x` time series. Defaults to 1.0.
838 window : str or tuple or array_like, optional
839 Desired window to use. If `window` is a string or tuple, it is
840 passed to `get_window` to generate the window values, which are
841 DFT-even by default. See `get_window` for a list of windows and
842 required parameters. If `window` is array_like it will be used
843 directly as the window and its length must be nperseg.
844 Defaults to a Tukey window with shape parameter of 0.25.
845 nperseg : int, optional
846 Length of each segment. Defaults to None, but if window is str or
847 tuple, is set to 256, and if window is array_like, is set to the
848 length of the window.
849 noverlap : int, optional
850 Number of points to overlap between segments. If `None`,
851 ``noverlap = nperseg // 8``. Defaults to `None`.
852 nfft : int, optional
853 Length of the FFT used, if a zero padded FFT is desired. If
854 `None`, the FFT length is `nperseg`. Defaults to `None`.
855 detrend : str or function or `False`, optional
856 Specifies how to detrend each segment. If `detrend` is a
857 string, it is passed as the `type` argument to the `detrend`
858 function. If it is a function, it takes a segment and returns a
859 detrended segment. If `detrend` is `False`, no detrending is
860 done. Defaults to 'constant'.
861 return_onesided : bool, optional
862 If `True`, return a one-sided spectrum for real data. If
863 `False` return a two-sided spectrum. Defaults to `True`, but for
864 complex data, a two-sided spectrum is always returned.
865 scaling : { 'density', 'spectrum' }, optional
866 Selects between computing the power spectral density ('density')
867 where `Sxx` has units of V**2/Hz and computing the power
868 spectrum ('spectrum') where `Sxx` has units of V**2, if `x`
869 is measured in V and `fs` is measured in Hz. Defaults to
870 'density'.
871 axis : int, optional
872 Axis along which the spectrogram is computed; the default is over
873 the last axis (i.e. ``axis=-1``).
874 mode : str, optional
875 Defines what kind of return values are expected. Options are
876 ['psd', 'complex', 'magnitude', 'angle', 'phase']. 'complex' is
877 equivalent to the output of `stft` with no padding or boundary
878 extension. 'magnitude' returns the absolute magnitude of the
879 STFT. 'angle' and 'phase' return the complex angle of the STFT,
880 with and without unwrapping, respectively.
881
882 Returns
883 -------
884 f : ndarray
885 Array of sample frequencies.
886 t : ndarray
887 Array of segment times.
888 Sxx : ndarray
889 Spectrogram of x. By default, the last axis of Sxx corresponds
890 to the segment times.
891
892 See Also
893 --------
894 periodogram: Simple, optionally modified periodogram
895 lombscargle: Lomb-Scargle periodogram for unevenly sampled data
896 welch: Power spectral density by Welch's method.
897 csd: Cross spectral density by Welch's method.
898 ShortTimeFFT: Newer STFT/ISTFT implementation providing more features,
899 which also includes a :meth:`~ShortTimeFFT.spectrogram`
900 method.
901
902 Notes
903 -----
904 An appropriate amount of overlap will depend on the choice of window
905 and on your requirements. In contrast to welch's method, where the
906 entire data stream is averaged over, one may wish to use a smaller
907 overlap (or perhaps none at all) when computing a spectrogram, to
908 maintain some statistical independence between individual segments.
909 It is for this reason that the default window is a Tukey window with
910 1/8th of a window's length overlap at each end.
911
912
913 .. versionadded:: 0.16.0
914
915 References
916 ----------
917 .. [1] Oppenheim, Alan V., Ronald W. Schafer, John R. Buck
918 "Discrete-Time Signal Processing", Prentice Hall, 1999.
919
920 Examples
921 --------
922 >>> import numpy as np
923 >>> from scipy import signal
924 >>> from scipy.fft import fftshift
925 >>> import matplotlib.pyplot as plt
926 >>> rng = np.random.default_rng()
927
928 Generate a test signal, a 2 Vrms sine wave whose frequency is slowly
929 modulated around 3kHz, corrupted by white noise of exponentially
930 decreasing magnitude sampled at 10 kHz.
931
932 >>> fs = 10e3
933 >>> N = 1e5
934 >>> amp = 2 * np.sqrt(2)
935 >>> noise_power = 0.01 * fs / 2
936 >>> time = np.arange(N) / float(fs)
937 >>> mod = 500*np.cos(2*np.pi*0.25*time)
938 >>> carrier = amp * np.sin(2*np.pi*3e3*time + mod)
939 >>> noise = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
940 >>> noise *= np.exp(-time/5)
941 >>> x = carrier + noise
942
943 Compute and plot the spectrogram.
944
945 >>> f, t, Sxx = signal.spectrogram(x, fs)
946 >>> plt.pcolormesh(t, f, Sxx, shading='gouraud')
947 >>> plt.ylabel('Frequency [Hz]')
948 >>> plt.xlabel('Time [sec]')
949 >>> plt.show()
950
951 Note, if using output that is not one sided, then use the following:
952
953 >>> f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False)
954 >>> plt.pcolormesh(t, fftshift(f), fftshift(Sxx, axes=0), shading='gouraud')
955 >>> plt.ylabel('Frequency [Hz]')
956 >>> plt.xlabel('Time [sec]')
957 >>> plt.show()
958
959 """
960 modelist = ['psd', 'complex', 'magnitude', 'angle', 'phase']
961 if mode not in modelist:
962 raise ValueError(f'unknown value for mode {mode}, must be one of {modelist}')
963
964 # need to set default for nperseg before setting default for noverlap below
965 window, nperseg = _triage_segments(window, nperseg,
966 input_length=x.shape[axis])
967
968 # Less overlap than welch, so samples are more statistically independent
969 if noverlap is None:
970 noverlap = nperseg // 8
971
972 if mode == 'psd':
973 freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg,
974 noverlap, nfft, detrend,
975 return_onesided, scaling, axis,
976 mode='psd')
977
978 else:
979 freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg,
980 noverlap, nfft, detrend,
981 return_onesided, scaling, axis,
982 mode='stft')
983
984 if mode == 'magnitude':
985 Sxx = np.abs(Sxx)
986 elif mode in ['angle', 'phase']:
987 Sxx = np.angle(Sxx)
988 if mode == 'phase':
989 # Sxx has one additional dimension for time strides
990 if axis < 0:
991 axis -= 1
992 Sxx = np.unwrap(Sxx, axis=axis)
993
994 # mode =='complex' is same as `stft`, doesn't need modification
995
996 return freqs, time, Sxx
997
998
999def check_COLA(window, nperseg, noverlap, tol=1e-10):
1000 r"""Check whether the Constant OverLap Add (COLA) constraint is met.
1001
1002 Parameters
1003 ----------
1004 window : str or tuple or array_like
1005 Desired window to use. If `window` is a string or tuple, it is
1006 passed to `get_window` to generate the window values, which are
1007 DFT-even by default. See `get_window` for a list of windows and
1008 required parameters. If `window` is array_like it will be used
1009 directly as the window and its length must be nperseg.
1010 nperseg : int
1011 Length of each segment.
1012 noverlap : int
1013 Number of points to overlap between segments.
1014 tol : float, optional
1015 The allowed variance of a bin's weighted sum from the median bin
1016 sum.
1017
1018 Returns
1019 -------
1020 verdict : bool
1021 `True` if chosen combination satisfies COLA within `tol`,
1022 `False` otherwise
1023
1024 See Also
1025 --------
1026 check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met
1027 stft: Short Time Fourier Transform
1028 istft: Inverse Short Time Fourier Transform
1029
1030 Notes
1031 -----
1032 In order to enable inversion of an STFT via the inverse STFT in
1033 `istft`, it is sufficient that the signal windowing obeys the constraint of
1034 "Constant OverLap Add" (COLA). This ensures that every point in the input
1035 data is equally weighted, thereby avoiding aliasing and allowing full
1036 reconstruction.
1037
1038 Some examples of windows that satisfy COLA:
1039 - Rectangular window at overlap of 0, 1/2, 2/3, 3/4, ...
1040 - Bartlett window at overlap of 1/2, 3/4, 5/6, ...
1041 - Hann window at 1/2, 2/3, 3/4, ...
1042 - Any Blackman family window at 2/3 overlap
1043 - Any window with ``noverlap = nperseg-1``
1044
1045 A very comprehensive list of other windows may be found in [2]_,
1046 wherein the COLA condition is satisfied when the "Amplitude
1047 Flatness" is unity.
1048
1049 .. versionadded:: 0.19.0
1050
1051 References
1052 ----------
1053 .. [1] Julius O. Smith III, "Spectral Audio Signal Processing", W3K
1054 Publishing, 2011,ISBN 978-0-9745607-3-1.
1055 .. [2] G. Heinzel, A. Ruediger and R. Schilling, "Spectrum and
1056 spectral density estimation by the Discrete Fourier transform
1057 (DFT), including a comprehensive list of window functions and
1058 some new at-top windows", 2002,
1059 http://hdl.handle.net/11858/00-001M-0000-0013-557A-5
1060
1061 Examples
1062 --------
1063 >>> from scipy import signal
1064
1065 Confirm COLA condition for rectangular window of 75% (3/4) overlap:
1066
1067 >>> signal.check_COLA(signal.windows.boxcar(100), 100, 75)
1068 True
1069
1070 COLA is not true for 25% (1/4) overlap, though:
1071
1072 >>> signal.check_COLA(signal.windows.boxcar(100), 100, 25)
1073 False
1074
1075 "Symmetrical" Hann window (for filter design) is not COLA:
1076
1077 >>> signal.check_COLA(signal.windows.hann(120, sym=True), 120, 60)
1078 False
1079
1080 "Periodic" or "DFT-even" Hann window (for FFT analysis) is COLA for
1081 overlap of 1/2, 2/3, 3/4, etc.:
1082
1083 >>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 60)
1084 True
1085
1086 >>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 80)
1087 True
1088
1089 >>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 90)
1090 True
1091
1092 """
1093 nperseg = int(nperseg)
1094
1095 if nperseg < 1:
1096 raise ValueError('nperseg must be a positive integer')
1097
1098 if noverlap >= nperseg:
1099 raise ValueError('noverlap must be less than nperseg.')
1100 noverlap = int(noverlap)
1101
1102 if isinstance(window, str) or type(window) is tuple:
1103 win = get_window(window, nperseg)
1104 else:
1105 win = np.asarray(window)
1106 if len(win.shape) != 1:
1107 raise ValueError('window must be 1-D')
1108 if win.shape[0] != nperseg:
1109 raise ValueError('window must have length of nperseg')
1110
1111 step = nperseg - noverlap
1112 binsums = sum(win[ii*step:(ii+1)*step] for ii in range(nperseg//step))
1113
1114 if nperseg % step != 0:
1115 binsums[:nperseg % step] += win[-(nperseg % step):]
1116
1117 deviation = binsums - np.median(binsums)
1118 return np.max(np.abs(deviation)) < tol
1119
1120
1121def check_NOLA(window, nperseg, noverlap, tol=1e-10):
1122 r"""Check whether the Nonzero Overlap Add (NOLA) constraint is met.
1123
1124 Parameters
1125 ----------
1126 window : str or tuple or array_like
1127 Desired window to use. If `window` is a string or tuple, it is
1128 passed to `get_window` to generate the window values, which are
1129 DFT-even by default. See `get_window` for a list of windows and
1130 required parameters. If `window` is array_like it will be used
1131 directly as the window and its length must be nperseg.
1132 nperseg : int
1133 Length of each segment.
1134 noverlap : int
1135 Number of points to overlap between segments.
1136 tol : float, optional
1137 The allowed variance of a bin's weighted sum from the median bin
1138 sum.
1139
1140 Returns
1141 -------
1142 verdict : bool
1143 `True` if chosen combination satisfies the NOLA constraint within
1144 `tol`, `False` otherwise
1145
1146 See Also
1147 --------
1148 check_COLA: Check whether the Constant OverLap Add (COLA) constraint is met
1149 stft: Short Time Fourier Transform
1150 istft: Inverse Short Time Fourier Transform
1151
1152 Notes
1153 -----
1154 In order to enable inversion of an STFT via the inverse STFT in
1155 `istft`, the signal windowing must obey the constraint of "nonzero
1156 overlap add" (NOLA):
1157
1158 .. math:: \sum_{t}w^{2}[n-tH] \ne 0
1159
1160 for all :math:`n`, where :math:`w` is the window function, :math:`t` is the
1161 frame index, and :math:`H` is the hop size (:math:`H` = `nperseg` -
1162 `noverlap`).
1163
1164 This ensures that the normalization factors in the denominator of the
1165 overlap-add inversion equation are not zero. Only very pathological windows
1166 will fail the NOLA constraint.
1167
1168 .. versionadded:: 1.2.0
1169
1170 References
1171 ----------
1172 .. [1] Julius O. Smith III, "Spectral Audio Signal Processing", W3K
1173 Publishing, 2011,ISBN 978-0-9745607-3-1.
1174 .. [2] G. Heinzel, A. Ruediger and R. Schilling, "Spectrum and
1175 spectral density estimation by the Discrete Fourier transform
1176 (DFT), including a comprehensive list of window functions and
1177 some new at-top windows", 2002,
1178 http://hdl.handle.net/11858/00-001M-0000-0013-557A-5
1179
1180 Examples
1181 --------
1182 >>> import numpy as np
1183 >>> from scipy import signal
1184
1185 Confirm NOLA condition for rectangular window of 75% (3/4) overlap:
1186
1187 >>> signal.check_NOLA(signal.windows.boxcar(100), 100, 75)
1188 True
1189
1190 NOLA is also true for 25% (1/4) overlap:
1191
1192 >>> signal.check_NOLA(signal.windows.boxcar(100), 100, 25)
1193 True
1194
1195 "Symmetrical" Hann window (for filter design) is also NOLA:
1196
1197 >>> signal.check_NOLA(signal.windows.hann(120, sym=True), 120, 60)
1198 True
1199
1200 As long as there is overlap, it takes quite a pathological window to fail