Module ButterworthTimeSeries.c

Detailed Description

Applies a low- or high-pass Butterworth filter to a time series.

Author

Creighton, T. D.

Synopsis

The Module ButterworthTimeSeries.c provides specific advice and guidelines for building a numerically stable time-domain filter, but the general procedure is as follows.

1. Decide on the desired filter response, and express it as a rational function of the frequency variable \(w=\tan(\pi f\Delta t)\) (which maps the Nyquist interval onto the positive real axis).
2. Factorize this rational function into zeros and poles, restricting oneself to the upper half of the \(w\) complex plane. Assign these to one or more objects of type <datatype>ZPGFilter
3. Use WToZ<datatype>ZPGFilter() to transform the filter to the more conventional \(z=\exp(2\pi if\Delta t)\) frequency variable.
4. Use the routines in Header IIRFilter.h to create IIR filters from the ZPG filters, and to apply them to data.

Description

These routines perform an in-place time-domain band-pass filtering of a data sequence *series, using a Butterworth filter generated from parameters *params. The routines construct a filter with the square root of the desired amplitude response, which it then applied to the data once forward and once in reverse. This gives the full amplitude response with little or no frequency-dependent phase shift.

The routine LALDButterworthREAL4TimeSeries() applies a double-precision filter to single-precision data, using LALDIIRFilterREAL4Vector() and LALDIIRFilterREAL4VectorR().

Algorithm

The frequency response of a Butterworth low-pass filter is easiest to express in terms of the transformed frequency variable \(w=\tan(\pi f\Delta t)\), where \(\Delta t\) is the sampling interval (i.e. series->deltaT). In this parameter, then, the power response (attenuation) of the filter is:

\[ |R|^2 = \sqrt{a} = \frac{1}{1+(w/w_c)^{2n}} \; , \]

where \(n\) is the filter order and \(w_c\) is the characteristic frequency. We have written the attenuation as \(\sqrt{a}\) to emphasize that the full attenuation \(a\) is achieved only after filtering twice (once forward, once in reverse). Similarly, a Butterworth high-pass filter is given by

\[ |R|^2 = \sqrt{a} = \frac{1}{1+(w_c/w)^{2n}} \; . \]

If one is given a filter order \(n\), then the characteristic frequency can be determined from the attenuation at some any given frequency. Alternatively, \(n\) and \(w_c\) can both be computed given attenuations at two different frequencies.

Frequencies in *params are assumed to be real frequencies \(f\) given in the inverse of the units used for the sampling interval series->deltaT. In order to be used, the pass band parameters must lie in the ranges given below; if a parameter lies outside of its range, then it is ignored and the filter is calculated from the remaining parameters. If too many parameters are missing, the routine will fail. The acceptable parameter ranges are:

params->nMax

= 1, 2, \(\ldots\)

params->f1, f2

\(\in(0,\{2\times<tt>series-\>deltaT</tt>\}^{-1}) \)

params->a1, a2

\(\in (0,1) \)

If both pairs of frequencies and amplitudes are given, then a1, a2 specify the minimal requirements on the attenuation of the filter at frequencies f1, f2. Whether the filter is a low- or high-pass filter is determined from the relative sizes of these parameters. In this case the nMax parameter is optional; if given, it specifies an upper limit on the filter order. If the desired attenuations would require a higher order, then the routine will sacrifice performance in the stop band in order to remain within the specified nMax.

If one of the frequency/attenuation pairs is missing, then the filter is computed using the remaining pair and nMax (which must be given). The filter is taken to be a low-pass filter if f1, a1 are given, and high-pass if f2, a2 are given. If only one frequency and no corresponding attenuation is specified, then it is taken to be the characteristic frequency (i.e. the corresponding attenuation is assumed to be \(\sqrt{a}=1/2\)). If none of these conditions are met, the routine will return an error.

The following table summarizes the decision algorithm. A \(\bullet\) symbol indicates that the parameter is specified in the range given above. A \(\circ\) symbol indicates that the parameter is not given, i.e. not specified in the valid range.

nMax	f1	a1	f2	a2	Procedure
\(\circ\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	Type of filter (low- or high-pass), \(w_c\), and \(n\) are computed from all four transition-band parameters.
\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	Ditto, but if the resulting \(n>\) nMax, \(w_c\) is computed from nMax and the (f,a) pair with the larger a.
\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\circ\)	\(\circ\)	Low-pass filter; \(w_c\) is computed from nMax, f1, and a1.
\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\circ\)	\(\bullet\)	Ditto; a2 is ignored.
\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	\(\circ\)	Ditto; f2 is ignored.
\(\bullet\)	\(\bullet\)	\(\circ\)	\(\circ\)	\(\circ\)	Low-pass filter; \(w_c\) is computed as above with a1 treated as 1/4.
\(\bullet\)	\(\bullet\)	\(\circ\)	\(\circ\)	\(\bullet\)	Ditto; a2 is ignored.
\(\bullet\)	\(\circ\)	\(\circ\)	\(\bullet\)	\(\bullet\)	High-pass filter; \(w_c\) is computed from nMax, f2, and a2.
\(\bullet\)	\(\circ\)	\(\bullet\)	\(\bullet\)	\(\bullet\)	Ditto; a1 is ignored.
\(\bullet\)	\(\bullet\)	\(\circ\)	\(\bullet\)	\(\bullet\)	Ditto; f1 is ignored.
\(\bullet\)	\(\circ\)	\(\circ\)	\(\bullet\)	\(\circ\)	High-pass filter; \(w_c\) is computed as above with a2 treated as 1/4.
\(\bullet\)	\(\circ\)	\(\bullet\)	\(\bullet\)	\(\circ\)	Ditto; a1 is ignored.
Other					Subroutine returns an error.

Once an order \(n\) and characteristic frequency \(w_c\) are known, the zeros and poles of a ZPG filter are readily determined. A stable, physically realizable Butterworth filter will have \(n\) poles evenly spaced on the upper half of a circle of radius \(w_c\); that is,

\[ R = \frac{(-iw_c)^n}{\prod_{k=0}^{n-1}(w - w_c e^{2\pi i(k+1/2)/n})} \]

for a low-pass filter, and

\[ R = \frac{w^n}{\prod_{k=0}^{n-1}(w - w_c e^{2\pi i(k+1/2)/n})} \]

for a high-pass filter. By choosing only poles on the upper-half plane, one ensures that after transforming to \(z\) the poles will have \(|z|<1\). Furthermore, the phase factor \((-i)^n\) in the numerator of the low-pass filter is chosen so that the DC response is purely real; this ensures that the response function in the \(z\)-plane will have a real gain factor, and the resulting IIR filter will be physically realizable. The high-pass filter has a purely real response at Nyquist ( \(w\rightarrow\infty\)), which similarly gives a physical IIR filter.

Although higher orders \(n\) would appear to produce better (i.e. sharper) filter responses, one rapidly runs into numerical errors, as one ends up adding and subtracting \(n\) large numbers to obtain small filter responses. One way around this is to break the filter up into several lower-order filters. The routines in this module do just that. Poles are paired up across the imaginary axis, (and combined with pairs of zeros at \(w=0\) for high-pass filters,) to form \([n/2]\) second-order filters. If \(n\) is odd, there will be an additional first-order filter, with one pole at \(w=iw_c\) (and one zero at \(w=0\) for a high-pass filter).

Each ZPG filter in the \(w\)-plane is first transformed to the \(z\)-plane by a bilinear transformation, and is then used to construct a time-domain IIR filter. Each filter is then applied to the time series. As mentioned in the description above, the filters are designed to give an overall amplitude response that is the square root of the desired attenuation; however, each time-domain filter is applied to the data stream twice: once in the normal sense, and once in the time-reversed sense. This gives the full attenuation with very little frequency-dependent phase shift.

Prototypes
static int	XLALParsePassBandParamStruc (PassBandParamStruc *params, INT4 *n, REAL8 *wc, REAL8 deltaT)
void	LALButterworthREAL4TimeSeries (LALStatus *stat, REAL4TimeSeries *series, PassBandParamStruc *params)
	Deprecated. More...
void	LALButterworthREAL8TimeSeries (LALStatus *stat, REAL8TimeSeries *series, PassBandParamStruc *params)
	Deprecated. More...
void	LALDButterworthREAL4TimeSeries (LALStatus *stat, REAL4TimeSeries *series, PassBandParamStruc *params)
	Deprecated. More...

Macros
#define	COMPLEX_DATA
#define	SINGLE_PRECISION
#define	SINGLE_PRECISION

Function Documentation

XLALParsePassBandParamStruc()

static int XLALParsePassBandParamStruc ( PassBandParamStruc *  params, INT4 *  n, REAL8 *  wc, REAL8  deltaT  )

static

Definition at line 447 of file ButterworthTimeSeries.c.

LALButterworthREAL4TimeSeries()

void LALButterworthREAL4TimeSeries	(	LALStatus *	stat,
		REAL4TimeSeries *	series,
		PassBandParamStruc *	params
	)

Deprecated.

Deprecated:

Use XLALButterworthREAL4TimeSeries() instead.

Definition at line 240 of file ButterworthTimeSeries.c.

LALButterworthREAL8TimeSeries()

void LALButterworthREAL8TimeSeries	(	LALStatus *	stat,
		REAL8TimeSeries *	series,
		PassBandParamStruc *	params
	)

Deprecated.

Deprecated:

Use XLALButterworthREAL8TimeSeries() instead.

Definition at line 394 of file ButterworthTimeSeries.c.

LALDButterworthREAL4TimeSeries()

void LALDButterworthREAL4TimeSeries	(	LALStatus *	stat,
		REAL4TimeSeries *	series,
		PassBandParamStruc *	params
	)

Deprecated.

Deprecated:

Definition at line 422 of file ButterworthTimeSeries.c.

Macro Definition Documentation

COMPLEX_DATA

#define COMPLEX_DATA

Definition at line 224 of file ButterworthTimeSeries.c.

SINGLE_PRECISION [1/2]

#define SINGLE_PRECISION

Definition at line 230 of file ButterworthTimeSeries.c.

SINGLE_PRECISION [2/2]

#define SINGLE_PRECISION

Definition at line 230 of file ButterworthTimeSeries.c.