Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Tools for spectral analysis. 

2""" 

3 

4import numpy as np 

5from scipy import fft as sp_fft 

6from . import signaltools 

7from .windows import get_window 

8from ._spectral import _lombscargle 

9from ._arraytools import const_ext, even_ext, odd_ext, zero_ext 

10import warnings 

11 

12 

13__all__ = ['periodogram', 'welch', 'lombscargle', 'csd', 'coherence', 

14 'spectrogram', 'stft', 'istft', 'check_COLA', 'check_NOLA'] 

15 

16 

17def lombscargle(x, 

18 y, 

19 freqs, 

20 precenter=False, 

21 normalize=False): 

22 """ 

23 lombscargle(x, y, freqs) 

24 

25 Computes the Lomb-Scargle periodogram. 

26 

27 The Lomb-Scargle periodogram was developed by Lomb [1]_ and further 

28 extended by Scargle [2]_ to find, and test the significance of weak 

29 periodic signals with uneven temporal sampling. 

30 

31 When *normalize* is False (default) the computed periodogram 

32 is unnormalized, it takes the value ``(A**2) * N/4`` for a harmonic 

33 signal with amplitude A for sufficiently large N. 

34 

35 When *normalize* is True the computed periodogram is normalized by 

36 the residuals of the data around a constant reference model (at zero). 

37 

38 Input arrays should be 1-D and will be cast to float64. 

39 

40 Parameters 

41 ---------- 

42 x : array_like 

43 Sample times. 

44 y : array_like 

45 Measurement values. 

46 freqs : array_like 

47 Angular frequencies for output periodogram. 

48 precenter : bool, optional 

49 Pre-center amplitudes by subtracting the mean. 

50 normalize : bool, optional 

51 Compute normalized periodogram. 

52 

53 Returns 

54 ------- 

55 pgram : array_like 

56 Lomb-Scargle periodogram. 

57 

58 Raises 

59 ------ 

60 ValueError 

61 If the input arrays `x` and `y` do not have the same shape. 

62 

63 Notes 

64 ----- 

65 This subroutine calculates the periodogram using a slightly 

66 modified algorithm due to Townsend [3]_ which allows the 

67 periodogram to be calculated using only a single pass through 

68 the input arrays for each frequency. 

69 

70 The algorithm running time scales roughly as O(x * freqs) or O(N^2) 

71 for a large number of samples and frequencies. 

72 

73 References 

74 ---------- 

75 .. [1] N.R. Lomb "Least-squares frequency analysis of unequally spaced 

76 data", Astrophysics and Space Science, vol 39, pp. 447-462, 1976 

77 

78 .. [2] J.D. Scargle "Studies in astronomical time series analysis. II - 

79 Statistical aspects of spectral analysis of unevenly spaced data", 

80 The Astrophysical Journal, vol 263, pp. 835-853, 1982 

81 

82 .. [3] R.H.D. Townsend, "Fast calculation of the Lomb-Scargle 

83 periodogram using graphics processing units.", The Astrophysical 

84 Journal Supplement Series, vol 191, pp. 247-253, 2010 

85 

86 See Also 

87 -------- 

88 istft: Inverse Short Time Fourier Transform 

89 check_COLA: Check whether the Constant OverLap Add (COLA) constraint is met 

90 welch: Power spectral density by Welch's method 

91 spectrogram: Spectrogram by Welch's method 

92 csd: Cross spectral density by Welch's method 

93 

94 Examples 

95 -------- 

96 >>> import matplotlib.pyplot as plt 

97 

98 First define some input parameters for the signal: 

99 

100 >>> A = 2. 

101 >>> w = 1. 

102 >>> phi = 0.5 * np.pi 

103 >>> nin = 1000 

104 >>> nout = 100000 

105 >>> frac_points = 0.9 # Fraction of points to select 

106 

107 Randomly select a fraction of an array with timesteps: 

108 

109 >>> r = np.random.rand(nin) 

110 >>> x = np.linspace(0.01, 10*np.pi, nin) 

111 >>> x = x[r >= frac_points] 

112 

113 Plot a sine wave for the selected times: 

114 

115 >>> y = A * np.sin(w*x+phi) 

116 

117 Define the array of frequencies for which to compute the periodogram: 

118 

119 >>> f = np.linspace(0.01, 10, nout) 

120 

121 Calculate Lomb-Scargle periodogram: 

122 

123 >>> import scipy.signal as signal 

124 >>> pgram = signal.lombscargle(x, y, f, normalize=True) 

125 

126 Now make a plot of the input data: 

127 

128 >>> plt.subplot(2, 1, 1) 

129 >>> plt.plot(x, y, 'b+') 

130 

131 Then plot the normalized periodogram: 

132 

133 >>> plt.subplot(2, 1, 2) 

134 >>> plt.plot(f, pgram) 

135 >>> plt.show() 

136 

137 """ 

138 

139 x = np.asarray(x, dtype=np.float64) 

140 y = np.asarray(y, dtype=np.float64) 

141 freqs = np.asarray(freqs, dtype=np.float64) 

142 

143 assert x.ndim == 1 

144 assert y.ndim == 1 

145 assert freqs.ndim == 1 

146 

147 if precenter: 

148 pgram = _lombscargle(x, y - y.mean(), freqs) 

149 else: 

150 pgram = _lombscargle(x, y, freqs) 

151 

152 if normalize: 

153 pgram *= 2 / np.dot(y, y) 

154 

155 return pgram 

156 

157 

158def periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', 

159 return_onesided=True, scaling='density', axis=-1): 

160 """ 

161 Estimate power spectral density using a periodogram. 

162 

163 Parameters 

164 ---------- 

165 x : array_like 

166 Time series of measurement values 

167 fs : float, optional 

168 Sampling frequency of the `x` time series. Defaults to 1.0. 

169 window : str or tuple or array_like, optional 

170 Desired window to use. If `window` is a string or tuple, it is 

171 passed to `get_window` to generate the window values, which are 

172 DFT-even by default. See `get_window` for a list of windows and 

173 required parameters. If `window` is array_like it will be used 

174 directly as the window and its length must be nperseg. Defaults 

175 to 'boxcar'. 

176 nfft : int, optional 

177 Length of the FFT used. If `None` the length of `x` will be 

178 used. 

179 detrend : str or function or `False`, optional 

180 Specifies how to detrend each segment. If `detrend` is a 

181 string, it is passed as the `type` argument to the `detrend` 

182 function. If it is a function, it takes a segment and returns a 

183 detrended segment. If `detrend` is `False`, no detrending is 

184 done. Defaults to 'constant'. 

185 return_onesided : bool, optional 

186 If `True`, return a one-sided spectrum for real data. If 

187 `False` return a two-sided spectrum. Defaults to `True`, but for 

188 complex data, a two-sided spectrum is always returned. 

189 scaling : { 'density', 'spectrum' }, optional 

190 Selects between computing the power spectral density ('density') 

191 where `Pxx` has units of V**2/Hz and computing the power 

192 spectrum ('spectrum') where `Pxx` has units of V**2, if `x` 

193 is measured in V and `fs` is measured in Hz. Defaults to 

194 'density' 

195 axis : int, optional 

196 Axis along which the periodogram is computed; the default is 

197 over the last axis (i.e. ``axis=-1``). 

198 

199 Returns 

200 ------- 

201 f : ndarray 

202 Array of sample frequencies. 

203 Pxx : ndarray 

204 Power spectral density or power spectrum of `x`. 

205 

206 Notes 

207 ----- 

208 .. versionadded:: 0.12.0 

209 

210 See Also 

211 -------- 

212 welch: Estimate power spectral density using Welch's method 

213 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

214 

215 Examples 

216 -------- 

217 >>> from scipy import signal 

218 >>> import matplotlib.pyplot as plt 

219 >>> np.random.seed(1234) 

220 

221 Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 

222 0.001 V**2/Hz of white noise sampled at 10 kHz. 

223 

224 >>> fs = 10e3 

225 >>> N = 1e5 

226 >>> amp = 2*np.sqrt(2) 

227 >>> freq = 1234.0 

228 >>> noise_power = 0.001 * fs / 2 

229 >>> time = np.arange(N) / fs 

230 >>> x = amp*np.sin(2*np.pi*freq*time) 

231 >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

232 

233 Compute and plot the power spectral density. 

234 

235 >>> f, Pxx_den = signal.periodogram(x, fs) 

236 >>> plt.semilogy(f, Pxx_den) 

237 >>> plt.ylim([1e-7, 1e2]) 

238 >>> plt.xlabel('frequency [Hz]') 

239 >>> plt.ylabel('PSD [V**2/Hz]') 

240 >>> plt.show() 

241 

242 If we average the last half of the spectral density, to exclude the 

243 peak, we can recover the noise power on the signal. 

244 

245 >>> np.mean(Pxx_den[25000:]) 

246 0.00099728892368242854 

247 

248 Now compute and plot the power spectrum. 

249 

250 >>> f, Pxx_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum') 

251 >>> plt.figure() 

252 >>> plt.semilogy(f, np.sqrt(Pxx_spec)) 

253 >>> plt.ylim([1e-4, 1e1]) 

254 >>> plt.xlabel('frequency [Hz]') 

255 >>> plt.ylabel('Linear spectrum [V RMS]') 

256 >>> plt.show() 

257 

258 The peak height in the power spectrum is an estimate of the RMS 

259 amplitude. 

260 

261 >>> np.sqrt(Pxx_spec.max()) 

262 2.0077340678640727 

263 

264 """ 

265 x = np.asarray(x) 

266 

267 if x.size == 0: 

268 return np.empty(x.shape), np.empty(x.shape) 

269 

270 if window is None: 

271 window = 'boxcar' 

272 

273 if nfft is None: 

274 nperseg = x.shape[axis] 

275 elif nfft == x.shape[axis]: 

276 nperseg = nfft 

277 elif nfft > x.shape[axis]: 

278 nperseg = x.shape[axis] 

279 elif nfft < x.shape[axis]: 

280 s = [np.s_[:]]*len(x.shape) 

281 s[axis] = np.s_[:nfft] 

282 x = x[tuple(s)] 

283 nperseg = nfft 

284 nfft = None 

285 

286 return welch(x, fs=fs, window=window, nperseg=nperseg, noverlap=0, 

287 nfft=nfft, detrend=detrend, return_onesided=return_onesided, 

288 scaling=scaling, axis=axis) 

289 

290 

291def welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, 

292 detrend='constant', return_onesided=True, scaling='density', 

293 axis=-1, average='mean'): 

294 r""" 

295 Estimate power spectral density using Welch's method. 

296 

297 Welch's method [1]_ computes an estimate of the power spectral 

298 density by dividing the data into overlapping segments, computing a 

299 modified periodogram for each segment and averaging the 

300 periodograms. 

301 

302 Parameters 

303 ---------- 

304 x : array_like 

305 Time series of measurement values 

306 fs : float, optional 

307 Sampling frequency of the `x` time series. Defaults to 1.0. 

308 window : str or tuple or array_like, optional 

309 Desired window to use. If `window` is a string or tuple, it is 

310 passed to `get_window` to generate the window values, which are 

311 DFT-even by default. See `get_window` for a list of windows and 

312 required parameters. If `window` is array_like it will be used 

313 directly as the window and its length must be nperseg. Defaults 

314 to a Hann window. 

315 nperseg : int, optional 

316 Length of each segment. Defaults to None, but if window is str or 

317 tuple, is set to 256, and if window is array_like, is set to the 

318 length of the window. 

319 noverlap : int, optional 

320 Number of points to overlap between segments. If `None`, 

321 ``noverlap = nperseg // 2``. Defaults to `None`. 

322 nfft : int, optional 

323 Length of the FFT used, if a zero padded FFT is desired. If 

324 `None`, the FFT length is `nperseg`. Defaults to `None`. 

325 detrend : str or function or `False`, optional 

326 Specifies how to detrend each segment. If `detrend` is a 

327 string, it is passed as the `type` argument to the `detrend` 

328 function. If it is a function, it takes a segment and returns a 

329 detrended segment. If `detrend` is `False`, no detrending is 

330 done. Defaults to 'constant'. 

331 return_onesided : bool, optional 

332 If `True`, return a one-sided spectrum for real data. If 

333 `False` return a two-sided spectrum. Defaults to `True`, but for 

334 complex data, a two-sided spectrum is always returned. 

335 scaling : { 'density', 'spectrum' }, optional 

336 Selects between computing the power spectral density ('density') 

337 where `Pxx` has units of V**2/Hz and computing the power 

338 spectrum ('spectrum') where `Pxx` has units of V**2, if `x` 

339 is measured in V and `fs` is measured in Hz. Defaults to 

340 'density' 

341 axis : int, optional 

342 Axis along which the periodogram is computed; the default is 

343 over the last axis (i.e. ``axis=-1``). 

344 average : { 'mean', 'median' }, optional 

345 Method to use when averaging periodograms. Defaults to 'mean'. 

346 

347 .. versionadded:: 1.2.0 

348 

349 Returns 

350 ------- 

351 f : ndarray 

352 Array of sample frequencies. 

353 Pxx : ndarray 

354 Power spectral density or power spectrum of x. 

355 

356 See Also 

357 -------- 

358 periodogram: Simple, optionally modified periodogram 

359 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

360 

361 Notes 

362 ----- 

363 An appropriate amount of overlap will depend on the choice of window 

364 and on your requirements. For the default Hann window an overlap of 

365 50% is a reasonable trade off between accurately estimating the 

366 signal power, while not over counting any of the data. Narrower 

367 windows may require a larger overlap. 

368 

369 If `noverlap` is 0, this method is equivalent to Bartlett's method 

370 [2]_. 

371 

372 .. versionadded:: 0.12.0 

373 

374 References 

375 ---------- 

376 .. [1] P. Welch, "The use of the fast Fourier transform for the 

377 estimation of power spectra: A method based on time averaging 

378 over short, modified periodograms", IEEE Trans. Audio 

379 Electroacoust. vol. 15, pp. 70-73, 1967. 

380 .. [2] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", 

381 Biometrika, vol. 37, pp. 1-16, 1950. 

382 

383 Examples 

384 -------- 

385 >>> from scipy import signal 

386 >>> import matplotlib.pyplot as plt 

387 >>> np.random.seed(1234) 

388 

389 Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 

390 0.001 V**2/Hz of white noise sampled at 10 kHz. 

391 

392 >>> fs = 10e3 

393 >>> N = 1e5 

394 >>> amp = 2*np.sqrt(2) 

395 >>> freq = 1234.0 

396 >>> noise_power = 0.001 * fs / 2 

397 >>> time = np.arange(N) / fs 

398 >>> x = amp*np.sin(2*np.pi*freq*time) 

399 >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

400 

401 Compute and plot the power spectral density. 

402 

403 >>> f, Pxx_den = signal.welch(x, fs, nperseg=1024) 

404 >>> plt.semilogy(f, Pxx_den) 

405 >>> plt.ylim([0.5e-3, 1]) 

406 >>> plt.xlabel('frequency [Hz]') 

407 >>> plt.ylabel('PSD [V**2/Hz]') 

408 >>> plt.show() 

409 

410 If we average the last half of the spectral density, to exclude the 

411 peak, we can recover the noise power on the signal. 

412 

413 >>> np.mean(Pxx_den[256:]) 

414 0.0009924865443739191 

415 

416 Now compute and plot the power spectrum. 

417 

418 >>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum') 

419 >>> plt.figure() 

420 >>> plt.semilogy(f, np.sqrt(Pxx_spec)) 

421 >>> plt.xlabel('frequency [Hz]') 

422 >>> plt.ylabel('Linear spectrum [V RMS]') 

423 >>> plt.show() 

424 

425 The peak height in the power spectrum is an estimate of the RMS 

426 amplitude. 

427 

428 >>> np.sqrt(Pxx_spec.max()) 

429 2.0077340678640727 

430 

431 If we now introduce a discontinuity in the signal, by increasing the 

432 amplitude of a small portion of the signal by 50, we can see the 

433 corruption of the mean average power spectral density, but using a 

434 median average better estimates the normal behaviour. 

435 

436 >>> x[int(N//2):int(N//2)+10] *= 50. 

437 >>> f, Pxx_den = signal.welch(x, fs, nperseg=1024) 

438 >>> f_med, Pxx_den_med = signal.welch(x, fs, nperseg=1024, average='median') 

439 >>> plt.semilogy(f, Pxx_den, label='mean') 

440 >>> plt.semilogy(f_med, Pxx_den_med, label='median') 

441 >>> plt.ylim([0.5e-3, 1]) 

442 >>> plt.xlabel('frequency [Hz]') 

443 >>> plt.ylabel('PSD [V**2/Hz]') 

444 >>> plt.legend() 

445 >>> plt.show() 

446 

447 """ 

448 

449 freqs, Pxx = csd(x, x, fs=fs, window=window, nperseg=nperseg, 

450 noverlap=noverlap, nfft=nfft, detrend=detrend, 

451 return_onesided=return_onesided, scaling=scaling, 

452 axis=axis, average=average) 

453 

454 return freqs, Pxx.real 

455 

456 

457def csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, 

458 detrend='constant', return_onesided=True, scaling='density', 

459 axis=-1, average='mean'): 

460 r""" 

461 Estimate the cross power spectral density, Pxy, using Welch's 

462 method. 

463 

464 Parameters 

465 ---------- 

466 x : array_like 

467 Time series of measurement values 

468 y : array_like 

469 Time series of measurement values 

470 fs : float, optional 

471 Sampling frequency of the `x` and `y` time series. Defaults 

472 to 1.0. 

473 window : str or tuple or array_like, optional 

474 Desired window to use. If `window` is a string or tuple, it is 

475 passed to `get_window` to generate the window values, which are 

476 DFT-even by default. See `get_window` for a list of windows and 

477 required parameters. If `window` is array_like it will be used 

478 directly as the window and its length must be nperseg. Defaults 

479 to a Hann window. 

480 nperseg : int, optional 

481 Length of each segment. Defaults to None, but if window is str or 

482 tuple, is set to 256, and if window is array_like, is set to the 

483 length of the window. 

484 noverlap: int, optional 

485 Number of points to overlap between segments. If `None`, 

486 ``noverlap = nperseg // 2``. Defaults to `None`. 

487 nfft : int, optional 

488 Length of the FFT used, if a zero padded FFT is desired. If 

489 `None`, the FFT length is `nperseg`. Defaults to `None`. 

490 detrend : str or function or `False`, optional 

491 Specifies how to detrend each segment. If `detrend` is a 

492 string, it is passed as the `type` argument to the `detrend` 

493 function. If it is a function, it takes a segment and returns a 

494 detrended segment. If `detrend` is `False`, no detrending is 

495 done. Defaults to 'constant'. 

496 return_onesided : bool, optional 

497 If `True`, return a one-sided spectrum for real data. If 

498 `False` return a two-sided spectrum. Defaults to `True`, but for 

499 complex data, a two-sided spectrum is always returned. 

500 scaling : { 'density', 'spectrum' }, optional 

501 Selects between computing the cross spectral density ('density') 

502 where `Pxy` has units of V**2/Hz and computing the cross spectrum 

503 ('spectrum') where `Pxy` has units of V**2, if `x` and `y` are 

504 measured in V and `fs` is measured in Hz. Defaults to 'density' 

505 axis : int, optional 

506 Axis along which the CSD is computed for both inputs; the 

507 default is over the last axis (i.e. ``axis=-1``). 

508 average : { 'mean', 'median' }, optional 

509 Method to use when averaging periodograms. Defaults to 'mean'. 

510 

511 .. versionadded:: 1.2.0 

512 

513 Returns 

514 ------- 

515 f : ndarray 

516 Array of sample frequencies. 

517 Pxy : ndarray 

518 Cross spectral density or cross power spectrum of x,y. 

519 

520 See Also 

521 -------- 

522 periodogram: Simple, optionally modified periodogram 

523 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

524 welch: Power spectral density by Welch's method. [Equivalent to 

525 csd(x,x)] 

526 coherence: Magnitude squared coherence by Welch's method. 

527 

528 Notes 

529 -------- 

530 By convention, Pxy is computed with the conjugate FFT of X 

531 multiplied by the FFT of Y. 

532 

533 If the input series differ in length, the shorter series will be 

534 zero-padded to match. 

535 

536 An appropriate amount of overlap will depend on the choice of window 

537 and on your requirements. For the default Hann window an overlap of 

538 50% is a reasonable trade off between accurately estimating the 

539 signal power, while not over counting any of the data. Narrower 

540 windows may require a larger overlap. 

541 

542 .. versionadded:: 0.16.0 

543 

544 References 

545 ---------- 

546 .. [1] P. Welch, "The use of the fast Fourier transform for the 

547 estimation of power spectra: A method based on time averaging 

548 over short, modified periodograms", IEEE Trans. Audio 

549 Electroacoust. vol. 15, pp. 70-73, 1967. 

550 .. [2] Rabiner, Lawrence R., and B. Gold. "Theory and Application of 

551 Digital Signal Processing" Prentice-Hall, pp. 414-419, 1975 

552 

553 Examples 

554 -------- 

555 >>> from scipy import signal 

556 >>> import matplotlib.pyplot as plt 

557 

558 Generate two test signals with some common features. 

559 

560 >>> fs = 10e3 

561 >>> N = 1e5 

562 >>> amp = 20 

563 >>> freq = 1234.0 

564 >>> noise_power = 0.001 * fs / 2 

565 >>> time = np.arange(N) / fs 

566 >>> b, a = signal.butter(2, 0.25, 'low') 

567 >>> x = np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

568 >>> y = signal.lfilter(b, a, x) 

569 >>> x += amp*np.sin(2*np.pi*freq*time) 

570 >>> y += np.random.normal(scale=0.1*np.sqrt(noise_power), size=time.shape) 

571 

572 Compute and plot the magnitude of the cross spectral density. 

573 

574 >>> f, Pxy = signal.csd(x, y, fs, nperseg=1024) 

575 >>> plt.semilogy(f, np.abs(Pxy)) 

576 >>> plt.xlabel('frequency [Hz]') 

577 >>> plt.ylabel('CSD [V**2/Hz]') 

578 >>> plt.show() 

579 """ 

580 

581 freqs, _, Pxy = _spectral_helper(x, y, fs, window, nperseg, noverlap, nfft, 

582 detrend, return_onesided, scaling, axis, 

583 mode='psd') 

584 

585 # Average over windows. 

586 if len(Pxy.shape) >= 2 and Pxy.size > 0: 

587 if Pxy.shape[-1] > 1: 

588 if average == 'median': 

589 Pxy = np.median(Pxy, axis=-1) / _median_bias(Pxy.shape[-1]) 

590 elif average == 'mean': 

591 Pxy = Pxy.mean(axis=-1) 

592 else: 

593 raise ValueError('average must be "median" or "mean", got %s' 

594 % (average,)) 

595 else: 

596 Pxy = np.reshape(Pxy, Pxy.shape[:-1]) 

597 

598 return freqs, Pxy 

599 

600 

601def spectrogram(x, fs=1.0, window=('tukey', .25), nperseg=None, noverlap=None, 

602 nfft=None, detrend='constant', return_onesided=True, 

603 scaling='density', axis=-1, mode='psd'): 

604 """ 

605 Compute a spectrogram with consecutive Fourier transforms. 

606 

607 Spectrograms can be used as a way of visualizing the change of a 

608 nonstationary signal's frequency content over time. 

609 

610 Parameters 

611 ---------- 

612 x : array_like 

613 Time series of measurement values 

614 fs : float, optional 

615 Sampling frequency of the `x` time series. Defaults to 1.0. 

616 window : str or tuple or array_like, optional 

617 Desired window to use. If `window` is a string or tuple, it is 

618 passed to `get_window` to generate the window values, which are 

619 DFT-even by default. See `get_window` for a list of windows and 

620 required parameters. If `window` is array_like it will be used 

621 directly as the window and its length must be nperseg. 

622 Defaults to a Tukey window with shape parameter of 0.25. 

623 nperseg : int, optional 

624 Length of each segment. Defaults to None, but if window is str or 

625 tuple, is set to 256, and if window is array_like, is set to the 

626 length of the window. 

627 noverlap : int, optional 

628 Number of points to overlap between segments. If `None`, 

629 ``noverlap = nperseg // 8``. Defaults to `None`. 

630 nfft : int, optional 

631 Length of the FFT used, if a zero padded FFT is desired. If 

632 `None`, the FFT length is `nperseg`. Defaults to `None`. 

633 detrend : str or function or `False`, optional 

634 Specifies how to detrend each segment. If `detrend` is a 

635 string, it is passed as the `type` argument to the `detrend` 

636 function. If it is a function, it takes a segment and returns a 

637 detrended segment. If `detrend` is `False`, no detrending is 

638 done. Defaults to 'constant'. 

639 return_onesided : bool, optional 

640 If `True`, return a one-sided spectrum for real data. If 

641 `False` return a two-sided spectrum. Defaults to `True`, but for 

642 complex data, a two-sided spectrum is always returned. 

643 scaling : { 'density', 'spectrum' }, optional 

644 Selects between computing the power spectral density ('density') 

645 where `Sxx` has units of V**2/Hz and computing the power 

646 spectrum ('spectrum') where `Sxx` has units of V**2, if `x` 

647 is measured in V and `fs` is measured in Hz. Defaults to 

648 'density'. 

649 axis : int, optional 

650 Axis along which the spectrogram is computed; the default is over 

651 the last axis (i.e. ``axis=-1``). 

652 mode : str, optional 

653 Defines what kind of return values are expected. Options are 

654 ['psd', 'complex', 'magnitude', 'angle', 'phase']. 'complex' is 

655 equivalent to the output of `stft` with no padding or boundary 

656 extension. 'magnitude' returns the absolute magnitude of the 

657 STFT. 'angle' and 'phase' return the complex angle of the STFT, 

658 with and without unwrapping, respectively. 

659 

660 Returns 

661 ------- 

662 f : ndarray 

663 Array of sample frequencies. 

664 t : ndarray 

665 Array of segment times. 

666 Sxx : ndarray 

667 Spectrogram of x. By default, the last axis of Sxx corresponds 

668 to the segment times. 

669 

670 See Also 

671 -------- 

672 periodogram: Simple, optionally modified periodogram 

673 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

674 welch: Power spectral density by Welch's method. 

675 csd: Cross spectral density by Welch's method. 

676 

677 Notes 

678 ----- 

679 An appropriate amount of overlap will depend on the choice of window 

680 and on your requirements. In contrast to welch's method, where the 

681 entire data stream is averaged over, one may wish to use a smaller 

682 overlap (or perhaps none at all) when computing a spectrogram, to 

683 maintain some statistical independence between individual segments. 

684 It is for this reason that the default window is a Tukey window with 

685 1/8th of a window's length overlap at each end. 

686 

687 .. versionadded:: 0.16.0 

688 

689 References 

690 ---------- 

691 .. [1] Oppenheim, Alan V., Ronald W. Schafer, John R. Buck 

692 "Discrete-Time Signal Processing", Prentice Hall, 1999. 

693 

694 Examples 

695 -------- 

696 >>> from scipy import signal 

697 >>> from scipy.fft import fftshift 

698 >>> import matplotlib.pyplot as plt 

699 

700 Generate a test signal, a 2 Vrms sine wave whose frequency is slowly 

701 modulated around 3kHz, corrupted by white noise of exponentially 

702 decreasing magnitude sampled at 10 kHz. 

703 

704 >>> fs = 10e3 

705 >>> N = 1e5 

706 >>> amp = 2 * np.sqrt(2) 

707 >>> noise_power = 0.01 * fs / 2 

708 >>> time = np.arange(N) / float(fs) 

709 >>> mod = 500*np.cos(2*np.pi*0.25*time) 

710 >>> carrier = amp * np.sin(2*np.pi*3e3*time + mod) 

711 >>> noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

712 >>> noise *= np.exp(-time/5) 

713 >>> x = carrier + noise 

714 

715 Compute and plot the spectrogram. 

716 

717 >>> f, t, Sxx = signal.spectrogram(x, fs) 

718 >>> plt.pcolormesh(t, f, Sxx, shading='gouraud') 

719 >>> plt.ylabel('Frequency [Hz]') 

720 >>> plt.xlabel('Time [sec]') 

721 >>> plt.show() 

722 

723 Note, if using output that is not one sided, then use the following: 

724 

725 >>> f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False) 

726 >>> plt.pcolormesh(t, fftshift(f), fftshift(Sxx, axes=0), shading='gouraud') 

727 >>> plt.ylabel('Frequency [Hz]') 

728 >>> plt.xlabel('Time [sec]') 

729 >>> plt.show() 

730 """ 

731 modelist = ['psd', 'complex', 'magnitude', 'angle', 'phase'] 

732 if mode not in modelist: 

733 raise ValueError('unknown value for mode {}, must be one of {}' 

734 .format(mode, modelist)) 

735 

736 # need to set default for nperseg before setting default for noverlap below 

737 window, nperseg = _triage_segments(window, nperseg, 

738 input_length=x.shape[axis]) 

739 

740 # Less overlap than welch, so samples are more statisically independent 

741 if noverlap is None: 

742 noverlap = nperseg // 8 

743 

744 if mode == 'psd': 

745 freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg, 

746 noverlap, nfft, detrend, 

747 return_onesided, scaling, axis, 

748 mode='psd') 

749 

750 else: 

751 freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg, 

752 noverlap, nfft, detrend, 

753 return_onesided, scaling, axis, 

754 mode='stft') 

755 

756 if mode == 'magnitude': 

757 Sxx = np.abs(Sxx) 

758 elif mode in ['angle', 'phase']: 

759 Sxx = np.angle(Sxx) 

760 if mode == 'phase': 

761 # Sxx has one additional dimension for time strides 

762 if axis < 0: 

763 axis -= 1 

764 Sxx = np.unwrap(Sxx, axis=axis) 

765 

766 # mode =='complex' is same as `stft`, doesn't need modification 

767 

768 return freqs, time, Sxx 

769 

770 

771def check_COLA(window, nperseg, noverlap, tol=1e-10): 

772 r""" 

773 Check whether the Constant OverLap Add (COLA) constraint is met 

774 

775 Parameters 

776 ---------- 

777 window : str or tuple or array_like 

778 Desired window to use. If `window` is a string or tuple, it is 

779 passed to `get_window` to generate the window values, which are 

780 DFT-even by default. See `get_window` for a list of windows and 

781 required parameters. If `window` is array_like it will be used 

782 directly as the window and its length must be nperseg. 

783 nperseg : int 

784 Length of each segment. 

785 noverlap : int 

786 Number of points to overlap between segments. 

787 tol : float, optional 

788 The allowed variance of a bin's weighted sum from the median bin 

789 sum. 

790 

791 Returns 

792 ------- 

793 verdict : bool 

794 `True` if chosen combination satisfies COLA within `tol`, 

795 `False` otherwise 

796 

797 See Also 

798 -------- 

799 check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met 

800 stft: Short Time Fourier Transform 

801 istft: Inverse Short Time Fourier Transform 

802 

803 Notes 

804 ----- 

805 In order to enable inversion of an STFT via the inverse STFT in 

806 `istft`, it is sufficient that the signal windowing obeys the constraint of 

807 "Constant OverLap Add" (COLA). This ensures that every point in the input 

808 data is equally weighted, thereby avoiding aliasing and allowing full 

809 reconstruction. 

810 

811 Some examples of windows that satisfy COLA: 

812 - Rectangular window at overlap of 0, 1/2, 2/3, 3/4, ... 

813 - Bartlett window at overlap of 1/2, 3/4, 5/6, ... 

814 - Hann window at 1/2, 2/3, 3/4, ... 

815 - Any Blackman family window at 2/3 overlap 

816 - Any window with ``noverlap = nperseg-1`` 

817 

818 A very comprehensive list of other windows may be found in [2]_, 

819 wherein the COLA condition is satisfied when the "Amplitude 

820 Flatness" is unity. 

821 

822 .. versionadded:: 0.19.0 

823 

824 References 

825 ---------- 

826 .. [1] Julius O. Smith III, "Spectral Audio Signal Processing", W3K 

827 Publishing, 2011,ISBN 978-0-9745607-3-1. 

828 .. [2] G. Heinzel, A. Ruediger and R. Schilling, "Spectrum and 

829 spectral density estimation by the Discrete Fourier transform 

830 (DFT), including a comprehensive list of window functions and 

831 some new at-top windows", 2002, 

832 http://hdl.handle.net/11858/00-001M-0000-0013-557A-5 

833 

834 Examples 

835 -------- 

836 >>> from scipy import signal 

837 

838 Confirm COLA condition for rectangular window of 75% (3/4) overlap: 

839 

840 >>> signal.check_COLA(signal.boxcar(100), 100, 75) 

841 True 

842 

843 COLA is not true for 25% (1/4) overlap, though: 

844 

845 >>> signal.check_COLA(signal.boxcar(100), 100, 25) 

846 False 

847 

848 "Symmetrical" Hann window (for filter design) is not COLA: 

849 

850 >>> signal.check_COLA(signal.hann(120, sym=True), 120, 60) 

851 False 

852 

853 "Periodic" or "DFT-even" Hann window (for FFT analysis) is COLA for 

854 overlap of 1/2, 2/3, 3/4, etc.: 

855 

856 >>> signal.check_COLA(signal.hann(120, sym=False), 120, 60) 

857 True 

858 

859 >>> signal.check_COLA(signal.hann(120, sym=False), 120, 80) 

860 True 

861 

862 >>> signal.check_COLA(signal.hann(120, sym=False), 120, 90) 

863 True 

864 

865 """ 

866 

867 nperseg = int(nperseg) 

868 

869 if nperseg < 1: 

870 raise ValueError('nperseg must be a positive integer') 

871 

872 if noverlap >= nperseg: 

873 raise ValueError('noverlap must be less than nperseg.') 

874 noverlap = int(noverlap) 

875 

876 if isinstance(window, str) or type(window) is tuple: 

877 win = get_window(window, nperseg) 

878 else: 

879 win = np.asarray(window) 

880 if len(win.shape) != 1: 

881 raise ValueError('window must be 1-D') 

882 if win.shape[0] != nperseg: 

883 raise ValueError('window must have length of nperseg') 

884 

885 step = nperseg - noverlap 

886 binsums = sum(win[ii*step:(ii+1)*step] for ii in range(nperseg//step)) 

887 

888 if nperseg % step != 0: 

889 binsums[:nperseg % step] += win[-(nperseg % step):] 

890 

891 deviation = binsums - np.median(binsums) 

892 return np.max(np.abs(deviation)) < tol 

893 

894 

895def check_NOLA(window, nperseg, noverlap, tol=1e-10): 

896 r""" 

897 Check whether the Nonzero Overlap Add (NOLA) constraint is met 

898 

899 Parameters 

900 ---------- 

901 window : str or tuple or array_like 

902 Desired window to use. If `window` is a string or tuple, it is 

903 passed to `get_window` to generate the window values, which are 

904 DFT-even by default. See `get_window` for a list of windows and 

905 required parameters. If `window` is array_like it will be used 

906 directly as the window and its length must be nperseg. 

907 nperseg : int 

908 Length of each segment. 

909 noverlap : int 

910 Number of points to overlap between segments. 

911 tol : float, optional 

912 The allowed variance of a bin's weighted sum from the median bin 

913 sum. 

914 

915 Returns 

916 ------- 

917 verdict : bool 

918 `True` if chosen combination satisfies the NOLA constraint within 

919 `tol`, `False` otherwise 

920 

921 See Also 

922 -------- 

923 check_COLA: Check whether the Constant OverLap Add (COLA) constraint is met 

924 stft: Short Time Fourier Transform 

925 istft: Inverse Short Time Fourier Transform 

926 

927 Notes 

928 ----- 

929 In order to enable inversion of an STFT via the inverse STFT in 

930 `istft`, the signal windowing must obey the constraint of "nonzero 

931 overlap add" (NOLA): 

932 

933 .. math:: \sum_{t}w^{2}[n-tH] \ne 0 

934 

935 for all :math:`n`, where :math:`w` is the window function, :math:`t` is the 

936 frame index, and :math:`H` is the hop size (:math:`H` = `nperseg` - 

937 `noverlap`). 

938 

939 This ensures that the normalization factors in the denominator of the 

940 overlap-add inversion equation are not zero. Only very pathological windows 

941 will fail the NOLA constraint. 

942 

943 .. versionadded:: 1.2.0 

944 

945 References 

946 ---------- 

947 .. [1] Julius O. Smith III, "Spectral Audio Signal Processing", W3K 

948 Publishing, 2011,ISBN 978-0-9745607-3-1. 

949 .. [2] G. Heinzel, A. Ruediger and R. Schilling, "Spectrum and 

950 spectral density estimation by the Discrete Fourier transform 

951 (DFT), including a comprehensive list of window functions and 

952 some new at-top windows", 2002, 

953 http://hdl.handle.net/11858/00-001M-0000-0013-557A-5 

954 

955 Examples 

956 -------- 

957 >>> from scipy import signal 

958 

959 Confirm NOLA condition for rectangular window of 75% (3/4) overlap: 

960 

961 >>> signal.check_NOLA(signal.boxcar(100), 100, 75) 

962 True 

963 

964 NOLA is also true for 25% (1/4) overlap: 

965 

966 >>> signal.check_NOLA(signal.boxcar(100), 100, 25) 

967 True 

968 

969 "Symmetrical" Hann window (for filter design) is also NOLA: 

970 

971 >>> signal.check_NOLA(signal.hann(120, sym=True), 120, 60) 

972 True 

973 

974 As long as there is overlap, it takes quite a pathological window to fail 

975 NOLA: 

976 

977 >>> w = np.ones(64, dtype="float") 

978 >>> w[::2] = 0 

979 >>> signal.check_NOLA(w, 64, 32) 

980 False 

981 

982 If there is not enough overlap, a window with zeros at the ends will not 

983 work: 

984 

985 >>> signal.check_NOLA(signal.hann(64), 64, 0) 

986 False 

987 >>> signal.check_NOLA(signal.hann(64), 64, 1) 

988 False 

989 >>> signal.check_NOLA(signal.hann(64), 64, 2) 

990 True 

991 """ 

992 

993 nperseg = int(nperseg) 

994 

995 if nperseg < 1: 

996 raise ValueError('nperseg must be a positive integer') 

997 

998 if noverlap >= nperseg: 

999 raise ValueError('noverlap must be less than nperseg') 

1000 if noverlap < 0: 

1001 raise ValueError('noverlap must be a nonnegative integer') 

1002 noverlap = int(noverlap) 

1003 

1004 if isinstance(window, str) or type(window) is tuple: 

1005 win = get_window(window, nperseg) 

1006 else: 

1007 win = np.asarray(window) 

1008 if len(win.shape) != 1: 

1009 raise ValueError('window must be 1-D') 

1010 if win.shape[0] != nperseg: 

1011 raise ValueError('window must have length of nperseg') 

1012 

1013 step = nperseg - noverlap 

1014 binsums = sum(win[ii*step:(ii+1)*step]**2 for ii in range(nperseg//step)) 

1015 

1016 if nperseg % step != 0: 

1017 binsums[:nperseg % step] += win[-(nperseg % step):]**2 

1018 

1019 return np.min(binsums) > tol 

1020 

1021 

1022def stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, 

1023 detrend=False, return_onesided=True, boundary='zeros', padded=True, 

1024 axis=-1): 

1025 r""" 

1026 Compute the Short Time Fourier Transform (STFT). 

1027 

1028 STFTs can be used as a way of quantifying the change of a 

1029 nonstationary signal's frequency and phase content over time. 

1030 

1031 Parameters 

1032 ---------- 

1033 x : array_like 

1034 Time series of measurement values 

1035 fs : float, optional 

1036 Sampling frequency of the `x` time series. Defaults to 1.0. 

1037 window : str or tuple or array_like, optional 

1038 Desired window to use. If `window` is a string or tuple, it is 

1039 passed to `get_window` to generate the window values, which are 

1040 DFT-even by default. See `get_window` for a list of windows and 

1041 required parameters. If `window` is array_like it will be used 

1042 directly as the window and its length must be nperseg. Defaults 

1043 to a Hann window. 

1044 nperseg : int, optional 

1045 Length of each segment. Defaults to 256. 

1046 noverlap : int, optional 

1047 Number of points to overlap between segments. If `None`, 

1048 ``noverlap = nperseg // 2``. Defaults to `None`. When 

1049 specified, the COLA constraint must be met (see Notes below). 

1050 nfft : int, optional 

1051 Length of the FFT used, if a zero padded FFT is desired. If 

1052 `None`, the FFT length is `nperseg`. Defaults to `None`. 

1053 detrend : str or function or `False`, optional 

1054 Specifies how to detrend each segment. If `detrend` is a 

1055 string, it is passed as the `type` argument to the `detrend` 

1056 function. If it is a function, it takes a segment and returns a 

1057 detrended segment. If `detrend` is `False`, no detrending is 

1058 done. Defaults to `False`. 

1059 return_onesided : bool, optional 

1060 If `True`, return a one-sided spectrum for real data. If 

1061 `False` return a two-sided spectrum. Defaults to `True`, but for 

1062 complex data, a two-sided spectrum is always returned. 

1063 boundary : str or None, optional 

1064 Specifies whether the input signal is extended at both ends, and 

1065 how to generate the new values, in order to center the first 

1066 windowed segment on the first input point. This has the benefit 

1067 of enabling reconstruction of the first input point when the 

1068 employed window function starts at zero. Valid options are 

1069 ``['even', 'odd', 'constant', 'zeros', None]``. Defaults to 

1070 'zeros', for zero padding extension. I.e. ``[1, 2, 3, 4]`` is 

1071 extended to ``[0, 1, 2, 3, 4, 0]`` for ``nperseg=3``. 

1072 padded : bool, optional 

1073 Specifies whether the input signal is zero-padded at the end to 

1074 make the signal fit exactly into an integer number of window 

1075 segments, so that all of the signal is included in the output. 

1076 Defaults to `True`. Padding occurs after boundary extension, if 

1077 `boundary` is not `None`, and `padded` is `True`, as is the 

1078 default. 

1079 axis : int, optional 

1080 Axis along which the STFT is computed; the default is over the 

1081 last axis (i.e. ``axis=-1``). 

1082 

1083 Returns 

1084 ------- 

1085 f : ndarray 

1086 Array of sample frequencies. 

1087 t : ndarray 

1088 Array of segment times. 

1089 Zxx : ndarray 

1090 STFT of `x`. By default, the last axis of `Zxx` corresponds 

1091 to the segment times. 

1092 

1093 See Also 

1094 -------- 

1095 istft: Inverse Short Time Fourier Transform 

1096 check_COLA: Check whether the Constant OverLap Add (COLA) constraint 

1097 is met 

1098 check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met 

1099 welch: Power spectral density by Welch's method. 

1100 spectrogram: Spectrogram by Welch's method. 

1101 csd: Cross spectral density by Welch's method. 

1102 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

1103 

1104 Notes 

1105 ----- 

1106 In order to enable inversion of an STFT via the inverse STFT in 

1107 `istft`, the signal windowing must obey the constraint of "Nonzero 

1108 OverLap Add" (NOLA), and the input signal must have complete 

1109 windowing coverage (i.e. ``(x.shape[axis] - nperseg) % 

1110 (nperseg-noverlap) == 0``). The `padded` argument may be used to 

1111 accomplish this. 

1112 

1113 Given a time-domain signal :math:`x[n]`, a window :math:`w[n]`, and a hop 

1114 size :math:`H` = `nperseg - noverlap`, the windowed frame at time index 

1115 :math:`t` is given by 

1116 

1117 .. math:: x_{t}[n]=x[n]w[n-tH] 

1118 

1119 The overlap-add (OLA) reconstruction equation is given by 

1120 

1121 .. math:: x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]} 

1122 

1123 The NOLA constraint ensures that every normalization term that appears 

1124 in the denomimator of the OLA reconstruction equation is nonzero. Whether a 

1125 choice of `window`, `nperseg`, and `noverlap` satisfy this constraint can 

1126 be tested with `check_NOLA`. 

1127 

1128 .. versionadded:: 0.19.0 

1129 

1130 References 

1131 ---------- 

1132 .. [1] Oppenheim, Alan V., Ronald W. Schafer, John R. Buck 

1133 "Discrete-Time Signal Processing", Prentice Hall, 1999. 

1134 .. [2] Daniel W. Griffin, Jae S. Lim "Signal Estimation from 

1135 Modified Short-Time Fourier Transform", IEEE 1984, 

1136 10.1109/TASSP.1984.1164317 

1137 

1138 Examples 

1139 -------- 

1140 >>> from scipy import signal 

1141 >>> import matplotlib.pyplot as plt 

1142 

1143 Generate a test signal, a 2 Vrms sine wave whose frequency is slowly 

1144 modulated around 3kHz, corrupted by white noise of exponentially 

1145 decreasing magnitude sampled at 10 kHz. 

1146 

1147 >>> fs = 10e3 

1148 >>> N = 1e5 

1149 >>> amp = 2 * np.sqrt(2) 

1150 >>> noise_power = 0.01 * fs / 2 

1151 >>> time = np.arange(N) / float(fs) 

1152 >>> mod = 500*np.cos(2*np.pi*0.25*time) 

1153 >>> carrier = amp * np.sin(2*np.pi*3e3*time + mod) 

1154 >>> noise = np.random.normal(scale=np.sqrt(noise_power), 

1155 ... size=time.shape) 

1156 >>> noise *= np.exp(-time/5) 

1157 >>> x = carrier + noise 

1158 

1159 Compute and plot the STFT's magnitude. 

1160 

1161 >>> f, t, Zxx = signal.stft(x, fs, nperseg=1000) 

1162 >>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud') 

1163 >>> plt.title('STFT Magnitude') 

1164 >>> plt.ylabel('Frequency [Hz]') 

1165 >>> plt.xlabel('Time [sec]') 

1166 >>> plt.show() 

1167 """ 

1168 

1169 freqs, time, Zxx = _spectral_helper(x, x, fs, window, nperseg, noverlap, 

1170 nfft, detrend, return_onesided, 

1171 scaling='spectrum', axis=axis, 

1172 mode='stft', boundary=boundary, 

1173 padded=padded) 

1174 

1175 return freqs, time, Zxx 

1176 

1177 

1178def istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, 

1179 input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2): 

1180 r""" 

1181 Perform the inverse Short Time Fourier transform (iSTFT). 

1182 

1183 Parameters 

1184 ---------- 

1185 Zxx : array_like 

1186 STFT of the signal to be reconstructed. If a purely real array 

1187 is passed, it will be cast to a complex data type. 

1188 fs : float, optional 

1189 Sampling frequency of the time series. Defaults to 1.0. 

1190 window : str or tuple or array_like, optional 

1191 Desired window to use. If `window` is a string or tuple, it is 

1192 passed to `get_window` to generate the window values, which are 

1193 DFT-even by default. See `get_window` for a list of windows and 

1194 required parameters. If `window` is array_like it will be used 

1195 directly as the window and its length must be nperseg. Defaults 

1196 to a Hann window. Must match the window used to generate the 

1197 STFT for faithful inversion. 

1198 nperseg : int, optional 

1199 Number of data points corresponding to each STFT segment. This 

1200 parameter must be specified if the number of data points per 

1201 segment is odd, or if the STFT was padded via ``nfft > 

1202 nperseg``. If `None`, the value depends on the shape of 

1203 `Zxx` and `input_onesided`. If `input_onesided` is `True`, 

1204 ``nperseg=2*(Zxx.shape[freq_axis] - 1)``. Otherwise, 

1205 ``nperseg=Zxx.shape[freq_axis]``. Defaults to `None`. 

1206 noverlap : int, optional 

1207 Number of points to overlap between segments. If `None`, half 

1208 of the segment length. Defaults to `None`. When specified, the 

1209 COLA constraint must be met (see Notes below), and should match 

1210 the parameter used to generate the STFT. Defaults to `None`. 

1211 nfft : int, optional 

1212 Number of FFT points corresponding to each STFT segment. This 

1213 parameter must be specified if the STFT was padded via ``nfft > 

1214 nperseg``. If `None`, the default values are the same as for 

1215 `nperseg`, detailed above, with one exception: if 

1216 `input_onesided` is True and 

1217 ``nperseg==2*Zxx.shape[freq_axis] - 1``, `nfft` also takes on 

1218 that value. This case allows the proper inversion of an 

1219 odd-length unpadded STFT using ``nfft=None``. Defaults to 

1220 `None`. 

1221 input_onesided : bool, optional 

1222 If `True`, interpret the input array as one-sided FFTs, such 

1223 as is returned by `stft` with ``return_onesided=True`` and 

1224 `numpy.fft.rfft`. If `False`, interpret the input as a a 

1225 two-sided FFT. Defaults to `True`. 

1226 boundary : bool, optional 

1227 Specifies whether the input signal was extended at its 

1228 boundaries by supplying a non-`None` ``boundary`` argument to 

1229 `stft`. Defaults to `True`. 

1230 time_axis : int, optional 

1231 Where the time segments of the STFT is located; the default is 

1232 the last axis (i.e. ``axis=-1``). 

1233 freq_axis : int, optional 

1234 Where the frequency axis of the STFT is located; the default is 

1235 the penultimate axis (i.e. ``axis=-2``). 

1236 

1237 Returns 

1238 ------- 

1239 t : ndarray 

1240 Array of output data times. 

1241 x : ndarray 

1242 iSTFT of `Zxx`. 

1243 

1244 See Also 

1245 -------- 

1246 stft: Short Time Fourier Transform 

1247 check_COLA: Check whether the Constant OverLap Add (COLA) constraint 

1248 is met 

1249 check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met 

1250 

1251 Notes 

1252 ----- 

1253 In order to enable inversion of an STFT via the inverse STFT with 

1254 `istft`, the signal windowing must obey the constraint of "nonzero 

1255 overlap add" (NOLA): 

1256 

1257 .. math:: \sum_{t}w^{2}[n-tH] \ne 0 

1258 

1259 This ensures that the normalization factors that appear in the denominator 

1260 of the overlap-add reconstruction equation 

1261 

1262 .. math:: x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]} 

1263 

1264 are not zero. The NOLA constraint can be checked with the `check_NOLA` 

1265 function. 

1266 

1267 An STFT which has been modified (via masking or otherwise) is not 

1268 guaranteed to correspond to a exactly realizible signal. This 

1269 function implements the iSTFT via the least-squares estimation 

1270 algorithm detailed in [2]_, which produces a signal that minimizes 

1271 the mean squared error between the STFT of the returned signal and 

1272 the modified STFT. 

1273 

1274 .. versionadded:: 0.19.0 

1275 

1276 References 

1277 ---------- 

1278 .. [1] Oppenheim, Alan V., Ronald W. Schafer, John R. Buck 

1279 "Discrete-Time Signal Processing", Prentice Hall, 1999. 

1280 .. [2] Daniel W. Griffin, Jae S. Lim "Signal Estimation from 

1281 Modified Short-Time Fourier Transform", IEEE 1984, 

1282 10.1109/TASSP.1984.1164317 

1283 

1284 Examples 

1285 -------- 

1286 >>> from scipy import signal 

1287 >>> import matplotlib.pyplot as plt 

1288 

1289 Generate a test signal, a 2 Vrms sine wave at 50Hz corrupted by 

1290 0.001 V**2/Hz of white noise sampled at 1024 Hz. 

1291 

1292 >>> fs = 1024 

1293 >>> N = 10*fs 

1294 >>> nperseg = 512 

1295 >>> amp = 2 * np.sqrt(2) 

1296 >>> noise_power = 0.001 * fs / 2 

1297 >>> time = np.arange(N) / float(fs) 

1298 >>> carrier = amp * np.sin(2*np.pi*50*time) 

1299 >>> noise = np.random.normal(scale=np.sqrt(noise_power), 

1300 ... size=time.shape) 

1301 >>> x = carrier + noise 

1302 

1303 Compute the STFT, and plot its magnitude 

1304 

1305 >>> f, t, Zxx = signal.stft(x, fs=fs, nperseg=nperseg) 

1306 >>> plt.figure() 

1307 >>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud') 

1308 >>> plt.ylim([f[1], f[-1]]) 

1309 >>> plt.title('STFT Magnitude') 

1310 >>> plt.ylabel('Frequency [Hz]') 

1311 >>> plt.xlabel('Time [sec]') 

1312 >>> plt.yscale('log') 

1313 >>> plt.show() 

1314 

1315 Zero the components that are 10% or less of the carrier magnitude, 

1316 then convert back to a time series via inverse STFT 

1317 

1318 >>> Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0) 

1319 >>> _, xrec = signal.istft(Zxx, fs) 

1320 

1321 Compare the cleaned signal with the original and true carrier signals. 

1322 

1323 >>> plt.figure() 

1324 >>> plt.plot(time, x, time, xrec, time, carrier) 

1325 >>> plt.xlim([2, 2.1]) 

1326 >>> plt.xlabel('Time [sec]') 

1327 >>> plt.ylabel('Signal') 

1328 >>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier']) 

1329 >>> plt.show() 

1330 

1331 Note that the cleaned signal does not start as abruptly as the original, 

1332 since some of the coefficients of the transient were also removed: 

1333 

1334 >>> plt.figure() 

1335 >>> plt.plot(time, x, time, xrec, time, carrier) 

1336 >>> plt.xlim([0, 0.1]) 

1337 >>> plt.xlabel('Time [sec]') 

1338 >>> plt.ylabel('Signal') 

1339 >>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier']) 

1340 >>> plt.show() 

1341 

1342 """ 

1343 

1344 # Make sure input is an ndarray of appropriate complex dtype 

1345 Zxx = np.asarray(Zxx) + 0j 

1346 freq_axis = int(freq_axis) 

1347 time_axis = int(time_axis) 

1348 

1349 if Zxx.ndim < 2: 

1350 raise ValueError('Input stft must be at least 2d!') 

1351 

1352 if freq_axis == time_axis: 

1353 raise ValueError('Must specify differing time and frequency axes!') 

1354 

1355 nseg = Zxx.shape[time_axis] 

1356 

1357 if input_onesided: 

1358 # Assume even segment length 

1359 n_default = 2*(Zxx.shape[freq_axis] - 1) 

1360 else: 

1361 n_default = Zxx.shape[freq_axis] 

1362 

1363 # Check windowing parameters 

1364 if nperseg is None: 

1365 nperseg = n_default 

1366 else: 

1367 nperseg = int(nperseg) 

1368 if nperseg < 1: 

1369 raise ValueError('nperseg must be a positive integer') 

1370 

1371 if nfft is None: 

1372 if (input_onesided) and (nperseg == n_default + 1): 

1373 # Odd nperseg, no FFT padding 

1374 nfft = nperseg 

1375 else: 

1376 nfft = n_default 

1377 elif nfft < nperseg: 

1378 raise ValueError('nfft must be greater than or equal to nperseg.') 

1379 else: 

1380 nfft = int(nfft) 

1381 

1382 if noverlap is None: 

1383 noverlap = nperseg//2 

1384 else: 

1385 noverlap = int(noverlap) 

1386 if noverlap >= nperseg: 

1387 raise ValueError('noverlap must be less than nperseg.') 

1388 nstep = nperseg - noverlap 

1389 

1390 # Rearrange axes if necessary 

1391 if time_axis != Zxx.ndim-1 or freq_axis != Zxx.ndim-2: 

1392 # Turn negative indices to positive for the call to transpose 

1393 if freq_axis < 0: 

1394 freq_axis = Zxx.ndim + freq_axis 

1395 if time_axis < 0: 

1396 time_axis = Zxx.ndim + time_axis 

1397 zouter = list(range(Zxx.ndim)) 

1398 for ax in sorted([time_axis, freq_axis], reverse=True): 

1399 zouter.pop(ax) 

1400 Zxx = np.transpose(Zxx, zouter+[freq_axis, time_axis]) 

1401 

1402 # Get window as array 

1403 if isinstance(window, str) or type(window) is tuple: 

1404 win = get_window(window, nperseg) 

1405 else: 

1406 win = np.asarray(window) 

1407 if len(win.shape) != 1: 

1408 raise ValueError('window must be 1-D') 

1409 if win.shape[0] != nperseg: 

1410 raise ValueError('window must have length of {0}'.format(nperseg)) 

1411 

1412 ifunc = sp_fft.irfft if input_onesided else sp_fft.ifft 

1413 xsubs = ifunc(Zxx, axis=-2, n=nfft)[..., :nperseg, :] 

1414 

1415 # Initialize output and normalization arrays 

1416 outputlength = nperseg + (nseg-1)*nstep 

1417 x = np.zeros(list(Zxx.shape[:-2])+[outputlength], dtype=xsubs.dtype) 

1418 norm = np.zeros(outputlength, dtype=xsubs.dtype) 

1419 

1420 if np.result_type(win, xsubs) != xsubs.dtype: 

1421 win = win.astype(xsubs.dtype) 

1422 

1423 xsubs *= win.sum() # This takes care of the 'spectrum' scaling 

1424 

1425 # Construct the output from the ifft segments 

1426 # This loop could perhaps be vectorized/strided somehow... 

1427 for ii in range(nseg): 

1428 # Window the ifft 

1429 x[..., ii*nstep:ii*nstep+nperseg] += xsubs[..., ii] * win 

1430 norm[..., ii*nstep:ii*nstep+nperseg] += win**2 

1431 

1432 # Remove extension points 

1433 if boundary: 

1434 x = x[..., nperseg//2:-(nperseg//2)] 

1435 norm = norm[..., nperseg//2:-(nperseg//2)] 

1436 

1437 # Divide out normalization where non-tiny 

1438 if np.sum(norm > 1e-10) != len(norm): 

1439 warnings.warn("NOLA condition failed, STFT may not be invertible") 

1440 x /= np.where(norm > 1e-10, norm, 1.0) 

1441 

1442 if input_onesided: 

1443 x = x.real 

1444 

1445 # Put axes back 

1446 if x.ndim > 1: 

1447 if time_axis != Zxx.ndim-1: 

1448 if freq_axis < time_axis: 

1449 time_axis -= 1 

1450 x = np.rollaxis(x, -1, time_axis) 

1451 

1452 time = np.arange(x.shape[0])/float(fs) 

1453 return time, x 

1454 

1455 

1456def coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, 

1457 nfft=None, detrend='constant', axis=-1): 

1458 r""" 

1459 Estimate the magnitude squared coherence estimate, Cxy, of 

1460 discrete-time signals X and Y using Welch's method. 

1461 

1462 ``Cxy = abs(Pxy)**2/(Pxx*Pyy)``, where `Pxx` and `Pyy` are power 

1463 spectral density estimates of X and Y, and `Pxy` is the cross 

1464 spectral density estimate of X and Y. 

1465 

1466 Parameters 

1467 ---------- 

1468 x : array_like 

1469 Time series of measurement values 

1470 y : array_like 

1471 Time series of measurement values 

1472 fs : float, optional 

1473 Sampling frequency of the `x` and `y` time series. Defaults 

1474 to 1.0. 

1475 window : str or tuple or array_like, optional 

1476 Desired window to use. If `window` is a string or tuple, it is 

1477 passed to `get_window` to generate the window values, which are 

1478 DFT-even by default. See `get_window` for a list of windows and 

1479 required parameters. If `window` is array_like it will be used 

1480 directly as the window and its length must be nperseg. Defaults 

1481 to a Hann window. 

1482 nperseg : int, optional 

1483 Length of each segment. Defaults to None, but if window is str or 

1484 tuple, is set to 256, and if window is array_like, is set to the 

1485 length of the window. 

1486 noverlap: int, optional 

1487 Number of points to overlap between segments. If `None`, 

1488 ``noverlap = nperseg // 2``. Defaults to `None`. 

1489 nfft : int, optional 

1490 Length of the FFT used, if a zero padded FFT is desired. If 

1491 `None`, the FFT length is `nperseg`. Defaults to `None`. 

1492 detrend : str or function or `False`, optional 

1493 Specifies how to detrend each segment. If `detrend` is a 

1494 string, it is passed as the `type` argument to the `detrend` 

1495 function. If it is a function, it takes a segment and returns a 

1496 detrended segment. If `detrend` is `False`, no detrending is 

1497 done. Defaults to 'constant'. 

1498 axis : int, optional 

1499 Axis along which the coherence is computed for both inputs; the 

1500 default is over the last axis (i.e. ``axis=-1``). 

1501 

1502 Returns 

1503 ------- 

1504 f : ndarray 

1505 Array of sample frequencies. 

1506 Cxy : ndarray 

1507 Magnitude squared coherence of x and y. 

1508 

1509 See Also 

1510 -------- 

1511 periodogram: Simple, optionally modified periodogram 

1512 lombscargle: Lomb-Scargle periodogram for unevenly sampled data 

1513 welch: Power spectral density by Welch's method. 

1514 csd: Cross spectral density by Welch's method. 

1515 

1516 Notes 

1517 -------- 

1518 An appropriate amount of overlap will depend on the choice of window 

1519 and on your requirements. For the default Hann window an overlap of 

1520 50% is a reasonable trade off between accurately estimating the 

1521 signal power, while not over counting any of the data. Narrower 

1522 windows may require a larger overlap. 

1523 

1524 .. versionadded:: 0.16.0 

1525 

1526 References 

1527 ---------- 

1528 .. [1] P. Welch, "The use of the fast Fourier transform for the 

1529 estimation of power spectra: A method based on time averaging 

1530 over short, modified periodograms", IEEE Trans. Audio 

1531 Electroacoust. vol. 15, pp. 70-73, 1967. 

1532 .. [2] Stoica, Petre, and Randolph Moses, "Spectral Analysis of 

1533 Signals" Prentice Hall, 2005 

1534 

1535 Examples 

1536 -------- 

1537 >>> from scipy import signal 

1538 >>> import matplotlib.pyplot as plt 

1539 

1540 Generate two test signals with some common features. 

1541 

1542 >>> fs = 10e3 

1543 >>> N = 1e5 

1544 >>> amp = 20 

1545 >>> freq = 1234.0 

1546 >>> noise_power = 0.001 * fs / 2 

1547 >>> time = np.arange(N) / fs 

1548 >>> b, a = signal.butter(2, 0.25, 'low') 

1549 >>> x = np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

1550 >>> y = signal.lfilter(b, a, x) 

1551 >>> x += amp*np.sin(2*np.pi*freq*time) 

1552 >>> y += np.random.normal(scale=0.1*np.sqrt(noise_power), size=time.shape) 

1553 

1554 Compute and plot the coherence. 

1555 

1556 >>> f, Cxy = signal.coherence(x, y, fs, nperseg=1024) 

1557 >>> plt.semilogy(f, Cxy) 

1558 >>> plt.xlabel('frequency [Hz]') 

1559 >>> plt.ylabel('Coherence') 

1560 >>> plt.show() 

1561 """ 

1562 

1563 freqs, Pxx = welch(x, fs=fs, window=window, nperseg=nperseg, 

1564 noverlap=noverlap, nfft=nfft, detrend=detrend, 

1565 axis=axis) 

1566 _, Pyy = welch(y, fs=fs, window=window, nperseg=nperseg, noverlap=noverlap, 

1567 nfft=nfft, detrend=detrend, axis=axis) 

1568 _, Pxy = csd(x, y, fs=fs, window=window, nperseg=nperseg, 

1569 noverlap=noverlap, nfft=nfft, detrend=detrend, axis=axis) 

1570 

1571 Cxy = np.abs(Pxy)**2 / Pxx / Pyy 

1572 

1573 return freqs, Cxy 

1574 

1575 

1576def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, 

1577 nfft=None, detrend='constant', return_onesided=True, 

1578 scaling='density', axis=-1, mode='psd', boundary=None, 

1579 padded=False): 

1580 """ 

1581 Calculate various forms of windowed FFTs for PSD, CSD, etc. 

1582 

1583 This is a helper function that implements the commonality between 

1584 the stft, psd, csd, and spectrogram functions. It is not designed to 

1585 be called externally. The windows are not averaged over; the result 

1586 from each window is returned. 

1587 

1588 Parameters 

1589 --------- 

1590 x : array_like 

1591 Array or sequence containing the data to be analyzed. 

1592 y : array_like 

1593 Array or sequence containing the data to be analyzed. If this is 

1594 the same object in memory as `x` (i.e. ``_spectral_helper(x, 

1595 x, ...)``), the extra computations are spared. 

1596 fs : float, optional 

1597 Sampling frequency of the time series. Defaults to 1.0. 

1598 window : str or tuple or array_like, optional 

1599 Desired window to use. If `window` is a string or tuple, it is 

1600 passed to `get_window` to generate the window values, which are 

1601 DFT-even by default. See `get_window` for a list of windows and 

1602 required parameters. If `window` is array_like it will be used 

1603 directly as the window and its length must be nperseg. Defaults 

1604 to a Hann window. 

1605 nperseg : int, optional 

1606 Length of each segment. Defaults to None, but if window is str or 

1607 tuple, is set to 256, and if window is array_like, is set to the 

1608 length of the window. 

1609 noverlap : int, optional 

1610 Number of points to overlap between segments. If `None`, 

1611 ``noverlap = nperseg // 2``. Defaults to `None`. 

1612 nfft : int, optional 

1613 Length of the FFT used, if a zero padded FFT is desired. If 

1614 `None`, the FFT length is `nperseg`. Defaults to `None`. 

1615 detrend : str or function or `False`, optional 

1616 Specifies how to detrend each segment. If `detrend` is a 

1617 string, it is passed as the `type` argument to the `detrend` 

1618 function. If it is a function, it takes a segment and returns a 

1619 detrended segment. If `detrend` is `False`, no detrending is 

1620 done. Defaults to 'constant'. 

1621 return_onesided : bool, optional 

1622 If `True`, return a one-sided spectrum for real data. If 

1623 `False` return a two-sided spectrum. Defaults to `True`, but for 

1624 complex data, a two-sided spectrum is always returned. 

1625 scaling : { 'density', 'spectrum' }, optional 

1626 Selects between computing the cross spectral density ('density') 

1627 where `Pxy` has units of V**2/Hz and computing the cross 

1628 spectrum ('spectrum') where `Pxy` has units of V**2, if `x` 

1629 and `y` are measured in V and `fs` is measured in Hz. 

1630 Defaults to 'density' 

1631 axis : int, optional 

1632 Axis along which the FFTs are computed; the default is over the 

1633 last axis (i.e. ``axis=-1``). 

1634 mode: str {'psd', 'stft'}, optional 

1635 Defines what kind of return values are expected. Defaults to 

1636 'psd'. 

1637 boundary : str or None, optional 

1638 Specifies whether the input signal is extended at both ends, and 

1639 how to generate the new values, in order to center the first 

1640 windowed segment on the first input point. This has the benefit 

1641 of enabling reconstruction of the first input point when the 

1642 employed window function starts at zero. Valid options are 

1643 ``['even', 'odd', 'constant', 'zeros', None]``. Defaults to 

1644 `None`. 

1645 padded : bool, optional 

1646 Specifies whether the input signal is zero-padded at the end to 

1647 make the signal fit exactly into an integer number of window 

1648 segments, so that all of the signal is included in the output. 

1649 Defaults to `False`. Padding occurs after boundary extension, if 

1650 `boundary` is not `None`, and `padded` is `True`. 

1651 Returns 

1652 ------- 

1653 freqs : ndarray 

1654 Array of sample frequencies. 

1655 t : ndarray 

1656 Array of times corresponding to each data segment 

1657 result : ndarray 

1658 Array of output data, contents dependent on *mode* kwarg. 

1659 

1660 Notes 

1661 ----- 

1662 Adapted from matplotlib.mlab 

1663 

1664 .. versionadded:: 0.16.0 

1665 """ 

1666 if mode not in ['psd', 'stft']: 

1667 raise ValueError("Unknown value for mode %s, must be one of: " 

1668 "{'psd', 'stft'}" % mode) 

1669 

1670 boundary_funcs = {'even': even_ext, 

1671 'odd': odd_ext, 

1672 'constant': const_ext, 

1673 'zeros': zero_ext, 

1674 None: None} 

1675 

1676 if boundary not in boundary_funcs: 

1677 raise ValueError("Unknown boundary option '{0}', must be one of: {1}" 

1678 .format(boundary, list(boundary_funcs.keys()))) 

1679 

1680 # If x and y are the same object we can save ourselves some computation. 

1681 same_data = y is x 

1682 

1683 if not same_data and mode != 'psd': 

1684 raise ValueError("x and y must be equal if mode is 'stft'") 

1685 

1686 axis = int(axis) 

1687 

1688 # Ensure we have np.arrays, get outdtype 

1689 x = np.asarray(x) 

1690 if not same_data: 

1691 y = np.asarray(y) 

1692 outdtype = np.result_type(x, y, np.complex64) 

1693 else: 

1694 outdtype = np.result_type(x, np.complex64) 

1695 

1696 if not same_data: 

1697 # Check if we can broadcast the outer axes together 

1698 xouter = list(x.shape) 

1699 youter = list(y.shape) 

1700 xouter.pop(axis) 

1701 youter.pop(axis) 

1702 try: 

1703 outershape = np.broadcast(np.empty(xouter), np.empty(youter)).shape 

1704 except ValueError: 

1705 raise ValueError('x and y cannot be broadcast together.') 

1706 

1707 if same_data: 

1708 if x.size == 0: 

1709 return np.empty(x.shape), np.empty(x.shape), np.empty(x.shape) 

1710 else: 

1711 if x.size == 0 or y.size == 0: 

1712 outshape = outershape + (min([x.shape[axis], y.shape[axis]]),) 

1713 emptyout = np.rollaxis(np.empty(outshape), -1, axis) 

1714 return emptyout, emptyout, emptyout 

1715 

1716 if x.ndim > 1: 

1717 if axis != -1: 

1718 x = np.rollaxis(x, axis, len(x.shape)) 

1719 if not same_data and y.ndim > 1: 

1720 y = np.rollaxis(y, axis, len(y.shape)) 

1721 

1722 # Check if x and y are the same length, zero-pad if necessary 

1723 if not same_data: 

1724 if x.shape[-1] != y.shape[-1]: 

1725 if x.shape[-1] < y.shape[-1]: 

1726 pad_shape = list(x.shape) 

1727 pad_shape[-1] = y.shape[-1] - x.shape[-1] 

1728 x = np.concatenate((x, np.zeros(pad_shape)), -1) 

1729 else: 

1730 pad_shape = list(y.shape) 

1731 pad_shape[-1] = x.shape[-1] - y.shape[-1] 

1732 y = np.concatenate((y, np.zeros(pad_shape)), -1) 

1733 

1734 if nperseg is not None: # if specified by user 

1735 nperseg = int(nperseg) 

1736 if nperseg < 1: 

1737 raise ValueError('nperseg must be a positive integer') 

1738 

1739 # parse window; if array like, then set nperseg = win.shape 

1740 win, nperseg = _triage_segments(window, nperseg, input_length=x.shape[-1]) 

1741 

1742 if nfft is None: 

1743 nfft = nperseg 

1744 elif nfft < nperseg: 

1745 raise ValueError('nfft must be greater than or equal to nperseg.') 

1746 else: 

1747 nfft = int(nfft) 

1748 

1749 if noverlap is None: 

1750 noverlap = nperseg//2 

1751 else: 

1752 noverlap = int(noverlap) 

1753 if noverlap >= nperseg: 

1754 raise ValueError('noverlap must be less than nperseg.') 

1755 nstep = nperseg - noverlap 

1756 

1757 # Padding occurs after boundary extension, so that the extended signal ends 

1758 # in zeros, instead of introducing an impulse at the end. 

1759 # I.e. if x = [..., 3, 2] 

1760 # extend then pad -> [..., 3, 2, 2, 3, 0, 0, 0] 

1761 # pad then extend -> [..., 3, 2, 0, 0, 0, 2, 3] 

1762 

1763 if boundary is not None: 

1764 ext_func = boundary_funcs[boundary] 

1765 x = ext_func(x, nperseg//2, axis=-1) 

1766 if not same_data: 

1767 y = ext_func(y, nperseg//2, axis=-1) 

1768 

1769 if padded: 

1770 # Pad to integer number of windowed segments 

1771 # I.e make x.shape[-1] = nperseg + (nseg-1)*nstep, with integer nseg 

1772 nadd = (-(x.shape[-1]-nperseg) % nstep) % nperseg 

1773 zeros_shape = list(x.shape[:-1]) + [nadd] 

1774 x = np.concatenate((x, np.zeros(zeros_shape)), axis=-1) 

1775 if not same_data: 

1776 zeros_shape = list(y.shape[:-1]) + [nadd] 

1777 y = np.concatenate((y, np.zeros(zeros_shape)), axis=-1) 

1778 

1779 # Handle detrending and window functions 

1780 if not detrend: 

1781 def detrend_func(d): 

1782 return d 

1783 elif not hasattr(detrend, '__call__'): 

1784 def detrend_func(d): 

1785 return signaltools.detrend(d, type=detrend, axis=-1) 

1786 elif axis != -1: 

1787 # Wrap this function so that it receives a shape that it could 

1788 # reasonably expect to receive. 

1789 def detrend_func(d): 

1790 d = np.rollaxis(d, -1, axis) 

1791 d = detrend(d) 

1792 return np.rollaxis(d, axis, len(d.shape)) 

1793 else: 

1794 detrend_func = detrend 

1795 

1796 if np.result_type(win, np.complex64) != outdtype: 

1797 win = win.astype(outdtype) 

1798 

1799 if scaling == 'density': 

1800 scale = 1.0 / (fs * (win*win).sum()) 

1801 elif scaling == 'spectrum': 

1802 scale = 1.0 / win.sum()**2 

1803 else: 

1804 raise ValueError('Unknown scaling: %r' % scaling) 

1805 

1806 if mode == 'stft': 

1807 scale = np.sqrt(scale) 

1808 

1809 if return_onesided: 

1810 if np.iscomplexobj(x): 

1811 sides = 'twosided' 

1812 warnings.warn('Input data is complex, switching to ' 

1813 'return_onesided=False') 

1814 else: 

1815 sides = 'onesided' 

1816 if not same_data: 

1817 if np.iscomplexobj(y): 

1818 sides = 'twosided' 

1819 warnings.warn('Input data is complex, switching to ' 

1820 'return_onesided=False') 

1821 else: 

1822 sides = 'twosided' 

1823 

1824 if sides == 'twosided': 

1825 freqs = sp_fft.fftfreq(nfft, 1/fs) 

1826 elif sides == 'onesided': 

1827 freqs = sp_fft.rfftfreq(nfft, 1/fs) 

1828 

1829 # Perform the windowed FFTs 

1830 result = _fft_helper(x, win, detrend_func, nperseg, noverlap, nfft, sides) 

1831 

1832 if not same_data: 

1833 # All the same operations on the y data 

1834 result_y = _fft_helper(y, win, detrend_func, nperseg, noverlap, nfft, 

1835 sides) 

1836 result = np.conjugate(result) * result_y 

1837 elif mode == 'psd': 

1838 result = np.conjugate(result) * result 

1839 

1840 result *= scale 

1841 if sides == 'onesided' and mode == 'psd': 

1842 if nfft % 2: 

1843 result[..., 1:] *= 2 

1844 else: 

1845 # Last point is unpaired Nyquist freq point, don't double 

1846 result[..., 1:-1] *= 2 

1847 

1848 time = np.arange(nperseg/2, x.shape[-1] - nperseg/2 + 1, 

1849 nperseg - noverlap)/float(fs) 

1850 if boundary is not None: 

1851 time -= (nperseg/2) / fs 

1852 

1853 result = result.astype(outdtype) 

1854 

1855 # All imaginary parts are zero anyways 

1856 if same_data and mode != 'stft': 

1857 result = result.real 

1858 

1859 # Output is going to have new last axis for time/window index, so a 

1860 # negative axis index shifts down one 

1861 if axis < 0: 

1862 axis -= 1 

1863 

1864 # Roll frequency axis back to axis where the data came from 

1865 result = np.rollaxis(result, -1, axis) 

1866 

1867 return freqs, time, result 

1868 

1869 

1870def _fft_helper(x, win, detrend_func, nperseg, noverlap, nfft, sides): 

1871 """ 

1872 Calculate windowed FFT, for internal use by 

1873 scipy.signal._spectral_helper 

1874 

1875 This is a helper function that does the main FFT calculation for 

1876 `_spectral helper`. All input validation is performed there, and the 

1877 data axis is assumed to be the last axis of x. It is not designed to 

1878 be called externally. The windows are not averaged over; the result 

1879 from each window is returned. 

1880 

1881 Returns 

1882 ------- 

1883 result : ndarray 

1884 Array of FFT data 

1885 

1886 Notes 

1887 ----- 

1888 Adapted from matplotlib.mlab 

1889 

1890 .. versionadded:: 0.16.0 

1891 """ 

1892 # Created strided array of data segments 

1893 if nperseg == 1 and noverlap == 0: 

1894 result = x[..., np.newaxis] 

1895 else: 

1896 # https://stackoverflow.com/a/5568169 

1897 step = nperseg - noverlap 

1898 shape = x.shape[:-1]+((x.shape[-1]-noverlap)//step, nperseg) 

1899 strides = x.strides[:-1]+(step*x.strides[-1], x.strides[-1]) 

1900 result = np.lib.stride_tricks.as_strided(x, shape=shape, 

1901 strides=strides) 

1902 

1903 # Detrend each data segment individually 

1904 result = detrend_func(result) 

1905 

1906 # Apply window by multiplication 

1907 result = win * result 

1908 

1909 # Perform the fft. Acts on last axis by default. Zero-pads automatically 

1910 if sides == 'twosided': 

1911 func = sp_fft.fft 

1912 else: 

1913 result = result.real 

1914 func = sp_fft.rfft 

1915 result = func(result, n=nfft) 

1916 

1917 return result 

1918 

1919 

1920def _triage_segments(window, nperseg, input_length): 

1921 """ 

1922 Parses window and nperseg arguments for spectrogram and _spectral_helper. 

1923 This is a helper function, not meant to be called externally. 

1924 

1925 Parameters 

1926 ---------- 

1927 window : string, tuple, or ndarray 

1928 If window is specified by a string or tuple and nperseg is not 

1929 specified, nperseg is set to the default of 256 and returns a window of 

1930 that length. 

1931 If instead the window is array_like and nperseg is not specified, then 

1932 nperseg is set to the length of the window. A ValueError is raised if 

1933 the user supplies both an array_like window and a value for nperseg but 

1934 nperseg does not equal the length of the window. 

1935 

1936 nperseg : int 

1937 Length of each segment 

1938 

1939 input_length: int 

1940 Length of input signal, i.e. x.shape[-1]. Used to test for errors. 

1941 

1942 Returns 

1943 ------- 

1944 win : ndarray 

1945 window. If function was called with string or tuple than this will hold 

1946 the actual array used as a window. 

1947 

1948 nperseg : int 

1949 Length of each segment. If window is str or tuple, nperseg is set to 

1950 256. If window is array_like, nperseg is set to the length of the 

1951 6 

1952 window. 

1953 """ 

1954 

1955 # parse window; if array like, then set nperseg = win.shape 

1956 if isinstance(window, str) or isinstance(window, tuple): 

1957 # if nperseg not specified 

1958 if nperseg is None: 

1959 nperseg = 256 # then change to default 

1960 if nperseg > input_length: 

1961 warnings.warn('nperseg = {0:d} is greater than input length ' 

1962 ' = {1:d}, using nperseg = {1:d}' 

1963 .format(nperseg, input_length)) 

1964 nperseg = input_length 

1965 win = get_window(window, nperseg) 

1966 else: 

1967 win = np.asarray(window) 

1968 if len(win.shape) != 1: 

1969 raise ValueError('window must be 1-D') 

1970 if input_length < win.shape[-1]: 

1971 raise ValueError('window is longer than input signal') 

1972 if nperseg is None: 

1973 nperseg = win.shape[0] 

1974 elif nperseg is not None: 

1975 if nperseg != win.shape[0]: 

1976 raise ValueError("value specified for nperseg is different" 

1977 " from length of window") 

1978 return win, nperseg 

1979 

1980 

1981def _median_bias(n): 

1982 """ 

1983 Returns the bias of the median of a set of periodograms relative to 

1984 the mean. 

1985 

1986 See arXiv:gr-qc/0509116 Appendix B for details. 

1987 

1988 Parameters 

1989 ---------- 

1990 n : int 

1991 Numbers of periodograms being averaged. 

1992 

1993 Returns 

1994 ------- 

1995 bias : float 

1996 Calculated bias. 

1997 """ 

1998 ii_2 = 2 * np.arange(1., (n-1) // 2 + 1) 

1999 return 1 + np.sum(1. / (ii_2 + 1) - 1. / ii_2)