From 5dc7e763b4058e24e24bec885e93e2da120edcd7 Mon Sep 17 00:00:00 2001 From: NickleDave Date: Thu, 1 Dec 2016 12:04:31 -0500 Subject: [PATCH] BUG: make _spectral_helper use window size if nperseg is unspecified This applies to cases where the window is specified by an array of values. This behavior was advertised in docstrings, but not actually implemented. --- THANKS.txt | 1 + scipy/signal/spectral.py | 122 +++++++++++++++++++++------- scipy/signal/tests/test_spectral.py | 96 ++++++++++++++++++---- 3 files changed, 176 insertions(+), 43 deletions(-) diff --git a/THANKS.txt b/THANKS.txt index e8e666964cb6..59786e080e21 100644 --- a/THANKS.txt +++ b/THANKS.txt @@ -169,6 +169,7 @@ John Draper for bug fixes. Alvaro Sanchez-Gonzalez for axis-dependent modes in multidimensional filters. Alessandro Pietro Bardelli for improvements to pdist/cdist and to related tests. Jonathan T. Siebert for bug fixes. +David Nicholson for bug fixes in spectral functions. Institutions ------------ diff --git a/scipy/signal/spectral.py b/scipy/signal/spectral.py index 6a9019322833..50e520b5ea99 100644 --- a/scipy/signal/spectral.py +++ b/scipy/signal/spectral.py @@ -141,7 +141,7 @@ def periodogram(x, fs=1.0, window=None, nfft=None, detrend='constant', scaling, axis) -def welch(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, +def welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1): """ Estimate power spectral density using Welch's method. @@ -162,7 +162,9 @@ def welch(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, directly as the window and its length will be used for nperseg. Defaults to 'hann'. nperseg : int, optional - Length of each segment. Defaults to 256. + Length of each segment. Defaults to None, but if window is str or + tuple, is set to 256, and if window is array_like, is set to the + length of the window. noverlap : int, optional Number of points to overlap between segments. If None, ``noverlap = nperseg // 2``. Defaults to None. @@ -275,7 +277,7 @@ def welch(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, return freqs, Pxx.real -def csd(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, +def csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1): """ Estimate the cross power spectral density, Pxy, using Welch's method. @@ -294,12 +296,14 @@ def csd(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, directly as the window and its length will be used for nperseg. Defaults to 'hann'. nperseg : int, optional - Length of each segment. Defaults to 256. + Length of each segment. Defaults to None, but if window is str or + tuple, is set to 256, and if window is array_like, is set to the + length of the window. noverlap: int, optional Number of points to overlap between segments. If None, ``noverlap = nperseg // 2``. Defaults to None. nfft : int, optional - Length of the FFT used, if a zero padded FFT is desired. If None, + Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is `nperseg`. Defaults to None. detrend : str or function or False, optional Specifies how to detrend each segment. If `detrend` is a string, @@ -400,7 +404,7 @@ def csd(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, return freqs, Pxy -def spectrogram(x, fs=1.0, window=('tukey',.25), nperseg=256, noverlap=None, +def spectrogram(x, fs=1.0, window=('tukey',.25), nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, mode='psd'): """ @@ -421,7 +425,9 @@ def spectrogram(x, fs=1.0, window=('tukey',.25), nperseg=256, noverlap=None, directly as the window and its length will be used for nperseg. Defaults to a Tukey window with shape parameter of 0.25. nperseg : int, optional - Length of each segment. Defaults to 256. + Length of each segment. Defaults to None, but if window is str or + tuple, is set to 256, and if window is array_like, is set to the + length of the window. noverlap : int, optional Number of points to overlap between segments. If None, ``noverlap = nperseg // 8``. Defaults to None. @@ -507,6 +513,11 @@ def spectrogram(x, fs=1.0, window=('tukey',.25), nperseg=256, noverlap=None, >>> plt.xlabel('Time [sec]') >>> plt.show() """ + + # need to set default for nperseg before setting default for noverlap below + window, nperseg = _triage_segments(window, nperseg, + input_length=x.shape[-1]) + # Less overlap than welch, so samples are more statisically independent if noverlap is None: noverlap = nperseg // 8 @@ -518,7 +529,7 @@ def spectrogram(x, fs=1.0, window=('tukey',.25), nperseg=256, noverlap=None, return freqs, time, Pxy -def coherence(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, +def coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', axis=-1): """ Estimate the magnitude squared coherence estimate, Cxy, of discrete-time @@ -542,7 +553,9 @@ def coherence(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, directly as the window and its length will be used for nperseg. Defaults to 'hann'. nperseg : int, optional - Length of each segment. Defaults to 256. + Length of each segment. Defaults to None, but if window is str or + tuple, is set to 256, and if window is array_like, is set to the + length of the window. noverlap: int, optional Number of points to overlap between segments. If None, ``noverlap = nperseg // 2``. Defaults to None. @@ -629,7 +642,7 @@ def coherence(x, y, fs=1.0, window='hann', nperseg=256, noverlap=None, return freqs, Cxy -def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=256, +def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='spectrum', axis=-1, mode='psd'): @@ -647,7 +660,7 @@ def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=256, Array or sequence containing the data to be analyzed. y : array_like Array or sequence containing the data to be analyzed. If this is - the same object in memoery as x (i.e. _spectral_helper(x, x, ...)), + the same object in memory as x (i.e. _spectral_helper(x, x, ...)), the extra computations are spared. fs : float, optional Sampling frequency of the time series. Defaults to 1.0. @@ -657,7 +670,9 @@ def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=256, directly as the window and its length will be used for nperseg. Defaults to 'hann'. nperseg : int, optional - Length of each segment. Defaults to 256. + Length of each segment. Defaults to None, but if window is str or + tuple, is set to 256, and if window is array_like, is set to the + length of the window. noverlap : int, optional Number of points to overlap between segments. If None, ``noverlap = nperseg // 2``. Defaults to None. @@ -766,15 +781,13 @@ def _spectral_helper(x, y, fs=1.0, window='hann', nperseg=256, pad_shape[-1] = x.shape[-1] - y.shape[-1] y = np.concatenate((y, np.zeros(pad_shape)), -1) - # X and Y are same length now, can test nperseg with either - if x.shape[-1] < nperseg: - warnings.warn('nperseg = {0:d}, is greater than input length = {1:d}, ' - 'using nperseg = {1:d}'.format(nperseg, x.shape[-1])) - nperseg = x.shape[-1] + if nperseg is not None: # if specified by user + nperseg = int(nperseg) + if nperseg < 1: + raise ValueError('nperseg must be a positive integer') - nperseg = int(nperseg) - if nperseg < 1: - raise ValueError('nperseg must be a positive integer') + # parse window; if array like, then set nperseg = win.shape + win, nperseg = _triage_segments(window, nperseg,input_length=x.shape[-1]) if nfft is None: nfft = nperseg @@ -807,15 +820,6 @@ def detrend_func(d): else: detrend_func = detrend - if isinstance(window, string_types) or type(window) is tuple: - win = get_window(window, nperseg) - else: - win = np.asarray(window) - if len(win.shape) != 1: - raise ValueError('window must be 1-D') - if win.shape[0] != nperseg: - raise ValueError('window must have length of nperseg') - if np.result_type(win,np.complex64) != outdtype: win = win.astype(outdtype) @@ -952,3 +956,63 @@ def _fft_helper(x, win, detrend_func, nperseg, noverlap, nfft): result = fftpack.fft(result, n=nfft) return result + +def _triage_segments(window, nperseg,input_length): + """ + Parses window and nperseg arguments for spectrogram and _spectral_helper. + This is a helper function, not meant to be called externally. + + Parameters + --------- + window : string, tuple, or ndarray + If window is specified by a string or tuple and nperseg is not + specified, nperseg is set to the default of 256 and returns a window of + that length. + If instead the window is array_like and nperseg is not specified, then + nperseg is set to the length of the window. A ValueError is raised if + the user supplies both an array_like window and a value for nperseg but + nperseg does not equal the length of the window. + + nperseg : int + Length of each segment + + input_length: int + Length of input signal, i.e. x.shape[-1]. Used to test for errors. + + Returns + ------- + win : ndarray + window. If function was called with string or tuple than this will hold + the actual array used as a window. + + nperseg : int + Length of each segment. If window is str or tuple, nperseg is set to + 256. If window is array_like, nperseg is set to the length of the + 6 + window. + """ + + #parse window; if array like, then set nperseg = win.shape + if isinstance(window, string_types) or isinstance(window, tuple): + # if nperseg not specified + if nperseg is None: + nperseg = 256 # then change to default + if nperseg > input_length: + warnings.warn('nperseg = {0:d} is greater than input length ' + ' = {1:d}, using nperseg = {1:d}' + .format(nperseg, input_length)) + nperseg = input_length + win = get_window(window, nperseg) + else: + win = np.asarray(window) + if len(win.shape) != 1: + raise ValueError('window must be 1-D') + if input_length < win.shape[-1]: + raise ValueError('window is longer than input signal') + if nperseg is None: + nperseg = win.shape[0] + elif nperseg is not None: + if nperseg != win.shape[0]: + raise ValueError("value specified for nperseg is different from" + " length of window") + return win, nperseg diff --git a/scipy/signal/tests/test_spectral.py b/scipy/signal/tests/test_spectral.py index a5035a9a8dea..916900675271 100644 --- a/scipy/signal/tests/test_spectral.py +++ b/scipy/signal/tests/test_spectral.py @@ -122,6 +122,9 @@ def test_window_external(self): fe, pe = periodogram(x, 10, win) assert_array_almost_equal_nulp(p, pe) assert_array_almost_equal_nulp(f, fe) + win_err = signal.get_window('hann', 32) + assert_raises(ValueError, periodogram, x, + 10, win_err) # win longer than signal def test_padded_fft(self): x = np.zeros(16) @@ -353,11 +356,18 @@ def test_window_external(self): x = np.zeros(16) x[0] = 1 x[8] = 1 - f, p = welch(x, 10, 'hann', 8) + f, p = welch(x, 10, 'hann', nperseg=8) win = signal.get_window('hann', 8) - fe, pe = welch(x, 10, win, 8) + fe, pe = welch(x, 10, win, nperseg=None) assert_array_almost_equal_nulp(p, pe) assert_array_almost_equal_nulp(f, fe) + assert_array_equal(fe.shape, (5,)) # because win length used as nperseg + assert_array_equal(pe.shape, (5,)) + assert_raises(ValueError, welch, x, + 10, win, nperseg=4) # because nperseg != win.shape[-1] + win_err = signal.get_window('hann', 32) + assert_raises(ValueError, welch, x, + 10, win_err, nperseg=None) # win longer than signal def test_empty_input(self): f, p = welch([]) @@ -377,13 +387,18 @@ def test_empty_input_other_axis(self): def test_short_data(self): x = np.zeros(8) x[0] = 1 + #for string-like window, input signal length < nperseg value gives + #UserWarning, sets nperseg to x.shape[-1] with warnings.catch_warnings(): warnings.simplefilter('ignore', UserWarning) - f, p = welch(x) - - f1, p1 = welch(x, nperseg=8) - assert_allclose(f, f1) - assert_allclose(p, p1) + f, p = welch(x,window='hann') # default nperseg + f1, p1 = welch(x,window='hann', + nperseg=256) # user-specified nperseg + f2, p2 = welch(x, nperseg=8) # valid nperseg, doesn't give warning + assert_allclose(f, f2) + assert_allclose(p, p2) + assert_allclose(f1, f2) + assert_allclose(p1, p2) def test_window_long_or_nd(self): with warnings.catch_warnings(): @@ -474,7 +489,7 @@ def test_padded_freqs(self): def test_window_correction(self): A = 20 fs = 1e4 - nperseg = fs//10 + nperseg = int(fs//10) fsig = 300 ii = int(fsig*nperseg//fs) # Freq index of fsig @@ -493,6 +508,7 @@ def test_window_correction(self): assert_allclose(np.sqrt(np.trapz(p_dens, freq)), A*np.sqrt(2)/2, rtol=1e-3) + class TestCSD: def test_pad_shorter_x(self): x = np.zeros(8) @@ -657,9 +673,16 @@ def test_window_external(self): x[8] = 1 f, p = csd(x, x, 10, 'hann', 8) win = signal.get_window('hann', 8) - fe, pe = csd(x, x, 10, win, 8) + fe, pe = csd(x, x, 10, win, nperseg=None) assert_array_almost_equal_nulp(p, pe) assert_array_almost_equal_nulp(f, fe) + assert_array_equal(fe.shape, (5,)) # because win length used as nperseg + assert_array_equal(pe.shape, (5,)) + assert_raises(ValueError, csd, x, x, + 10, win, nperseg=256) # because nperseg != win.shape[-1] + win_err = signal.get_window('hann', 32) + assert_raises(ValueError, csd, x, x, + 10, win_err, nperseg=None) # because win longer than signal def test_empty_input(self): f, p = csd([],np.zeros(10)) @@ -700,13 +723,19 @@ def test_empty_input_other_axis(self): def test_short_data(self): x = np.zeros(8) x[0] = 1 + + #for string-like window, input signal length < nperseg value gives + #UserWarning, sets nperseg to x.shape[-1] with warnings.catch_warnings(): warnings.simplefilter('ignore', UserWarning) - f, p = csd(x, x) - - f1, p1 = csd(x, x, nperseg=8) - assert_allclose(f, f1) - assert_allclose(p, p1) + f, p = csd(x, x, window='hann') # default nperseg + f1, p1 = csd(x, x, window='hann', + nperseg=256) # user-specified nperseg + f2, p2 = csd(x, x, nperseg=8) # valid nperseg, doesn't give warning + assert_allclose(f, f2) + assert_allclose(p, p2) + assert_allclose(f1, f2) + assert_allclose(p1, p2) def test_window_long_or_nd(self): with warnings.catch_warnings(): @@ -797,6 +826,7 @@ def test_padded_freqs(self): assert_allclose(f, fodd) assert_allclose(f, feven) + class TestCoherence: def test_identical_input(self): x = np.random.randn(20) @@ -835,6 +865,44 @@ def test_average_all_segments(self): assert_allclose(f, fw) assert_allclose(np.mean(P, axis=-1), Pw) + def test_window_external(self): + x = np.random.randn(1024) + + fs = 1.0 + window = ('tukey', 0.25) + nperseg = 16 + noverlap = 2 + f, _, P = spectrogram(x, fs, window, nperseg, noverlap) + + win = signal.get_window(('tukey', 0.25), 16) + fe, _, Pe = spectrogram(x, fs, win, nperseg=None, noverlap=2) + assert_array_equal(fe.shape, (9,)) # because win length used as nperseg + assert_array_equal(Pe.shape, (9,73)) + assert_raises(ValueError, spectrogram, x, + fs, win, nperseg=8) # because nperseg != win.shape[-1] + win_err = signal.get_window(('tukey', 0.25), 2048) + assert_raises(ValueError, spectrogram, x, + fs, win_err, nperseg=None) # win longer than signal + + def test_short_data(self): + x = np.random.randn(1024) + fs = 1.0 + + #for string-like window, input signal length < nperseg value gives + #UserWarning, sets nperseg to x.shape[-1] + with warnings.catch_warnings(): + warnings.simplefilter('ignore', UserWarning) + f, _, p = spectrogram(x, fs, + window=('tukey',0.25)) # default nperseg + f1, _, p1 = spectrogram(x, fs, window=('tukey',0.25), + nperseg=1025) # user-specified nperseg + f2, _, p2 = spectrogram(x, fs, nperseg=256) # to compare w/default + f3, _, p3 = spectrogram(x, fs, nperseg=1024) # compare w/user-spec'd + assert_allclose(f, f2) + assert_allclose(p, p2) + assert_allclose(f1, f3) + assert_allclose(p1, p3) + class TestLombscargle: def test_frequency(self): """Test if frequency location of peak corresponds to frequency of