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.ZeroPoleGainType
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[m]!)}{(x - \verb!p[1]!) \ldots (x - \verb!p[n]!)}\]

source
DSP.Filters.PolynomialRatioType
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^{m-1} + \ldots + \verb!b[m]!}{\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.

source
DSP.Filters.BiquadType
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}}\]

source
DSP.Filters.SecondOrderSectionsType
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.

source

These filter coefficient objects support the following arithmetic operations: inversion (inv), multiplication (*) for series connection, and integral power (^) for repeated mutlpilcation with itself. For example:

julia> H = PolynomialRatio([1.0], [1.0, 0.3])
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.3*z⁻¹ + 1.0))

julia> inv(H)
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(0.3*z⁻¹ + 1.0), Polynomials.LaurentPolynomial(1.0))

julia> H * H
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))

julia> H^2
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))

julia> H^-2
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0), Polynomials.LaurentPolynomial(1.0))

Stateful filter objects

DSP.Filters.DF2TFilterType
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.

source

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.FIRFilterType
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).

source
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}}$. is an optional parameter which specifies the number of phases created from h. defaults to 32.

source

Filter application

DSP.filtFunction
filt(b, a, x, [si])

Apply filter described by vectors a and b to vector x, with an optional initial filter state vector si (defaults to zeros).

source
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.

source
DSP.filt!Function
filt!(out, b, a, x, [si])

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

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

Same as filt() but writes the result into the out argument. Output array out may not be an alias of x, i.e. filtering may not be done in place.

source
DSP.Filters.filtfiltFunction
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.

source
DSP.Filters.fftfiltFunction
fftfilt(h, x)

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

source
DSP.Filters.tdfiltFunction
tdfilt(h, x)

Apply filter or filter coefficients h along the first dimension of array x using a naïve time-domain algorithm

source
DSP.Filters.tdfilt!Function
tdfilt!(out, h, x)

Like tdfilt, but writes the result into array out. Output array out may not be an alias of x, i.e. filtering may not be done in place.

source
DSP.Filters.resampleFunction
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 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.

source
resample(x::AbstractMatrix, rate::Real, h::Vector = resample_filter(rate); dims)

Resample a matrix x along dimension dims.

source

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. The response type is Lowpass, Highpass, Bandpass or Bandstop and includes the edges of the bands. The design method is Butterworth, Chebyshev1, Chebyshev2, Elliptic, or FIRWindow, and includes any necessary parameters for the method that affect the shape of the response, such as filter order, ripple, and attenuation.

DSP.Filters.analogfilterFunction
analogfilter(responsetype, designmethod)

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

source
DSP.Filters.digitalfilterFunction
digitalfilter(responsetype, designmethod)

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

source

For some filters, the design method is more general or inherently implies a response type; these direct design methods include remez which designs equiripple FIR filters of all types, and iirnotch which designs a 2nd order "biquad" IIR notch filter.

Filter response types

DSP.Filters.LowpassType
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.

source
DSP.Filters.HighpassType
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.

source
DSP.Filters.BandpassType
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.

source
DSP.Filters.BandstopType
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.

source

Filter design methods

IIR filter design methods

DSP.Filters.EllipticFunction
Elliptic(n, rp, rs)

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

source

FIR filter design methods

DSP.Filters.FIRWindowType
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.
source
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.

source

Direct filter design methods

DSP.Filters.remezFunction
remez(numtaps::Integer, band_defs;
      Hz::Real=1.0,
      neg::Bool=false,
      maxiter::Integer=25,
      grid_density::Integer=16)

Calculate the minimax optimal filter using the Remez exchange algorithm [McClellan1973a] [McClellan1973b].

This is the simplified API that accepts just 2 required arguments (numtaps, band_defs). For a scipy compatible version see the 3 arguments version (numtaps, bands, desired).

Calculate the filter-coefficients for the finite impulse response (FIR) filter whose transfer function minimizes the maximum error between the desired gain and the realized gain in the specified frequency bands using the Remez exchange algorithm.

Arguments

  • numtaps::Integer: The desired number of taps in the filter. The number of taps is the number of terms in the filter, or the filter order plus one.
  • bands_defs: A sequence of band definitions. This sequence defines the bands. Each entry is a pair. The pair's first item is a tuple of band edges (low, high). The pair's second item defines the desired response and weight in that band. The weight is optional and defaults to 1.0. Both the desired response and weight may be either scalars or functions. If a function, the function should accept a real frequency and return the real desired response or real weight. Examples:
    • LPF with unity weights. [(0, 0.475) => 1, (0.5, 1.0) => 0]
    • LPF with weight of 2 in the stop band. [(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]
    • BPF with unity weights. [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]
    • Hilbert transformer. [(0.1, 0.95) => 1]; neg=true
    • Differentiator. [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true
  • Hz::Real: The sampling frequency in Hz. Default is 1.
  • neg::Bool: Whether the filter has negative symmetry or not. Default is false. If false, the filter is even-symmetric. If true, the filter is odd-symmetric. neg=true means that h[n]=-h[end+1-n]; neg=false means that h[n]=h[end+1-n].
  • maxiter::Integer: (optional) Maximum number of iterations of the algorithm. Default is 25.
  • grid_density:Integer: (optional) Grid density. The dense grid used in remez is of size (numtaps + 1) * grid_density. Default is 16.

Returns

  • h::Array{Float64,1}: A rank-1 array containing the coefficients of the optimal (in a minimax sense) filter.

J. H. McClellan and T. W. Parks, A unified approach to the design of optimum FIR linear phase digital filters, IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.

J. H. McClellan, T. W. Parks and L. R. Rabiner, A Computer Program for Designing Optimum FIR Linear Phase Digital Filters, IEEE Trans. Audio Electroacoust., vol. AU-21, pp. 506-525, 1973.

Examples

Construct a length 35 filter with a passband at 0.15-0.4 Hz (desired response of 1), and stop bands at 0-0.1 Hz and 0.45-0.5 Hz (desired response of 0). Note: the behavior in the frequency ranges between those bands - the transition bands - is unspecified.

julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);

You can trade-off maximum error achieved for transition bandwidth. The wider the transition bands, the lower the maximum error in the bands specified. Here is a bandpass filter with the same passband, but wider transition bands.

julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);

Here we compute the frequency responses and plot them in dB.

julia> using PyPlot
julia> b = DSP.Filters.PolynomialRatio(bpass, [1.0])
julia> b2 = DSP.Filters.PolynomialRatio(bpass2, [1.0])
julia> f = range(0, stop=0.5, length=1000)
julia> plot(f, 20*log10.(abs.(freqresp(b,f,1.0))))
julia> plot(f, 20*log10.(abs.(freqresp(b2,f,1.0))))
julia> grid()

Examples from the unittests - standard (even) symmetry.

Length 151 LPF (Low Pass Filter).

julia> h = remez(151, [(0, 0.475) => 1, (0.5, 1.0) => 0]; Hz=2.0);

Length 152 LPF. Non-default "weight" input.

julia> h = remez(152, [(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]; Hz=2.0);

Length 51 HPF (High Pass Filter).

julia> h = remez(51, [(0, 0.75) => 0, (0.8, 1.0) => 1]; Hz=2.0);

Length 180 BPF (Band Pass Filter).

julia> h = remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0, maxiter=30);

Examples from the unittests - Odd-symmetric filters - hilbert and differentiators type.

Even length - has a much better approximation since the response is not constrained to 0 at the nyquist frequency. Length 20 Hilbert transformer.

julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);

Length 21 Hilbert transformer.

julia> h = remez(21, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);

Length 200 differentiator.

julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);

Length 201 differentiator.

julia> h = remez(201, [(0.05, 0.95) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);

Inverse sinc filter - custom response function

julia> L = 64; Fs = 4800*L;
julia> passband_response_function = f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sin.(π*f/4800));
julia> h = remez(201, [(    0.0, 2880.0) => (passband_response_function, 1.0),
                (10000.0,   Fs/2) => (0.0, 100.0)]; Hz=Fs);
source
remez(numtaps::Integer, 
      bands::Vector, 
      desired::Vector; 
      weight::Vector=[], 
      Hz::Real=1.0, 
      filter_type::RemezFilterType=filter_type_bandpass,
      maxiter::Integer=25, 
      grid_density::Integer=16)

This is the scipy compatible version that requires 3 arguments (numtaps, bands, desired). For a simplified API, see the 2 argument version (numtaps, band_defs). The filters designed are equivalent, the inputs are just specified in a different way. Below the arguments and examples are described that differ from the simplified API version.

Arguments

  • bands::Vector: A monotonic sequence containing the band edges in Hz. All elements must be non-negative and less than half the sampling frequency as given by Hz.
  • desired::Vector:A sequence half the size of bands containing the desired gain in each of the specified bands.
  • weight::Vector: (optional) A relative weighting to give to each band region. The length of weight has to be half the length of bands.
  • filter_type::RemezFilterType: Default is filter_type_bandpass. The type of filter:
    • filter_type_bandpass : flat response in bands. This is the default.
    • filter_type_differentiator : frequency proportional response in bands. Odd symetric as in filter_type_hilbert case, but with a linear sloping desired response.
    • filter_type_hilbert : filter with odd symmetry, that is, type III (for even order) or type IV (for odd order) linear phase filters.

Examples

Compare the examples with the simplified API and the Scipy API. Each of the following blocks first designs a filter using the simplified (recommended) API, and then designs the same filter using the Scipy-compatible API.

julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);

julia> bpass = remez(35, [0, 0.1, 0.15, 0.4, 0.45, 0.5], [0, 1, 0]);
julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);

julia> bpass2 = remez(35, [0, 0.08, 0.15, 0.4, 0.47, 0.5], [0, 1, 0]);
julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);

julia> h = remez(20, [0.1, 0.95], [1]; filter_type=filter_type_hilbert, Hz=2.0);
julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);

julia> h = remez(200, [0.01, 0.99], [1]; filter_type=filter_type_differentiator, Hz=2.0);
source
DSP.Filters.iirnotchFunction
iirnotch(Wn, bandwidth[; fs])

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

source

Filter response

DSP.Filters.freqrespFunction
H, w = freqresp(filter)

Frequency response H of a filter at (normalized) frequencies w in radians/sample for a digital filter or radians/second for an analog filter chosen as a reasonable default.

source
freqresp(filter::FilterCoefficients{:z}, w)

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

source
freqresp(filter::FilterCoefficients{:s}, w)

Frequency response of analog filter at frequency or frequencies w in radians/second.

source
DSP.Filters.phaserespFunction
phi, w = phaseresp(filter)

Phase response phi of a filter at (normalized) frequencies w in radians/sample for a digital filter or radians/second for an analog filter chosen as a reasonable default.

source
phaseresp(filter, w)

Phase response of a filter at (normalized) frequency or frequencies w in radians/sample for a digital filter or radians/second for an analog filter.

source
DSP.Filters.grpdelayFunction
tau, w = grpdelay(filter)

Group delay tau of a filter at (normalized) frequencies w in radians/sample for a digital filter or radians/second for an analog filter chosen as a reasonable default.

source
grpdelay(fliter, w)

Group delay of a digital 'filter' at normalized frequency or frequencies 'w' in radians/sample.

source

Miscellaneous

DSP.Filters.coefbFunction
coefb(f)

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

source
DSP.Filters.coefaFunction
coefa(f)

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

source

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)