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
— TypeZeroPoleGain(z::Vector, p::Vector, k::Number)
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
— TypePolynomialRatio(b::Union{Number, Vector{<:Number}}, a::Union{Number, Vector{<:Number}})
Filter representation in terms of the coefficients of the numerator b
and denominator a
in the z or s domain where b
and a
are vectors ordered from highest power to lowest.
Inputs that are Number
s are treated as one-element Vector
s.
Filter with:
- Transfer function in z domain (zero & negative z powers):
\[H(z) = \frac{\verb!b[1]! + \ldots + \verb!b[m]! z^{-m+1}}{\verb!a[1]! + \ldots + \verb!a[n]! z^{-n+1}}\]
returns PolynomialRatio
object with a[1] = 1
and other specified coefficients divided by a[1]
.
julia> PolynomialRatio([1,1],[1,2])
PolynomialRatio{:z, Float64}(LaurentPolynomial(1.0*z⁻¹ + 1.0), LaurentPolynomial(2.0*z⁻¹ + 1.0))
julia> PolynomialRatio{:z}([1,2,3],[2,3,4])
PolynomialRatio{:z, Float64}(LaurentPolynomial(1.5*z⁻² + 1.0*z⁻¹ + 0.5), LaurentPolynomial(2.0*z⁻² + 1.5*z⁻¹ + 1.0))
- Transfer function in s domain (zero & positive s powers):
\[H(s) = \frac{\verb!b[1]! s^{m-1} + \ldots + \verb!b[m]!}{\verb!a[1]! s^{n-1} + \ldots + \verb!a[n]!}\]
returns PolynomialRatio
object with specified b
and a
coefficients.
julia> PolynomialRatio{:s}([1,2,3],[2,3,4])
PolynomialRatio{:s, Int64}(LaurentPolynomial(3 + 2*s + s²), LaurentPolynomial(4 + 3*s + 2*s²))
DSP.Filters.Biquad
— TypeBiquad(b0::T, b1::T, b2::T, a1::T, a2::T) where T <: Number
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
— TypeSecondOrderSections(biquads::Vector{<:Biquad}, gain::Number)
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 multiplication with itself. For example:
julia> H = PolynomialRatio([1.0], [1.0, 0.3])
PolynomialRatio{:z, Float64}(LaurentPolynomial(1.0), LaurentPolynomial(0.3*z⁻¹ + 1.0))
julia> inv(H)
PolynomialRatio{:z, Float64}(LaurentPolynomial(0.3*z⁻¹ + 1.0), LaurentPolynomial(1.0))
julia> H * H
PolynomialRatio{:z, Float64}(LaurentPolynomial(1.0), LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))
julia> H^2
PolynomialRatio{:z, Float64}(LaurentPolynomial(1.0), LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))
julia> H^-2
PolynomialRatio{:z, Float64}(LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0), LaurentPolynomial(1.0))
Stateful filter objects
DSP.Filters.DF2TFilter
— TypeDF2TFilter(coef::FilterCoefficients{:z})
DF2TFilter(coef::FilterCoefficients{:z}, coldims::Tuple)
DF2TFilter(coef::FilterCoefficients{:z}, sitype::Type, coldims::Tuple = ())
DF2TFilter(coef::FilterCoefficients{:z}, si)
Construct a stateful direct form II transposed filter with coefficients coef
.
The initial filter state defaults to zeros (of a type derived from coef
) suitable for vector input. Another element type of the state can be specified with sitype
.
To allow column-wise filtering of higher-dimensional input, the size of the extra dimensions have to be given in coldims
. To e.g. column-wise filter an input with size (L, N1, N2)
, set coldims
to (N1, N2)
.
Alternatively, an array representing the initial filter state can be passed as si
.
If coef
is a PolynomialRatio
, Biquad
, or SecondOrderSections
, filtering is implemented directly. If coef
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
— TypeFIRFilter(h::Vector[, ratio::Union{Integer,Rational}])
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::Vector, rate::AbstractFloat[, Nϕ::Integer=32])
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
— Functionfilt(b::Union{AbstractVector,Number},
a::Union{AbstractVector,Number},
x::AbstractArray)
Apply filter described by vectors a
and b
to vector x
.
Inputs that are Number
s are treated as one-element Vector
s.
filt(f::FilterCoefficients{:z}, x::AbstractArray)
Apply filter or filter coefficients f
along the first dimension of array x
. 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.
filt(f::DF2TFilter{<:FilterCoefficients{:z},<:Array{T}}, x::AbstractArray{V}) where {T,V}
Apply the stateful filter f
on x
.
The output array has eltype promote_type(T, V)
, where T
is the eltype of the filter state.
For more control over the output type, provide a preallocated output array out
to filt!(out, f, x)
.
DSP.filt!
— Functionfilt!(out, b, a, x)
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)
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
— Functionfiltfilt(coef::FilterCoefficients, x::AbstractArray)
filtfilt(b::AbstractVector, x::AbstractArray)
filtfilt(b::AbstractVector, a::AbstractVector, x::AbstractArray)
Filter x
in the forward and reverse directions using either a FilterCoefficients
object coef
, or the coefficients b
and optionally a
as in filt
. 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 min(3*(max(length(b), length(a))-1), size(x, 1) - 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
— Functionfftfilt(h::AbstractVector{<:Real}, x::AbstractArray{<:Real})
Apply FIR filter taps h
along the first dimension of array x
using an FFT-based overlap-save algorithm.
DSP.Filters.fftfilt!
— Functionfftfilt!(out::AbstractArray, h::AbstractVector, x::AbstractArray)
Like fftfilt
but writes result into out array.
DSP.Filters.tdfilt
— Functiontdfilt(h::AbstractVector, x::AbstractArray)
Apply filter or filter coefficients h
along the first dimension of array x
using a naïve time-domain algorithm
DSP.Filters.tdfilt!
— Functiontdfilt!(out::AbstractArray, h::AbstractVector, x::AbstractArray)
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
— Functionresample(x::AbstractVector, rate::Real[, coef::Vector])
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 than the other.
resample(x::AbstractArray, rate::Real, h::Vector = resample_filter(rate); dims)
Resample an array 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.
Filter order estimation methods are available in buttord
, cheb1ord
, cheb2ord
, and ellipord
if the corner frequencies for different IIR filter types are known. remezord
can be used for an initial FIR filter order estimate.
DSP.Filters.analogfilter
— Functionanalogfilter(responsetype::FilterType, designmethod::FilterCoefficients)
Construct an analog filter. See below for possible response and filter types.
DSP.Filters.digitalfilter
— Functiondigitalfilter(responsetype::FilterType, designmethod::FilterCoefficients[; fs::Real])
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.
For a more general application of creating a digital filter from s-domain representation of an analog filter, one can use bilinear
:
DSP.Filters.bilinear
— Functionbilinear(f::FilterCoefficients{:s}, fs::Real)
Calculate the digital filter (z-domain) ZPK representation of an analog filter defined in s-domain using bilinear transform with sampling frequency fs
. The s-domain representation is first converted to a ZPK representation in s-domain and then transformed to z-domain using bilinear transform.
bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K}
Calculate the digital filter (z-domain) ZPK representation of an analog filter defined as a ZPK representation in s-domain using bilinear transform with sampling frequency fs
.
Input s-domain representation must be a ZeroPoleGain{:s, Z, P, K}
object:
\[H(s) = f.k\frac{(s - \verb!f.z[1]!) \ldots (s - \verb!f.z[m]!)}{(s - \verb!f.p[1]!) \ldots (s - \verb!f.p[n]!)}\]
Output z-domain representation is a ZeroPoleGain{:z, Z, P, K}
object:
\[H(z) = K\frac{(z - \verb!Z[1]!) \ldots (z - \verb!Z[m]!)}{(z - \verb!P[1]!) \ldots (z - \verb!P[n]!)}\]
where Z, P, K
are calculated as:
\[Z[i] = \frac{(2 + \verb!f.z[i]!/\verb!fs!)}{(2 - \verb!f.z[i]!/\verb!fs!)} \quad \text{for } i = 1, \ldots, m\]
\[P[i] = \frac{(2 + \verb!f.p[i]!/\verb!fs!)}{(2 - \verb!f.p[i]!/\verb!fs!)} \quad \text{for } i = 1, \ldots, n\]
\[K = f.k \ \mathcal{Re} \left[ \frac{\prod_{i=1}^m (2*fs - f.z[i])}{\prod_{i=1}^n (2*fs - f.p[i])} \right]\]
Here, m
and n
are respectively the numbers of zeros and poles in the s-domain representation. If m < n
, then additional n-m
zeros are added at z = -1
.
Filter response types
DSP.Filters.Lowpass
— TypeLowpass(Wn::Real)
Low pass filter with cutoff frequency Wn
.
DSP.Filters.Highpass
— TypeHighpass(Wn::Real)
High pass filter with cutoff frequency Wn
.
DSP.Filters.Bandpass
— TypeBandpass(Wn1::Real, Wn2::Real)
Band pass filter with pass band frequencies (Wn1
, Wn2
).
DSP.Filters.ComplexBandpass
— TypeComplexBandpass(Wn1, Wn2)
Complex band pass filter with pass band frequencies (Wn1
, Wn2
).
DSP.Filters.Bandstop
— TypeBandstop(Wn1::Real, Wn2::Real)
Band stop filter with stop band frequencies (Wn1
, Wn2
).
The interpretation of the frequencies Wn
, Wn1
and Wn2
depends on wether an analog or a digital filter is designed.
- If an analog filter is designed using
analogfilter
, the frequencies are interpreted as analog frequencies in radians/second. - If a digital filter is designed using
digitalfilter
and the sampling frequencyfs
is specified, the frequencies of the filter response type are normalized tofs
. This requires that the sampling frequency and the filter response type use the same frequency unit (Hz, radians/second, ...). Iffs
is not specified, the frequencies of the filter response type are interpreted as normalized frequencies in half-cycles/sample.
Filter design methods
IIR filter design methods
DSP.Filters.Butterworth
— FunctionButterworth(n::Integer)
$n$ pole Butterworth filter.
DSP.Filters.Chebyshev1
— FunctionChebyshev1(n::Integer, ripple::Real)
n
pole Chebyshev type I filter with ripple
dB ripple in the passband.
DSP.Filters.Chebyshev2
— FunctionChebyshev2(n::Integer, ripple::Real)
n
pole Chebyshev type II filter with ripple
dB ripple in the stopband.
DSP.Filters.Elliptic
— FunctionElliptic(n::Integer, rp::Real, rs::Real)
n
pole elliptic (Cauer) filter with rp
dB ripple in the passband and rs
dB attentuation in the stopband.
Filter order estimation methods
IIR filter order estimation methods
DSP.Filters.buttord
— Function(N, ωn) = buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain=:z)
Butterworth order estimate for bandpass and bandstop filter types. Wp
and Ws
are 2-element pass and stopband frequency edges, with no more than Rp
dB passband ripple and at least Rs
dB stopband attenuation. Based on the ordering of passband and bandstop edges, the Bandstop or Bandpass filter type is inferred. N
is an integer indicating the lowest estimated filter order, with ωn
specifying the cutoff or "-3 dB" frequencies.
If a domain of :s
is specified, the passband and stopband frequencies are interpreted as radians/second, giving an order and natural frequencies for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
(N, ωn) = buttord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain=:z)
LPF/HPF Butterworth filter order and -3 dB frequency approximation. Wp
and Ws
are the passband and stopband frequencies, whereas Rp
and Rs
are the passband and stopband ripple attenuations in dB. If the passband is greater than stopband, the filter type is inferred to be for estimating the order of a highpass filter. N
specifies the lowest possible integer filter order, whereas ωn
is the cutoff or "-3 dB" frequency.
If a domain of :s
is specified, the passband and stopband edges are interpreted as radians/second, giving an order and natural frequency result for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
DSP.Filters.cheb1ord
— Function(N, ωn) = cheb1ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
(N, ωn) = cheb1ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Integer and natural frequency order estimation for Chebyshev Type I Filters. Wp
and Ws
indicate the passband and stopband frequency edges, and Rp
and Rs
indicate the maximum loss in the passband and the minimum attenuation in the stopband, (in dB.)
Based on the ordering of passband and stopband edges, the Lowpass or Highpass filter type is inferred. N
indicates the smallest integer filter order that achieves the desired specifications, and ωn
contains the natural frequency of the filter, (in this case, simply the passband edge.)
If a domain of :s
is specified, the passband and stopband edges are interpreted as radians/second, giving the order and natural frequencies for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
Wp
and Ws
can also be 2-element frequency edges for Bandpass/Bandstop cases.
N
is an integer indicating the lowest estimated filter order, with ωn
specifying the cutoff or "-3 dB" frequencies.
DSP.Filters.cheb2ord
— Function(N, ωn) = cheb2ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
Integer and natural frequency order estimation for Chebyshev Type II (inverse) Filters. Wp
and Ws
are the frequency edges for Bandpass/Bandstop cases, with Rp
and Rs
representing the ripple maximum loss in the passband and minimum ripple attenuation in the stopband in dB.
Based on the ordering of the passband and stopband edges, the Lowpass or Highpass filter type is inferred. N
is an integer indicating the lowest estimated filter order, with ωn
specifying the "-3 dB" cutoff frequency.
If a domain of :s
is specified, the passband and stopband frequencies are interpreted as radians/second, giving the order and natural frequencies for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
(N, ωn) = cheb2ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Integer and natural frequency order estimation for Chebyshev Type II (inverse) Filters. Wp
and Ws
are 2-element frequency edges for Bandpass/Bandstop cases, with Rp
and Rs
representing the ripple maximum loss in the passband and minimum ripple attenuation in the stopband in dB.
Based on the ordering of the passband and bandstop edges, the Bandstop or Bandpass filter type is inferred. N
is an integer indicating the lowest estimated filter order, with ωn
specifying the "-3 dB" cutoff frequency.
If a domain of :s
is specified, the passband and stopband frequencies are interpreted as radians/second, giving the order and natural frequencies for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
DSP.Filters.ellipord
— Function(N, ωn) = ellipord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
(N, ωn) = ellipord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Integer and natural frequency order estimation for Elliptic (Cauer) Filters. Wp
and Ws
indicate the passband and stopband frequency edges, and Rp
and Rs
indicate the maximum loss in the passband and the minimum attenuation in the stopband, (in dB.)
Based on the ordering of passband and stopband edges, the Lowpass or Highpass filter type is inferred. N
indicates the smallest integer filter order that achieves the desired specifications, and ωn
contains the natural frequency of the filter, (in this case, simply the passband edge.)
If a domain of :s
is specified, the passband and stopband edges are interpreted as radians/second, giving the order and natural frequencies for an analog filter. The default domain is :z
, interpreting the input frequencies as normalized from 0 to 1, where 1 corresponds to π radians/sample.
Wp
and Ws
can also be 2-element frequency edges for Bandpass/Bandstop cases.
N
is an integer indicating the lowest estimated filter order, with ωn
specifying the cutoff or "-3 dB" frequencies.
FIR filter order estimation methods
DSP.Filters.remezord
— FunctionN = remezord(Wp::Real, Ws::Real, Rp::Real, Rs::Real)
Order estimation for lowpass digital filter cases based on the equations and coefficients in [Rabiner]. The original equation returned the minimum filter length, whereas this implementation returns the order (N=L-1). Wp
and Ws
are the normalized passband and stopband frequencies, with Rp
indicating the passband ripple and Rs
is the stopband attenuation (linear.)
NOTE: The value of N is an initial approximation. If under-estimated, the order should be increased until the design specifications are met.
FIR filter design methods
DSP.Filters.FIRWindow
— TypeFIRWindow(window::Vector; 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::Real, attenuation::Real=60, scale::Bool=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
— Functionremez(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.band_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 inremez
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.
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 = PolynomialRatio(bpass, [1.0])
julia> b2 = 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 byHz
.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 ofweight
has to be half the length ofbands
.filter_type::RemezFilterType
: Default isfilter_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 infilter_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
— Functioniirnotch(Wn::Real, bandwidth::Real[; fs=2])
Second-order digital IIR notch filter [Orfanidis] 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
— FunctionH, w = freqresp(filter::FilterCoefficients)
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
— Functionphi, w = phaseresp(filter::FilterCoefficients)
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::FilterCoefficients, 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
— Functiontau, w = grpdelay(filter::FilterCoefficients)
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(filter::FilterCoefficients{:z}, w)
Group delay of a digital filter
at normalized frequency or frequencies w
in radians/sample.
DSP.Filters.impresp
— Functionimpresp(filter::FilterCoefficients{:z}, n=100)
Impulse response of a digital filter
with n
points.
DSP.Filters.stepresp
— Functionstepresp(filter::FilterCoefficients{:z}, n=100)
Step response of a digital filter
with n
points.
Miscellaneous
DSP.Filters.coefb
— Functioncoefb(f::PolynomialRatio)
Coefficients of the numerator of a PolynomialRatio
object, highest power first, i.e., the b
passed to filt()
DSP.Filters.coefa
— Functioncoefa(f::PolynomialRatio)
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)
designmethod = Butterworth(4)
filt(digitalfilter(responsetype, designmethod; fs=1000), 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)
designmethod = FIRWindow(hanning(64))
filt(digitalfilter(responsetype, designmethod; fs=50), x)
Estimate a Lowpass Elliptic filter order with a normalized passband cutoff frequency of 0.2, a stopband cutoff frequency of 0.4, 3 dB of passband ripple, and 40 dB attenuation in the stopband:
(N, ωn) = ellipord(0.2, 0.4, 3, 40)
- RabinerHerrmann, O., Lawrence R. Rabiner, and D. S. K. Chan. "Practical design rules for optimum finite impulse response lowpass digital filters." Bell System Technical Journal 52.6 (1973): 769-799.
- McClellan1973aJ. 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.
- McClellan1973bJ. 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.
- OrfanidisOrfanidis, S. J. (1996). Introduction to signal processing. Englewood Cliffs, N.J: Prentice Hall, p. 370.