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(self[, fftlength, overlap, …])

Calculate the PSD FrequencySeries for this TimeSeries

TimeSeries.asd(self[, fftlength, overlap, …])

Calculate the ASD FrequencySeries of this TimeSeries

TimeSeries.spectrogram(self, stride[, …])

Calculate the average power spectrogram of this TimeSeries using the specified average spectrum method.

TimeSeries.spectrogram2(self, 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(flow, fhigh, sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)[source]

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)

../_images/index-12.png
filter_design.lowpass(frequency, sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)[source]

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)

../_images/index-22.png
filter_design.highpass(frequency, sample_rate, fstop=None, gpass=2, gstop=30, type='iir', **kwargs)[source]

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)

../_images/index-32.png
filter_design.notch(frequency, sample_rate, type='iir', **kwargs)[source]

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

**kwargs

other keyword arguments to pass to scipy.signal.iirdesign

Returns

zpk : tuple of complex or float

the filter components in digital zero-pole-gain format

See also

scipy.signal.iirdesign

for details on the IIR filter design method

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)

../_images/index-42.png
filter_design.concatenate_zpks(*zpks)[source]

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)

../_images/index-52.png