From 0fa4999d7b83e9e73535fa0548aab9277cdc40c8 Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 12:20:02 +0000 Subject: [PATCH 1/6] signal.fft.lal: use CreateNamed{}Window from lal easier UI --- gwpy/signal/fft/lal.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/gwpy/signal/fft/lal.py b/gwpy/signal/fft/lal.py index ab5f93d21..c66ad3aef 100644 --- a/gwpy/signal/fft/lal.py +++ b/gwpy/signal/fft/lal.py @@ -31,6 +31,7 @@ import numpy from ...frequencyseries import FrequencySeries +from ..window import canonical_name from .utils import scale_timeseries_unit from . import registry as fft_registry @@ -121,14 +122,13 @@ def generate_window(length, window=('kaiser', 24), dtype='float64'): except KeyError: # parse window as name and arguments, e.g. ('kaiser', 24) if isinstance(window, (list, tuple)): - args = window[1:] - window = str(window[0]) + window, beta = window else: - args = [] - window = window.title() if window.islower() else window + beta = 0 + window = canonical_name(window) # create window - create = getattr(lal, 'Create%s%sWindow' % (window, laltype)) - LAL_WINDOWS[key] = create(length, *args) + create = getattr(lal, 'CreateNamed{}Window'.format(laltype)) + LAL_WINDOWS[key] = create(window, beta, length) return LAL_WINDOWS[key] From 32b49bb717fcf4cef55cc4cdc4ee2be2bc5d0d52 Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 12:22:58 +0000 Subject: [PATCH 2/6] signal.fft.lal: handle window as array --- gwpy/signal/fft/lal.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/gwpy/signal/fft/lal.py b/gwpy/signal/fft/lal.py index c66ad3aef..773763256 100644 --- a/gwpy/signal/fft/lal.py +++ b/gwpy/signal/fft/lal.py @@ -132,6 +132,21 @@ def generate_window(length, window=('kaiser', 24), dtype='float64'): return LAL_WINDOWS[key] +def window_from_array(array): + """Convert a `numpy.ndarray` into a LAL `Window` object + """ + import lal + from ...utils.lal import LAL_TYPE_STR_FROM_NUMPY + laltype = LAL_TYPE_STR_FROM_NUMPY[array.dtype.type] + + # create sequence + seq = getattr(lal, 'Create{}Sequence'.format(laltype))(array.size) + seq.data = array + + # create window from sequence + return getattr(lal, 'Create{}WindowFromSequence'.format(laltype))(seq) + + # -- spectrumm methods ------------------------------------------------------ def _lal_spectrum(timeseries, segmentlength, noverlap=None, method='welch', @@ -177,6 +192,9 @@ def _lal_spectrum(timeseries, segmentlength, noverlap=None, method='welch', elif isinstance(window, (tuple, str)): window = generate_window(segmentlength, window=window, dtype=timeseries.dtype) + elif isinstance(window, numpy.ndarray): + window = window_from_array(window) + # get FFT plan if plan is None: plan = generate_fft_plan(segmentlength, dtype=timeseries.dtype) From 7d1620828f8a7702654bca048dd7d9920f8b65e6 Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 13:10:50 +0000 Subject: [PATCH 3/6] signal.fft.lal: modified default window allow `window=None` and set to kaiser --- gwpy/signal/fft/lal.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gwpy/signal/fft/lal.py b/gwpy/signal/fft/lal.py index 773763256..b0ae75553 100644 --- a/gwpy/signal/fft/lal.py +++ b/gwpy/signal/fft/lal.py @@ -87,7 +87,7 @@ def generate_fft_plan(length, level=None, dtype='float64', forward=True): return LAL_FFTPLANS[key] -def generate_window(length, window=('kaiser', 24), dtype='float64'): +def generate_window(length, window=None, dtype='float64'): """Generate a time-domain window for use in a LAL FFT Parameters @@ -111,6 +111,9 @@ def generate_window(length, window=('kaiser', 24), dtype='float64'): import lal from ...utils.lal import LAL_TYPE_STR_FROM_NUMPY + if window is None: + window = ('kaiser', 24) + # generate key for caching window laltype = LAL_TYPE_STR_FROM_NUMPY[numpy.dtype(dtype).type] key = (length, str(window), laltype) From d0410bb63ef1fb827e9f455df0f91560ade46e17 Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 13:11:38 +0000 Subject: [PATCH 4/6] signal.fft.ui: parse window during normalize func allows single place to generate windows for all libraries --- gwpy/signal/fft/lal.py | 11 ++----- gwpy/signal/fft/ui.py | 66 ++++++++++++++++++++++-------------------- 2 files changed, 38 insertions(+), 39 deletions(-) diff --git a/gwpy/signal/fft/lal.py b/gwpy/signal/fft/lal.py index b0ae75553..cca0f9a4d 100644 --- a/gwpy/signal/fft/lal.py +++ b/gwpy/signal/fft/lal.py @@ -170,10 +170,10 @@ def _lal_spectrum(timeseries, segmentlength, noverlap=None, method='welch', noverlap : `int` number of samples to overlap between segments, defaults to 50%. - window : `tuple`, `str`, optional - window parameters to apply to timeseries prior to FFT + window : `lal.REAL8Window`, optional + window to apply to timeseries prior to FFT - plan : `REAL8FFTPlan`, optional + plan : `lal.REAL8FFTPlan`, optional LAL FFT plan to use when generating average spectrum Returns @@ -192,11 +192,6 @@ def _lal_spectrum(timeseries, segmentlength, noverlap=None, method='welch', # get window if window is None: window = generate_window(segmentlength, dtype=timeseries.dtype) - elif isinstance(window, (tuple, str)): - window = generate_window(segmentlength, window=window, - dtype=timeseries.dtype) - elif isinstance(window, numpy.ndarray): - window = window_from_array(window) # get FFT plan if plan is None: diff --git a/gwpy/signal/fft/ui.py b/gwpy/signal/fft/ui.py index 26c670cdf..924621a41 100644 --- a/gwpy/signal/fft/ui.py +++ b/gwpy/signal/fft/ui.py @@ -70,7 +70,7 @@ def seconds_to_samples(x, rate): return int((Quantity(x, 's') * rate).decompose().value) -def normalize_fft_params(series, kwargs=None): +def normalize_fft_params(series, kwargs=None, library=None): """Normalize a set of FFT parameters for processing This method reads the ``fftlength`` and ``overlap`` keyword arguments @@ -89,6 +89,10 @@ def normalize_fft_params(series, kwargs=None): kwargs : `dict` the dict of keyword arguments passed by the user + library: `str`, optional + the name of the library that provides the FFT methods, e.g. + 'scipy' + Examples -------- >>> from numpy.random import normal @@ -104,10 +108,11 @@ def normalize_fft_params(series, kwargs=None): samp = series.sample_rate fftlength = kwargs.pop('fftlength', None) overlap = kwargs.pop('overlap', None) + window = kwargs.pop('window', None) # get canonical window name - if isinstance(kwargs.get('window', None), str): - kwargs['window'] = canonical_name(kwargs['window']) + if isinstance(window, str): + window = canonical_name(window) # fftlength -> nfft if fftlength is None: @@ -115,13 +120,32 @@ def normalize_fft_params(series, kwargs=None): nfft = seconds_to_samples(fftlength, samp) # overlap -> noverlap - if overlap is None and isinstance(kwargs.get('window', None), str): - noverlap = recommended_overlap(kwargs['window'], nfft) + if overlap is None and isinstance(window, str): + noverlap = recommended_overlap(window, nfft) elif overlap is None: noverlap = 0 else: noverlap = seconds_to_samples(overlap, samp) + # create window + if library == 'lal' and isinstance(window, numpy.ndarray): + from .lal import window_from_array + window = window_from_array(window) + elif library == 'lal': + from .lal import generate_window + window = generate_window(nfft, window=window, dtype=series.dtype) + elif isinstance(window, (str, tuple)): + window = get_window(window, nfft) + + # allow FFT methods to use their own defaults + if window is not None: + kwargs['window'] = window + + # create FFT plan for LAL + if library == 'lal' and kwargs.get('plan', None): + from .lal import generate_fft_plan + kwargs['plan'] = generate_fft_plan(nfft, dtype=series.dtype) + kwargs.update({ 'nfft': nfft, 'noverlap': noverlap, @@ -141,8 +165,9 @@ def wrapped_func(series, method_func, *args, **kwargs): else: data = series - # extract parameters in seconds, setting recommended default overlap - normalize_fft_params(data, kwargs) + # normalise FFT parmeters for all libraries + library = method_func.__module__.rsplit('.', 1)[-1] + normalize_fft_params(data, kwargs=kwargs, library=library) return func(series, method_func, *args, **kwargs) @@ -217,10 +242,10 @@ def average_spectrogram(timeseries, method_func, stride, *args, **kwargs): epoch = timeseries.t0.value nstride = seconds_to_samples(stride, timeseries.sample_rate) kwargs['fftlength'] = kwargs.pop('fftlength', stride) or stride - normalize_fft_params(timeseries, kwargs) + normalize_fft_params(timeseries, kwargs=kwargs, + library=method_func.__module__.rsplit('.', 1)[-1]) nfft = kwargs['nfft'] noverlap = kwargs['noverlap'] - window = kwargs.pop('window', None) # sanity check parameters if nstride > timeseries.size: @@ -231,21 +256,6 @@ def average_spectrogram(timeseries, method_func, stride, *args, **kwargs): if noverlap >= nfft: raise ValueError("overlap must be less than fftlength") - # generate windows and FFT plans up-front - if method_func.__module__.endswith('.lal'): - from .lal import (generate_fft_plan, generate_window) - if isinstance(window, (str, tuple)): - window = generate_window(nfft, window=window, - dtype=timeseries.dtype) - if kwargs.get('plan', None) is None: - kwargs['plan'] = generate_fft_plan(nfft, dtype=timeseries.dtype) - else: - if isinstance(window, (str, tuple)): - window = get_window(window, nfft) - # don't operate on None, let the method_func work out its own defaults - if window is not None: - kwargs['window'] = window - # set up single process Spectrogram method def _psd(series): """Calculate a single PSD for a spectrogram @@ -293,18 +303,12 @@ def spectrogram(timeseries, method_func, **kwargs): if noverlap >= nfft: raise ValueError("overlap must be less than fftlength") - # get window once (if given) - window = kwargs.pop('window', None) or 'hann' - if isinstance(window, (str, tuple)): - window = get_window(window, nfft) - # set up single process Spectrogram method def _psd(series): """Calculate a single PSD for a spectrogram """ try: - return method_func(series, nfft=nfft, window=window, - **kwargs)[1] + return method_func(series, nfft=nfft, **kwargs)[1] except Exception as exc: # pylint: disable=broad-except if nproc == 1: raise From 27813c83f08d42c515b9405c3b4c285550fd62d5 Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 15:29:04 +0000 Subject: [PATCH 5/6] TimeSeries.spectrogram: default to overlap=None so that it gets set to the recommended overlap for the window --- gwpy/timeseries/timeseries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwpy/timeseries/timeseries.py b/gwpy/timeseries/timeseries.py index 9d57e2d85..63b142782 100644 --- a/gwpy/timeseries/timeseries.py +++ b/gwpy/timeseries/timeseries.py @@ -368,7 +368,7 @@ def csd(self, other, fftlength=None, overlap=None, window='hann', overlap=overlap, window=window, **kwargs) @_update_doc_with_fft_methods - def spectrogram(self, stride, fftlength=None, overlap=0, + def spectrogram(self, stride, fftlength=None, overlap=None, window='hann', method='scipy-welch', nproc=1, **kwargs): """Calculate the average power spectrogram of this `TimeSeries` using the specified average spectrum method. From 804f208d4815318830bbf9f1a7b734712b7840bf Mon Sep 17 00:00:00 2001 From: Duncan Macleod Date: Wed, 21 Feb 2018 15:29:38 +0000 Subject: [PATCH 6/6] tests: improved tests of TimeSeries spectral methods --- gwpy/tests/test_timeseries.py | 205 +++++++++++++++++++++------------- 1 file changed, 129 insertions(+), 76 deletions(-) diff --git a/gwpy/tests/test_timeseries.py b/gwpy/tests/test_timeseries.py index 51b7fe15c..e5e83a5cf 100644 --- a/gwpy/tests/test_timeseries.py +++ b/gwpy/tests/test_timeseries.py @@ -19,9 +19,11 @@ """Unit test for timeseries module """ +import importlib import os import pytest import tempfile +from itertools import (chain, product) from six.moves.urllib.request import urlopen from six.moves.urllib.error import URLError @@ -844,57 +846,77 @@ def test_average_fft(self, losc): fs = losc.average_fft(fftlength=0.4, overlap=0.2) - def test_psd_simple(self, losc): - # test all defaults - fs = losc.psd() - assert isinstance(fs, FrequencySeries) - assert fs.size == losc.size // 2 + 1 - assert fs.f0 == 0 * units.Hz - assert fs.df == 1 / losc.duration - assert fs.channel is losc.channel - assert fs.unit == losc.unit ** 2 / units.Hz - - # test fftlength - fs = losc.psd(fftlength=0.5) - assert fs.size == 0.5 * losc.sample_rate.value // 2 + 1 - assert fs.df == 2 * units.Hz - - # test overlap - fs = losc.psd(fftlength=0.4, overlap=0.2) - - # test default overlap - fs2 = losc.psd(fftlength=.4) - utils.assert_quantity_sub_equal(fs, fs2) - - @pytest.mark.parametrize('library', ( - pytest.param('pycbc', - marks=utils.skip_missing_dependency('pycbc.psd')), - pytest.param('lal', marks=utils.skip_missing_dependency('lal')), + @pytest.mark.parametrize('library, method', chain( + product([None], ['welch', 'bartlett']), + product(['scipy'], ['welch', 'bartlett']), + product(['pycbc'], ['welch', 'bartlett', 'median', 'median_mean']), + product(['lal'], ['welch', 'bartlett', 'median', 'median_mean']), )) @pytest.mark.parametrize( - 'method', ('welch', 'bartlett', 'median', 'median_mean'), + 'window', (None, 'hann', ('kaiser', 24), 'array'), ) - def test_psd_library(self, losc, library, method): - method = '{}_{}'.format(library, method) + def test_psd(self, losc, library, method, window): + if library: + try: + importlib.import_module(library) + except ImportError as exc: + pytest.skip(str(exc)) + + method = '{}_{}'.format(library, method) - def _psd(): - if method == 'lal_median_mean': + def _psd(fftlength, overlap, **kwargs): + # check number of segments for LAL median-mean + if fftlength is None: + nfft = losc.size + else: + nfft = int(fftlength * losc.sample_rate.value) + if overlap is None: + noverlap = int(nfft/2.) + else: + noverlap = int(overlap * losc.sample_rate.value) + nstep = nfft - noverlap + req = (int((losc.size - nfft) / nstep) * nstep + nfft) + if method == 'lal_median_mean' and req != losc.size: # LAL should warn about the data being the wrong length warnctx = pytest.warns(UserWarning) + elif library is None and method is not None: + # not specifying library should print warning + warnctx = pytest.warns(UserWarning) else: warnctx = null_context() + + # create window of the correct length + if window == 'array': + _window = signal.get_window('hamming', nfft) + else: + _window = window + + # generate PSD with warnctx: - return losc.psd(fftlength=.5, overlap=.25, method=method) + return losc.psd(fftlength=fftlength, overlap=overlap, + method=method, window=_window) - # check simple - psd = _psd() - assert isinstance(psd, FrequencySeries) - assert psd.f0 == 0 * units.Hz + # test basic + if method.endswith('median_mean'): + ctx = pytest.raises(ValueError) + else: + ctx = null_context() + with ctx: + fs = _psd(None, None) + assert fs.size == losc.size // 2 + 1 + assert fs.f0 == 0 * units.Hz + assert fs.df == 1 / losc.duration + assert fs.channel is losc.channel + assert fs.unit == losc.unit ** 2 / units.Hz + + # test fftlength maps to yindex + psd = _psd(.5, None) + assert psd.size == 0.5 * losc.sample_rate.value // 2 + 1 assert psd.df == 2 * units.Hz - # check window selection - if library != 'pycbc': - _psd() + # test default overlap + if library is None and method == 'welch' and window == 'hann': + utils.assert_quantity_sub_equal(psd, _psd(.5, .25)) def test_asd(self, losc): fs = losc.asd() @@ -914,11 +936,42 @@ def test_csd(self, losc): # test overlap losc.csd(losc, fftlength=0.4, overlap=0.2) - def test_spectrogram(self, losc): + @pytest.mark.parametrize('library, method', chain( + product([None], ['welch', 'bartlett']), + product(['scipy'], ['welch', 'bartlett']), + product(['pycbc'], ['welch', 'bartlett', 'median', 'median_mean']), + product(['lal'], ['welch', 'bartlett', 'median', 'median_mean']), + )) + @pytest.mark.parametrize( + 'window', (None, 'hann', ('kaiser', 24), 'array'), + ) + def test_spectrogram(self, losc, library, method, window): + if library: + try: + importlib.import_module(library) + except ImportError as exc: + pytest.skip(str(exc)) + method = '{}_{}'.format(library, method) + ctx = null_context + else: + ctx = lambda: pytest.warns(UserWarning) + + def _spectrogram(*args, **kwargs): + kwargs.setdefault('method', method) + if window == 'array': + nfft = int(losc.sample_rate.value * ( + kwargs.get('fftlength', args[0]) or args[0])) + w = signal.get_window('hamming', nfft) + else: + w = window + kwargs.setdefault('window', w) + with ctx(): + return losc.spectrogram(*args, **kwargs) + # test defaults - sg = losc.spectrogram(1) # defaults to 50% overlap for 'hann' window + sg = _spectrogram(1) assert isinstance(sg, Spectrogram) - assert sg.shape == (4, losc.sample_rate.value // 2 + 1) + assert sg.shape == (abs(losc.span), losc.sample_rate.value // 2 + 1) assert sg.f0 == 0 * units.Hz assert sg.df == 1 * units.Hz assert sg.channel is losc.channel @@ -927,56 +980,56 @@ def test_spectrogram(self, losc): assert sg.span == losc.span # check the same result as PSD - psd = losc[:int(losc.sample_rate.value)].psd() + with ctx(): + if window == 'array': + win = signal.get_window('hamming', int(losc.sample_rate.value)) + else: + win = window + n = int(losc.sample_rate.value) + overlap = 0 + if window in {'hann'}: + overlap = .5 + n += int(overlap * losc.sample_rate.value) + psd = losc[:n].psd(fftlength=1, overlap=overlap, + method=method, window=win) # FIXME: epoch should not be excluded here (probably) utils.assert_quantity_sub_equal(sg[0], psd, exclude=['epoch']) # test fftlength - sg = losc.spectrogram(1, fftlength=0.5) - assert sg.shape == (4, 0.5 * losc.sample_rate.value // 2 + 1) + sg = _spectrogram(1, fftlength=0.5) + assert sg.shape == (abs(losc.span), + 0.5 * losc.sample_rate.value // 2 + 1) assert sg.df == 2 * units.Hertz assert sg.dt == 1 * units.second + # test overlap - sg = losc.spectrogram(0.5, fftlength=0.25, overlap=0.125) - assert sg.shape == (8, 0.25 * losc.sample_rate.value // 2 + 1) - assert sg.df == 4 * units.Hertz - assert sg.dt == 0.5 * units.second + if window == 'hann': + sg2 = _spectrogram(1, fftlength=0.5, overlap=.25) + utils.assert_quantity_sub_equal(sg, sg2) + # test multiprocessing - sg2 = losc.spectrogram(0.5, fftlength=0.25, overlap=0.125, nproc=2) + sg2 = _spectrogram(1, fftlength=0.5, nproc=2) utils.assert_quantity_sub_equal(sg, sg2) - # test a couple of methods - try: - with pytest.warns(UserWarning): - sg = losc.spectrogram(0.5, fftlength=0.25, method='welch') - except TypeError: # old pycbc doesn't accept window as array - pass - else: - assert sg.shape == (8, 0.25 * losc.sample_rate.value // 2 + 1) - assert sg.df == 4 * units.Hertz - assert sg.dt == 0.5 * units.second - - sg = losc.spectrogram(0.5, fftlength=0.25, method='scipy-bartlett') - assert sg.shape == (8, 0.25 * losc.sample_rate.value // 2 + 1) - assert sg.df == 4 * units.Hertz - assert sg.dt == 0.5 * units.second - # check that `cross` keyword gets deprecated properly # TODO: removed before 1.0 release - with pytest.warns(DeprecationWarning) as wng: - try: - out = losc.spectrogram(0.5, fftlength=.25, cross=losc) - except AttributeError: - return # scipy is too old - assert '`cross` keyword argument has been deprecated' in \ - wng[0].message.args[0] - utils.assert_quantity_sub_equal( - out, losc.csd_spectrogram(losc, 0.5, fftlength=.25)) + if method == 'scipy_welch' and window is None: + with pytest.warns(DeprecationWarning) as wng: + try: + out = _spectrogram(0.5, fftlength=.25, cross=losc) + except AttributeError: + return # scipy is too old + assert '`cross` keyword argument has been deprecated' in \ + wng[0].message.args[0] + utils.assert_quantity_sub_equal( + out, losc.csd_spectrogram(losc, 0.5, fftlength=.25)) def test_spectrogram2(self, losc): # test defaults sg = losc.spectrogram2(1) - utils.assert_quantity_sub_equal(sg, losc.spectrogram(1)) + utils.assert_quantity_sub_equal( + sg, losc.spectrogram(1, fftlength=1, overlap=0, + method='scipy-welch', window='boxcar')) # test fftlength sg = losc.spectrogram2(0.5)