Enable Analog Filter Design (#458)

* Shift frequency normalization from construction of FilterType to filter design

When designing a digital filter using function digitalfilter() the
sampling frequency is not passed to the FilterType definition anymore,
but to function digitalfilter(). That means frequency normalization is
not done during construction of FilterType anymore, but during filter
design inside function digitalfilter().

For this the sampling frequency has to be passed to the internal
functions prewarp(), firprototype() and scalefactor() as well.

The intended side effect is, that an analog filter can now be designed
by specifying analog frequencies directly without any frequency
normalization and without a sampling frequency.

* Update documentation

* Extend filter design tests

Extend test case 'digital IIR' to test digital filter design with
non-normalized frequency values in the filter response definition and
specified sampling frequency.

* Fix test set "freq. scaling"

This is an amendment to commit 54bc9d5
where test sets have been introduced. The code following the current
test set "freq. scaling" is actually a "lowpass" test. It seems that a
corresponding comment on this low pass filter test was already missing
before this commit. And the comment "Freqency scaling" seems to be meant
for all four following test sets "lowpass", "highpass", "bandpass" and
"bandstop".

* Rename test set "freq. scaling" to "analogfilter"

The new test set name relates to the name of the tested function
`analogfilter()` and not to what it is doing.

* Remove unnecessary conversion of filter representation

In the current version the function `analogfilter()` already generates
a `ZeroPoleGain` filter representation. It seems, that in earlier
versions this wasn't the case (see e.g. commit
7fd2102)

* Update analog filter design tests using non-normalized frequency values

For providing new test data I had no MATLAB available and therefore
switched to Octave. Unfortunatelly Octave's elliptic filter design
function `ellip()` generates different filter coefficients than Julia
when designing a 20th order analog filter. Maybe this is a bug in
Octave's 'signal' package. I changed the test to design a Chebyshev
type II filter instead.

---------

Co-authored-by: Christian Gruber <Christian.Gruber-e6t@rub.de>
Co-authored-by: Viral B. Shah <ViralBShah@users.noreply.github.com>

Author: 3 people

docs/src/filters.md

@@ -115,6 +115,18 @@ Bandpass
 Bandstop
 ```
 
+The interpretation of the frequencies `Wn`, `Wn1` and `Wn2` depends on wether an analog
+or a digital filter is designed.
+1. If an analog filter is designed using [`analogfilter`](@ref), the frequencies are
+   interpreted as analog frequencies in radians/second.
+1. If a digital filter is designed using [`digitalfilter`](@ref) and the sampling
+   frequency `fs` is specified, the frequencies of the filter response type are
+   normalized to `fs`. This requires that the sampling frequency and the filter response
+   type use the same frequency unit (Hz, radians/second, ...). If `fs` is not specified,
+   the frequencies of the filter response type are interpreted as normalized frequencies
+   in half-cycles/sample.
+
+
 ### [Filter design methods](@id design-methods)
 
 #### IIR filter design methods
@@ -192,18 +204,18 @@ Filter the data in `x`, sampled at 1000 Hz, with a 4th order
 Butterworth bandpass filter between 10 and 40 Hz:
 
 ```julia
-responsetype = Bandpass(10, 40; fs=1000)
+responsetype = Bandpass(10, 40)
 designmethod = Butterworth(4)
-filt(digitalfilter(responsetype, designmethod), x)
+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:
 
 ```julia
-responsetype = Lowpass(5; fs=50)
+responsetype = Lowpass(5)
 designmethod = FIRWindow(hanning(64))
-filt(digitalfilter(responsetype, designmethod), x)
+filt(digitalfilter(responsetype, designmethod; fs=50), x)
 ```
 
 Estimate a Lowpass Elliptic filter order with a normalized

src/Filters/design.jl

@@ -245,42 +245,36 @@ struct Lowpass{T} <: FilterType
 end
 
 """
-    Lowpass(Wn[; fs])
+    Lowpass(Wn)
 
-Low pass filter with cutoff frequency `Wn`. If `fs` is not
-specified, `Wn` is interpreted as a normalized frequency in
-half-cycles/sample.
+Low pass filter with cutoff frequency `Wn`.
 """
-Lowpass(w::Real; fs::Real=2) = Lowpass{typeof(w/1)}(normalize_freq(w, fs))
+Lowpass(w::Real) = Lowpass{typeof(w/1)}(w)
 
 struct Highpass{T} <: FilterType
     w::T
 end
 
 """
-    Highpass(Wn[; fs])
+    Highpass(Wn)
 
-High pass filter with cutoff frequency `Wn`. If `fs` is not
-specified, `Wn` is interpreted as a normalized frequency in
-half-cycles/sample.
+High pass filter with cutoff frequency `Wn`.
 """
-Highpass(w::Real; fs::Real=2) = Highpass{typeof(w/1)}(normalize_freq(w, fs))
+Highpass(w::Real) = Highpass{typeof(w/1)}(w)
 
 struct Bandpass{T} <: FilterType
     w1::T
     w2::T
 end
 
 """
-    Bandpass(Wn1, Wn2[; fs])
+    Bandpass(Wn1, Wn2)
 
-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.
+Band pass filter with pass band frequencies (`Wn1`, `Wn2`).
 """
-function Bandpass(w1::Real, w2::Real; fs::Real=2)
+function Bandpass(w1::Real, w2::Real)
     w1 < w2 || error("w1 must be less than w2")
-    Bandpass{Base.promote_typeof(w1/1, w2/1)}(normalize_freq(w1, fs), normalize_freq(w2, fs))
+    Bandpass{Base.promote_typeof(w1/1, w2/1)}(w1, w2)
 end
 
 struct Bandstop{T} <: FilterType
@@ -289,15 +283,13 @@ struct Bandstop{T} <: FilterType
 end
 
 """
-    Bandstop(Wn1, Wn2[; fs])
+    Bandstop(Wn1, Wn2)
 
-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.
+Band stop filter with stop band frequencies (`Wn1`, `Wn2`).
 """
-function Bandstop(w1::Real, w2::Real; fs::Real=2)
+function Bandstop(w1::Real, w2::Real)
     w1 < w2 || error("w1 must be less than w2")
-    Bandstop{Base.promote_typeof(w1/1, w2/1)}(normalize_freq(w1, fs), normalize_freq(w2, fs))
+    Bandstop{Base.promote_typeof(w1/1, w2/1)}(w1, w2)
 end
 
 #
@@ -444,20 +436,20 @@ function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K}
 end
 
 # Pre-warp filter frequencies for digital filtering
-prewarp(ftype::Union{Lowpass, Highpass}) = (typeof(ftype))(prewarp(ftype.w))
-prewarp(ftype::Union{Bandpass, Bandstop}) = (typeof(ftype))(prewarp(ftype.w1), prewarp(ftype.w2))
+prewarp(ftype::Union{Lowpass, Highpass}, fs::Real) = (typeof(ftype))(prewarp(normalize_freq(ftype.w, fs)))
+prewarp(ftype::Union{Bandpass, Bandstop}, fs::Real) = (typeof(ftype))(prewarp(normalize_freq(ftype.w1, fs)), prewarp(normalize_freq(ftype.w2, fs)))
 # freq in half-samples per cycle
 prewarp(f::Real) = 4*tan(pi*f/2)
 
 # Digital filter design
 """
-    digitalfilter(responsetype, designmethod)
+    digitalfilter(responsetype, designmethod[; fs])
 
 Construct a digital filter. See below for possible response and
 filter types.
 """
-digitalfilter(ftype::FilterType, proto::FilterCoefficients) =
-    bilinear(transform_prototype(prewarp(ftype), proto), 2)
+digitalfilter(ftype::FilterType, proto::FilterCoefficients; fs::Real=2) =
+    bilinear(transform_prototype(prewarp(ftype, fs), proto), 2)
 
 #
 # Special filter types
@@ -541,61 +533,61 @@ FIRWindow(; transitionwidth::Real=throw(ArgumentError("must specify transitionwi
     FIRWindow(kaiser(kaiserord(transitionwidth, attenuation)...), scale)
 
 # Compute coefficients for FIR prototype with specified order
-function firprototype(n::Integer, ftype::Lowpass)
-    w = ftype.w
+function firprototype(n::Integer, ftype::Lowpass, fs::Real)
+    w = normalize_freq(ftype.w, fs)
 
     [w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
 end
 
-function firprototype(n::Integer, ftype::Bandpass)
-    w1 = ftype.w1
-    w2 = ftype.w2
+function firprototype(n::Integer, ftype::Bandpass, fs::Real)
+    w1 = normalize_freq(ftype.w1, fs)
+    w2 = normalize_freq(ftype.w2, fs)
 
     [w2*sinc(w2*(k-(n-1)/2)) - w1*sinc(w1*(k-(n-1)/2)) for k = 0:(n-1)]
 end
 
-function firprototype(n::Integer, ftype::Highpass)
-    w = ftype.w
+function firprototype(n::Integer, ftype::Highpass, fs::Real)
+    w = normalize_freq(ftype.w, fs)
     isodd(n) || throw(ArgumentError("FIRWindow highpass filters must have an odd number of coefficients"))
 
     out = [-w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
     out[div(n, 2)+1] += 1
     out
 end
 
-function firprototype(n::Integer, ftype::Bandstop)
-    w1 = ftype.w1
-    w2 = ftype.w2
+function firprototype(n::Integer, ftype::Bandstop, fs::Real)
+    w1 = normalize_freq(ftype.w1, fs)
+    w2 = normalize_freq(ftype.w2, fs)
     isodd(n) || throw(ArgumentError("FIRWindow bandstop filters must have an odd number of coefficients"))
 
     out = [w1*sinc(w1*(k-(n-1)/2)) - w2*sinc(w2*(k-(n-1)/2)) for k = 0:(n-1)]
     out[div(n, 2)+1] += 1
     out
 end
 
-scalefactor(coefs::Vector, ::Union{Lowpass, Bandstop}) = sum(coefs)
-function scalefactor(coefs::Vector, ::Highpass)
+scalefactor(coefs::Vector, ::Union{Lowpass, Bandstop}, fs::Real) = sum(coefs)
+function scalefactor(coefs::Vector, ::Highpass, fs::Real)
     c = zero(coefs[1])
     for k = 1:length(coefs)
         c += ifelse(isodd(k), coefs[k], -coefs[k])
     end
     c
 end
-function scalefactor(coefs::Vector, ftype::Bandpass)
+function scalefactor(coefs::Vector, ftype::Bandpass, fs::Real)
     n = length(coefs)
-    freq = middle(ftype.w1, ftype.w2)
+    freq = normalize_freq(middle(ftype.w1, ftype.w2), fs)
     c = zero(coefs[1])
     for k = 0:n-1
         c += coefs[k+1]*cospi(freq*(k-(n-1)/2))
     end
     c
 end
 
-function digitalfilter(ftype::FilterType, proto::FIRWindow)
-    coefs = firprototype(length(proto.window), ftype)
+function digitalfilter(ftype::FilterType, proto::FIRWindow; fs::Real=2)
+    coefs = firprototype(length(proto.window), ftype, fs)
     @assert length(proto.window) == length(coefs)
     out = coefs .* proto.window
-    proto.scale ? rmul!(out, 1/scalefactor(out, ftype)) : out
+    proto.scale ? rmul!(out, 1/scalefactor(out, ftype, fs)) : out
 end
 
 

test/filt_stream.jl

@@ -284,7 +284,7 @@ end
 function test_arbitrary(Th, x, resampleRate, numFilters)
     cutoffFreq      = 0.45
     transitionWidth = 0.05
-    h               = digitalfilter(Lowpass(cutoffFreq, fs=numFilters), FIRWindow(transitionwidth=transitionWidth/numFilters)) .* numFilters
+    h               = digitalfilter(Lowpass(cutoffFreq), FIRWindow(transitionwidth=transitionWidth/numFilters); fs=numFilters) .* numFilters
     h               = convert(Vector{Th}, h)
     myfilt          = FIRFilter(h, resampleRate, numFilters)
     xLen            = length(x)

test/filter_conversion.jl

@@ -133,7 +133,7 @@ using Polynomials.PolyCompat
     1.0000000000000000e+00  -1.9021224191804869e+00   1.0000000000000000e+00   1.0000000000000000e+00  -1.8964983429993663e+00   9.9553672990017417e-01
     1.0000000000000000e+00  -1.9021224191804869e+00   1.0000000000000000e+00   1.0000000000000000e+00  -1.8992956433548462e+00   9.9559721515078736e-01
     ]
-    f = convert(SecondOrderSections, digitalfilter(Bandstop(49.5, 50.5; fs=1000), DSP.Butterworth(2)))
+    f = convert(SecondOrderSections, digitalfilter(Bandstop(49.5, 50.5), DSP.Butterworth(2); fs=1000))
     @test m_sos_butterworth_bs ≈ sosfilter_to_matrix(f)
     @test f.g ≈ 0.995566972017647
 
@@ -349,9 +349,9 @@ end
     bp1 = 0.75
     bp2 = 10.0
     Fsamp = 180
-    responsetype = Bandpass(bp1, bp2; fs = Fsamp)
+    responsetype = Bandpass(bp1, bp2)
     designmethod = Elliptic(11, 0.25, 40)
-    H = digitalfilter(responsetype, designmethod)
+    H = digitalfilter(responsetype, designmethod; fs = Fsamp)
     H′ = ZeroPoleGain(SecondOrderSections(H))
     @test sort(H.p, by=z->(real(z), imag(z))) ≈ sort(H′.p, by=z->(real(z), imag(z)))
     @test sort(H.z, by=z->(real(z), imag(z))) ≈ sort(H′.z, by=z->(real(z), imag(z)))