.. sectionauthor:: Duncan Macleod <duncan.macleod@ligo.org>

.. currentmodule:: gwpy.timeseries

Filtering a `TimeSeries` with a ZPK filter
##########################################

Several data streams read from the LIGO detectors are whitened before being
recorded to prevent numerical errors when using single-precision data
storage.
In this example we read such `channel <gwpy.detector.Channel>` and undo the
whitening to show the physical content of these data.




First, we import the `TimeSeries` and :meth:`~TimeSeries.get` the data:

.. code-block:: python

   from gwpy.timeseries import TimeSeries
   white = TimeSeries.get(
       'L1:OAF-CAL_DARM_DQ', 'March 2 2015 12:00', 'March 2 2015 12:30')


Now, we can re-calibrate these data into displacement units by first applying
a `highpass <TimeSeries.highpass>` filter to remove the low-frequency noise,
and then applying our de-whitening filter in `ZPK <TimeSeries.zpk>` format
with five zeros at 100 Hz and five poles at 1 Hz (giving an overall DC
gain of 10 :sup:`-10`:

.. code-block:: python

   hp = white.highpass(4)
   displacement = hp.zpk([100]*5, [1]*5, 1e-10)


We can visualise the impact of the whitening by calculating the ASD
`~gwpy.frequencyseries.FrequencySeries` before and after the filter,


.. code-block:: python

   whiteasd = white.asd(8, 4)
   dispasd = displacement.asd(8, 4)


and plotting:


.. code-block:: python

   from gwpy.plotter import FrequencySeriesPlot
   plot = FrequencySeriesPlot(whiteasd, dispasd, sep=True, sharex=True, label=None)


Here we have passed the two
`spectra <gwpy.frequencyseries.FrequencySeries>` in order,
then `sep=True` to display them on separate Axes, `sharex=True` to tie
the `~matplotlib.axis.XAxis` of each of the `~gwpy.plotter.FrequencySeriesAxes`
together, and `label=None` to remove any unwanted legends.

Finally, we prettify our plot with some limits, and some labels:

.. code-block:: python

   plot.axes[0].set_ylabel('ASD [whitened]')
   plot.axes[1].set_ylabel(r'ASD [m/\rtHz]')
   plot.axes[1].set_xlabel('Frequency [Hz]')
   plot.axes[1].set_ylim(1e-20, 1e-15)
   plot.axes[1].set_xlim(5, 4000)
   plot.show()

.. plot:: ../examples/timeseries/filter.py