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

`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^{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.

`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`

.

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

`DSP.filt`

— Function`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).

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

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

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

`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.fftfilt!`

— Function`fftfilt!(out, h, x)`

Like `fftfilt`

but writes result into out array.

`DSP.Filters.tdfilt`

— Function`tdfilt(h, x)`

Apply filter or filter coefficients `h`

along the first dimension of array `x`

using a naïve time-domain algorithm

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

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

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

Resample a matrix `x`

along dimension `dims`

.

## 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.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 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.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:

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

### Direct filter design methods

`DSP.Filters.remez`

— Function```
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`

- LPF with unity weights.
`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);
```

```
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);
```

`DSP.Filters.iirnotch`

— Function`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.

## Filter response

`DSP.Filters.freqresp`

— Function`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.

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

Frequency response of digital `filter`

at normalized frequency or frequencies `w`

in radians/sample.

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

Frequency response of analog `filter`

at frequency or frequencies `w`

in radians/second.

`DSP.Filters.phaseresp`

— Function`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.

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

`DSP.Filters.grpdelay`

— Function`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.

`grpdelay(fliter, w)`

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

`DSP.Filters.impresp`

— Function`impresp(filter, n=100)`

Impulse response of a digital `filter`

with `n`

points.

`DSP.Filters.stepresp`

— Function`stepresp(filter, n=100)`

Step response of a digital `filter`

with `n`

points.

## 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)
```

- McClellan1973a
- McClellan1973b
- OrfandisOrfanidis, S. J. (1996). Introduction to signal processing. Englewood Cliffs, N.J: Prentice Hall, p. 370.