`Filters` - filter design and filtering

DSP.jl differentiates between filter coefficients and stateful filters. Filter coefficient objects specify the response of the filter in one of several standard forms. Stateful filter objects carry the state of the filter together with filter coefficients in an implementable form (PolynomialRatio, Biquad, or SecondOrderSections). When invoked on a filter coefficient object, filt does not preserve state.

Filter coefficient objects

DSP.jl supports common filter representations. Filter coefficients can be converted from one type to another using convert.

DSP.Filters.ZeroPoleGain — Type.

ZeroPoleGain(z, p, k)

Filter representation in terms of zeros z, poles p, and gain k:

\[H(x) = k\frac{(x - \verb!z[1]!) \ldots (x - \verb!z[end]!)}{(x - \verb!p[1]!) \ldots (x - \verb!p[end]!)}\]

DSP.Filters.PolynomialRatio — Type.

PolynomialRatio(b, a)

Filter representation in terms of the coefficients of the numerator b and denominator a of the transfer function:

\[H(s) = \frac{\verb!b[1]! s^{n-1} + \ldots + \verb!b[n]!}{\verb!a[1]! s^{n-1} + \ldots + \verb!a[n]!}\]

or equivalently:

\[H(z) = \frac{\verb!b[1]! + \ldots + \verb!b[n]! z^{-n+1}}{\verb!a[1]! + \ldots + \verb!a[n]! z^{-n+1}}\]

b and a may be specified as Polynomial objects or vectors ordered from highest power to lowest.

DSP.Filters.Biquad — Type.

Biquad(b0, b1, b2, a1, a2)

Filter representation in terms of the transfer function of a single second-order section given by:

\[H(s) = \frac{\verb!b0! s^2+\verb!b1! s+\verb!b2!}{s^2+\verb!a1! s + \verb!a2!}\]

or equivalently:

\[H(z) = \frac{\verb!b0!+\verb!b1! z^{-1}+\verb!b2! z^{-2}}{1+\verb!a1! z^{-1} + \verb!a2! z^{-2}}\]

DSP.Filters.SecondOrderSections — Type.

SecondOrderSections(biquads, gain)

Filter representation in terms of a cascade of second-order sections and gain. biquads must be specified as a vector of Biquads.

Stateful filter objects

DSP.Filters.DF2TFilter — Type.

DF2TFilter(coef[, si])

Construct a stateful direct form II transposed filter with coefficients coef. si is an optional array representing the initial filter state (defaults to zeros). If f is a PolynomialRatio, Biquad, or SecondOrderSections, filtering is implemented directly. If f is a ZeroPoleGain object, it is first converted to a SecondOrderSections object.

DSP.jl's FIRFilter type maintains state between calls to filt, allowing you to filter a signal of indefinite length in RAM-friendly chunks. FIRFilter contains nothing more that the state of the filter, and a FIRKernel. There are five different kinds of FIRKernel for single rate, up-sampling, down-sampling, rational resampling, and arbitrary sample-rate conversion. You need not specify the type of kernel. The FIRFilter constructor selects the correct kernel based on input parameters.

DSP.Filters.FIRFilter — Type.

FIRFilter(h[, ratio])

Construct a stateful FIRFilter object from the vector of filter taps h. ratio is an optional rational integer which specifies the input to output sample rate relationship (e.g. 147//160 for converting recorded audio from 48 KHz to 44.1 KHz).

FIRFilter(h, rate[, Nϕ])

Returns a polyphase FIRFilter object from the vector of filter taps h. rate is a floating point number that specifies the input to output sample-rate relationship $\frac{fs_{out}}{fs_{in}}$. Nϕ is an optional parameter which specifies the number of phases created from h. Nϕ defaults to 32.

Filter application

Base.DSP.filt — Function.

filt(f, x[, si])

Apply filter or filter coefficients f along the first dimension of array x. If f is a filter coefficient object, si is an optional array representing the initial filter state (defaults to zeros). If f is a PolynomialRatio, Biquad, or SecondOrderSections, filtering is implemented directly. If f is a ZeroPoleGain object, it is first converted to a SecondOrderSections object. If f is a Vector, it is interpreted as an FIR filter, and a naïve or FFT-based algorithm is selected based on the data and filter length.

Base.DSP.filt! — Function.

filt!(out, f, x[, si])

Same as filt() but writes the result into the out argument, which may alias the input x to modify it in-place.

DSP.Filters.filtfilt — Function.

filtfilt(coef, x)

Filter x in the forward and reverse directions using filter coefficients coef. The initial state of the filter is computed so that its response to a step function is steady state. Before filtering, the data is extrapolated at both ends with an odd-symmetric extension of length 3*(max(length(b), length(a))-1).

Because filtfilt applies the given filter twice, the effective filter order is twice the order of coef. The resulting signal has zero phase distortion.

DSP.Filters.fftfilt — Function.

fftfilt(h, x)

Apply FIR filter taps h along the first dimension of array x using an FFT-based overlap-save algorithm.

DSP.Filters.resample — Function.

resample(x, rate[, coef])

Resample x at rational or arbitrary rate. coef is an optional vector of FIR filter taps. If coef is not provided, the taps will be computed using a Kaiser window.

Internally, resample uses a polyphase FIRFilter object, but performs additional operations to make resampling a signal easier. It compensates for for the FIRFilter's delay (ramp-up), and appends zeros to x. The result is that when the input and output signals are plotted on top of each other, they correlate very well, but one signal will have more samples that the other.

Filter design

Most analog and digital filters are constructed by composing response types, which determine the frequency response of the filter, with design methods, which determine how the filter is constructed.

DSP.Filters.analogfilter — Function.

analogfilter(responsetype, designmethod)

Construct an analog filter. See below for possible response and filter types.

DSP.Filters.digitalfilter — Function.

digitalfilter(responsetype, designmethod)

Construct a digital filter. See below for possible response and filter types.

For some filters, the design method inherently implies a response type. Such filters are documented below.

DSP.Filters.iirnotch — Function.

iirnotch(Wn, bandwidth[; fs])

Second-order digital IIR notch filter at frequency Wn with bandwidth bandwidth. If fs is not specified, Wn is interpreted as a normalized frequency in half-cycles/sample.

Filter response types

DSP.Filters.Lowpass — Type.

Lowpass(Wn[; fs])

Low pass filter with cutoff frequency Wn. If fs is not specified, Wn is interpreted as a normalized frequency in half-cycles/sample.

DSP.Filters.Highpass — Type.

Highpass(Wn[; fs])

High pass filter with cutoff frequency Wn. If fs is not specified, Wn is interpreted as a normalized frequency in half-cycles/sample.

DSP.Filters.Bandpass — Type.

Bandpass(Wn1, Wn2[; fs])

Band pass filter with normalized pass band (Wn1, Wn2). If fs is not specified, Wn1 and Wn2 are interpreted as normalized frequencies in half-cycles/sample.

DSP.Filters.Bandstop — Type.

Bandstop(Wn1, Wn2[; fs])

Band stop filter with normalized stop band (Wn1, Wn2). If fs is not specified, Wn1 and Wn2 are interpreted as normalized frequencies in half-cycles/sample.

Filter design methods

IIR filter design methods

DSP.Filters.Butterworth — Function.

Butterworth(n)

$n$ pole Butterworth filter.

DSP.Filters.Chebyshev1 — Function.

Chebyshev1(n, ripple)

n pole Chebyshev type I filter with ripple dB ripple in the passband.

DSP.Filters.Chebyshev2 — Function.

Chebyshev2(n, ripple)

n pole Chebyshev type II filter with ripple dB ripple in the stopband.

DSP.Filters.Elliptic — Function.

Elliptic(n, rp, rs)

n pole elliptic (Cauer) filter with rp dB ripple in the passband and rs dB attentuation in the stopband.

FIR filter design methods

DSP.Filters.FIRWindow — Type.

FIRWindow(window; scale=true)

FIR filter design using window window, a vector whose length matches the number of taps in the resulting filter.

If scale is true (default), the designed FIR filter is scaled so that the following holds:

For Lowpass and Bandstop filters, the frequency response is unity at 0 (DC).
For Highpass filters, the frequency response is unity at the Nyquist frequency.
For Bandpass filters, the frequency response is unity in the center of the passband.

FIRWindow(; transitionwidth, attenuation=60, scale=true)

Kaiser window FIR filter design. The required number of taps is calculated based on transitionwidth (in half-cycles/sample) and stopband attenuation (in dB). attenuation defaults to 60 dB.

Filter response

DSP.Filters.freqz — Function.

freqz(filter, w = range(0, stop=π, length=250))

Frequency response of a digital filter at normalised frequency or frequencies w in radians/sample.

freqz(filter, hz, fs)

Frequency response of a digital filter at frequency or frequencies hz with sampling rate fs.

DSP.Filters.phasez — Function.

phasez(filter, w = range(0, stop=π, length=250))

Phase response of a digital filter at normalised frequency or frequencies w in radians/sample.

DSP.Filters.impz — Function.

impz(filter, n=100)

Impulse response of a digital filter with n points.

DSP.Filters.stepz — Function.

stepz(filter, n=100)

Step response of a digital filter with n points.

DSP.Filters.freqs — Function.

freqs(filter, w)

Frequency response of an analog filter at normalised frequency or frequencies w in radians/sample.

freqs(filter, hz, fs)

Frequency response of an analog filter at frequency or frequencies hz with sampling rate fs.

Miscellaneous

DSP.Filters.coefb — Function.

coefb(f)

Coefficients of the numerator of a PolynomialRatio object, highest power first, i.e., the b passed to filt()

DSP.Filters.coefa — Function.

coefa(f)

Coefficients of the denominator of a PolynomialRatio object, highest power first, i.e., the a passed to filt()

Examples

Construct a 4th order elliptic lowpass filter with normalized cutoff frequency 0.2, 0.5 dB of passband ripple, and 30 dB attentuation in the stopband and extract the coefficients of the numerator and denominator of the transfer function:

responsetype = Lowpass(0.2)
designmethod = Elliptic(4, 0.5, 30)
tf = convert(PolynomialRatio, digitalfilter(responsetype, designmethod))
numerator_coefs = coefb(tf)
denominator_coefs = coefa(tf)

Filter the data in x, sampled at 1000 Hz, with a 4th order Butterworth bandpass filter between 10 and 40 Hz:

responsetype = Bandpass(10, 40; fs=1000)
designmethod = Butterworth(4)
filt(digitalfilter(responsetype, designmethod), x)

Filter the data in x, sampled at 50 Hz, with a 64 tap Hanning window FIR lowpass filter at 5 Hz:

responsetype = Lowpass(5; fs=50)
designmethod = FIRWindow(hanning(64))
filt(digitalfilter(responsetype, designmethod), x)