Signal processing

In a wide-array of applications, the original data recorded from a digital system must be manipulated in order to extract the greatest amount of information. GWpy provides a suite of functions to simplify and extend the excellent digital signal processing suite in scipy.signal.

Spectral density estimation

Spectral density estimation is a common way of investigating the frequency-domain content of a time-domain signal. GWpy provides wrappers of power spectral density (PSD) estimation methods from scipy.signal to simplify calculating a FrequencySeries from a TimeSeries.

The gwpy.signal.spectral sub-package provides the following PSD estimation averaging methods:

• 'bartlett' - mean average of non-overlapping periodograms

• 'median' - median average of overlapping periodograms

• 'welch' - mean average of overlapping periodograms

Each of these can be specified by passing the function name as the method keyword argument to any of the relevant TimeSeries instance methods:

TimeSeries.psd([fftlength, overlap, window, ...])	Calculate the PSD FrequencySeries for this TimeSeries
TimeSeries.asd([fftlength, overlap, window, ...])	Calculate the ASD FrequencySeries of this TimeSeries
TimeSeries.spectrogram(stride[, fftlength, ...])	Calculate the average power spectrogram of this TimeSeries using the specified average spectrum method.
TimeSeries.spectrogram2(fftlength[, ...])	Calculate the non-averaged power Spectrogram of this TimeSeries

e.g, TimeSeries.psd():

>>> ts = TimeSeries(...)
>>> psd = ts.psd(..., method='median', ...)

See scipy.signal.welch() for more detailed documentation on the PSD estimation method used.

Time-domain filtering

The TimeSeries object comes with a number of instance methods that should make filtering data trivial for a number of common use cases. Available methods include:

TimeSeries.highpass	Filter this TimeSeries with a high-pass filter.
TimeSeries.lowpass	Filter this TimeSeries with a Butterworth low-pass filter.
TimeSeries.bandpass	Filter this TimeSeries with a band-pass filter.
TimeSeries.zpk	Filter this TimeSeries by applying a zero-pole-gain filter
TimeSeries.whiten	Whiten this TimeSeries using inverse spectrum truncation
TimeSeries.filter	Filter this TimeSeries with an IIR or FIR filter

Each of the above methods eventually calls out to TimeSeries.filter() to apply a digital linear filter, normally via cascaded second-order-sections (requires scipy >= 0.16).

For a worked example of how to filter LIGO data to discover a gravitational-wave signal, see the example Filtering a TimeSeries to detect gravitational waves.

Frequency-domain filtering

Additionally, the TimeSeries object includes a number of instance methods to generate frequency-domain information for some data. Available methods include:

TimeSeries.psd	Calculate the PSD FrequencySeries for this TimeSeries
TimeSeries.asd	Calculate the ASD FrequencySeries of this TimeSeries
TimeSeries.spectrogram	Calculate the average power spectrogram of this TimeSeries using the specified average spectrum method.
TimeSeries.q_transform	Scan a TimeSeries using the multi-Q transform and return an interpolated high-resolution spectrogram
TimeSeries.rayleigh_spectrum	Calculate the Rayleigh FrequencySeries for this TimeSeries.
TimeSeries.rayleigh_spectrogram	Calculate the Rayleigh statistic spectrogram of this TimeSeries

For a worked example of how to load data and calculate the Amplitude Spectral Density FrequencySeries, see the example Calculating and plotting a FrequencySeries.

Filter design

The gwpy.signal provides a number of filter design methods which, when combined with the BodePlot visualisation, can be used to create a number of common filters:

lowpass	Design a low-pass filter for the given cutoff frequency
highpass	Design a high-pass filter for the given cutoff frequency
bandpass	Design a band-pass filter for the given cutoff frequencies
notch	Design a ZPK notch filter for the given frequency and sampling rate
concatenate_zpks	Concatenate a list of zero-pole-gain (ZPK) filters

Each of these will return filter coefficients that can be passed directly into zpk (default for analogue filters) or filter (default for digital filters).

For a worked example of how to filter LIGO data to discover a gravitational-wave signal, see the example Filtering a TimeSeries to detect gravitational waves.

Cross-channel correlations:

TimeSeries.coherence	Calculate the frequency-coherence between this TimeSeries and another.
TimeSeries.coherence_spectrogram	Calculate the coherence spectrogram between this TimeSeries and other.

For a worked example of how to compare channels like this, see the example Calculating the coherence between two channels.

Reference/API

filter_design.bandpass(fhigh, sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)

Design a band-pass filter for the given cutoff frequencies

Parameters

flow : float

lower corner frequency of pass band

fhigh : float

upper corner frequency of pass band

sample_rate : float

sampling rate of target data

fstop : tuple of float, optional

(low, high) edge-frequencies of stop band

gpass : float, optional, default: 2

the maximum loss in the passband (dB)

gstop : float, optional, default: 30

the minimum attenuation in the stopband (dB)

type : str, optional, default: 'iir'

the filter type, either 'iir' or 'fir'

**kwargs

other keyword arguments are passed directly to iirdesign() or firwin()

Returns

filter

the formatted filter. the output format for an IIR filter depends on the input arguments, default is a tuple of (zeros, poles, gain)

Notes

By default a digital filter is returned, meaning the zeros and poles are given in the Z-domain in units of radians/sample.

Examples

To create a band-pass filter for 100-1000 Hz for 4096 Hz-sampled data:

>>> from gwpy.signal.filter_design import bandpass
>>> bp = bandpass(100, 1000, 4096)

To view the filter, you can use the BodePlot:

>>> from gwpy.plot import BodePlot
>>> plot = BodePlot(bp, sample_rate=4096)
>>> plot.show()

(png)

filter_design.lowpass(sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)

Design a low-pass filter for the given cutoff frequency

Parameters

frequency : float

corner frequency of low-pass filter (Hertz)

sample_rate : float

sampling rate of target data (Hertz)

fstop : float, optional

edge-frequency of stop-band (Hertz)

gpass : float, optional, default: 2

the maximum loss in the passband (dB)

gstop : float, optional, default: 30

the minimum attenuation in the stopband (dB)

type : str, optional, default: 'iir'

the filter type, either 'iir' or 'fir'

**kwargs

other keyword arguments are passed directly to iirdesign() or firwin()

Returns

filter

the formatted filter. the output format for an IIR filter depends on the input arguments, default is a tuple of (zeros, poles, gain)

Notes

By default a digital filter is returned, meaning the zeros and poles are given in the Z-domain in units of radians/sample.

Examples

To create a low-pass filter at 1000 Hz for 4096 Hz-sampled data:

>>> from gwpy.signal.filter_design import lowpass
>>> lp = lowpass(1000, 4096)

To view the filter, you can use the BodePlot:

>>> from gwpy.plot import BodePlot
>>> plot = BodePlot(lp, sample_rate=4096)
>>> plot.show()

(png)

filter_design.highpass(sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)

Design a high-pass filter for the given cutoff frequency

Parameters

frequency : float

corner frequency of high-pass filter

sample_rate : float

sampling rate of target data

fstop : float, optional

edge-frequency of stop-band

gpass : float, optional, default: 2

the maximum loss in the passband (dB)

gstop : float, optional, default: 30

the minimum attenuation in the stopband (dB)

type : str, optional, default: 'iir'

the filter type, either 'iir' or 'fir'

**kwargs

other keyword arguments are passed directly to iirdesign() or firwin()

Returns

filter

the formatted filter. the output format for an IIR filter depends on the input arguments, default is a tuple of (zeros, poles, gain)

Notes

By default a digital filter is returned, meaning the zeros and poles are given in the Z-domain in units of radians/sample.

Examples

To create a high-pass filter at 100 Hz for 4096 Hz-sampled data:

>>> from gwpy.signal.filter_design import highpass
>>> hp = highpass(100, 4096)

To view the filter, you can use the BodePlot:

>>> from gwpy.plot import BodePlot
>>> plot = BodePlot(hp, sample_rate=4096)
>>> plot.show()

(png)

filter_design.notch(sample_rate, type='iir', output='zpk', **kwargs)

Design a ZPK notch filter for the given frequency and sampling rate

Parameters

frequency : float, Quantity

frequency (default in Hertz) at which to apply the notch

sample_rate : float, Quantity

number of samples per second for TimeSeries to which this notch filter will be applied

type : str, optional, default: ‘iir’

type of filter to apply, currently only ‘iir’ is supported

output : str, optional, default: ‘zpk’

output format for notch

**kwargs

other keyword arguments to pass to scipy.signal.iirdesign

Returns

filter

the formatted filter; the output format for an IIR filter depends on the input arguments, default is a tuple of (zeros, poles, gain)

See also

scipy.signal.iirdesign

for details on the IIR filter design method and the output formats

Notes

By default a digital filter is returned, meaning the zeros and poles are given in the Z-domain in units of radians/sample.

Examples

To create a low-pass filter at 1000 Hz for 4096 Hz-sampled data:

>>> from gwpy.signal.filter_design import notch
>>> n = notch(100, 4096)

To view the filter, you can use the BodePlot:

>>> from gwpy.plot import BodePlot
>>> plot = BodePlot(n, sample_rate=4096)
>>> plot.show()

(png)

filter_design.concatenate_zpks()

Concatenate a list of zero-pole-gain (ZPK) filters

Parameters

*zpks

one or more zero-pole-gain format, each one should be a 3-tuple containing an array of zeros, an array of poles, and a gain float

Returns

zeros : numpy.ndarray

the concatenated array of zeros

poles : numpy.ndarray

the concatenated array of poles

gain : float

the overall gain

Examples

Create a lowpass and a highpass filter, and combine them:

>>> from gwpy.signal.filter_design import (
...     highpass, lowpass, concatenate_zpks)
>>> hp = highpass(100, 4096)
>>> lp = lowpass(1000, 4096)
>>> zpk = concatenate_zpks(hp, lp)

Plot the filter:

>>> from gwpy.plot import BodePlot
>>> plot = BodePlot(zpk, sample_rate=4096)
>>> plot.show()

(png)