diff --git a/docs/src/filters.md b/docs/src/filters.md index 7230aa800..8ae6f1510 100644 --- a/docs/src/filters.md +++ b/docs/src/filters.md @@ -119,6 +119,7 @@ bilinear Lowpass Highpass Bandpass +ComplexBandpass Bandstop ``` diff --git a/src/Filters/Filters.jl b/src/Filters/Filters.jl index 96b295ff9..557f4581b 100644 --- a/src/Filters/Filters.jl +++ b/src/Filters/Filters.jl @@ -37,6 +37,7 @@ export FilterType, Lowpass, Highpass, Bandpass, + ComplexBandpass, Bandstop, analogfilter, digitalfilter, diff --git a/src/Filters/design.jl b/src/Filters/design.jl index 69d06d415..df485b3fe 100644 --- a/src/Filters/design.jl +++ b/src/Filters/design.jl @@ -234,13 +234,14 @@ Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs) # returns frequency in half-cycles per sample ∈ (0, 1) function normalize_freq(w::Real, fs::Real) w <= 0 && throw(DomainError(w, "frequencies must be positive")) - f = 2*w/fs + f = 2 * w / fs f >= 1 && throw(DomainError(f, "frequencies must be less than the Nyquist frequency $(fs/2)")) f end - -struct Lowpass{T} <: FilterType - w::T +function normalize_complex_freq(w::Real, fs::Real) + f = 2 * w / fs + f >= 2 && throw(DomainError(f, "frequencies must be less than the sampling frequency $(fs)")) + f end """ @@ -248,10 +249,10 @@ end Low pass filter with cutoff frequency `Wn`. """ -Lowpass(w::Real) = Lowpass{typeof(w/1)}(w) - -struct Highpass{T} <: FilterType +struct Lowpass{T<:Real} <: FilterType w::T + Lowpass{T}(w::Real) where {T<:Real} = new{T}(w) + Lowpass(w::Real) = Lowpass{typeof(w / 1)}(w) end """ @@ -259,11 +260,10 @@ end High pass filter with cutoff frequency `Wn`. """ -Highpass(w::Real) = Highpass{typeof(w/1)}(w) - -struct Bandpass{T} <: FilterType - w1::T - w2::T +struct Highpass{T<:Real} <: FilterType + w::T + Highpass{T}(w::Real) where {T<:Real} = new{T}(w) + Highpass(w::Real) = Highpass{typeof(w / 1)}(w) end """ @@ -271,14 +271,31 @@ end Band pass filter with pass band frequencies (`Wn1`, `Wn2`). """ -function Bandpass(w1::Real, w2::Real) - w1 < w2 || throw(ArgumentError("w1 must be less than w2")) - Bandpass{Base.promote_typeof(w1/1, w2/1)}(w1, w2) +struct Bandpass{T<:Real} <: FilterType + w1::T + w2::T + function Bandpass{T}(w1::Real, w2::Real) where {T<:Real} + w1 < w2 || throw(ArgumentError("w1 must be less than w2")) + new{T}(w1, w2) + end + Bandpass(w1::T, w2::V) where {T<:Real,V<:Real} = + Bandpass{typeof(one(promote_type(T, V)) / 1)}(w1, w2) end -struct Bandstop{T} <: FilterType +""" + ComplexBandpass(Wn1, Wn2) + +Complex band pass filter with pass band frequencies (`Wn1`, `Wn2`). +""" +struct ComplexBandpass{T<:Real} <: FilterType w1::T w2::T + function ComplexBandpass{T}(w1::Real, w2::Real) where {T<:Real} + w1 < w2 || throw(ArgumentError("w1 must be less than w2")) + new{T}(w1, w2) + end + ComplexBandpass(w1::T, w2::V) where {T,V} = + ComplexBandpass{typeof(one(promote_type(T, V)) / 1)}(w1, w2) end """ @@ -286,9 +303,15 @@ end Band stop filter with stop band frequencies (`Wn1`, `Wn2`). """ -function Bandstop(w1::Real, w2::Real) - w1 < w2 || throw(ArgumentError("w1 must be less than w2")) - Bandstop{Base.promote_typeof(w1/1, w2/1)}(w1, w2) +struct Bandstop{T<:Real} <: FilterType + w1::T + w2::T + function Bandstop{T}(w1::Real, w2::Real) where {T<:Real} + w1 < w2 || throw(ArgumentError("w1 must be less than w2")) + new{T}(w1, w2) + end + Bandstop(w1::T, w2::V) where {T,V} = + Bandstop{typeof(one(promote_type(T, V)) / 1)}(w1, w2) end # @@ -472,8 +495,10 @@ 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}, 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))) +prewarp(ftype::Lowpass, fs::Real) = Lowpass(prewarp(normalize_freq(ftype.w, fs))) +prewarp(ftype::Highpass, fs::Real) = Highpass(prewarp(normalize_freq(ftype.w, fs))) +prewarp(ftype::Bandpass, fs::Real) = Bandpass(prewarp(normalize_freq(ftype.w1, fs)), prewarp(normalize_freq(ftype.w2, fs))) +prewarp(ftype::Bandstop, fs::Real) = Bandstop(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) @@ -569,12 +594,15 @@ 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, fs::Real) +function _firprototype(n::Integer, ftype::Lowpass, fs::Real, ::Type{T}) where {T<:Number} w = normalize_freq(ftype.w, fs) - [w*sinc(w*(k-(n+1)/2)) for k = 1:n] + promote_type(typeof(w), T)[w*sinc(w*(k-(n+1)/2)) for k = 1:n] end +firprototype(n::Integer, ftype::Lowpass, fs::Real) = + _firprototype(n, ftype, fs, typeof(fs)) + function firprototype(n::Integer, ftype::Bandpass, fs::Real) w1 = normalize_freq(ftype.w1, fs) w2 = normalize_freq(ftype.w2, fs) @@ -582,6 +610,15 @@ function firprototype(n::Integer, ftype::Bandpass, fs::Real) [w2*sinc(w2*(k-(n+1)/2)) - w1*sinc(w1*(k-(n+1)/2)) for k = 1:n] end +function firprototype(n::Integer, ftype::ComplexBandpass, fs::T) where {T<:Real} + w1 = normalize_complex_freq(ftype.w1, fs) + w2 = normalize_complex_freq(ftype.w2, fs) + w_center = (w2 + w1) / 2 + w_cutoff = (w2 - w1) / 2 + lp = Lowpass(w_cutoff) + _firprototype(n, lp, 2, Complex{T}) .*= cispi.(w_center * (0:(n-1))) +end + 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")) @@ -602,14 +639,14 @@ function firprototype(n::Integer, ftype::Bandstop, fs::Real) end scalefactor(coefs::Vector, ::Union{Lowpass, Bandstop}, fs::Real) = sum(coefs) -function scalefactor(coefs::Vector{T}, ::Highpass, fs::Real) where T +function scalefactor(coefs::Vector{T}, ::Highpass, fs::Real) where {T<:Number} c = zero(T) for k = 1:length(coefs) c += ifelse(isodd(k), coefs[k], -coefs[k]) end c end -function scalefactor(coefs::Vector{T}, ftype::Bandpass, fs::Real) where T +function scalefactor(coefs::Vector{T}, ftype::Bandpass, fs::Real) where {T<:Number} n = length(coefs) freq = normalize_freq(middle(ftype.w1, ftype.w2), fs) c = zero(T) @@ -618,11 +655,20 @@ function scalefactor(coefs::Vector{T}, ftype::Bandpass, fs::Real) where T end c end +function scalefactor(coefs::Vector{T}, ftype::ComplexBandpass, fs::Real) where T<:Number + n = length(coefs) + freq = normalize_complex_freq(middle(ftype.w1, ftype.w2), fs) + c = zero(T) + for k = 1:n + c = muladd(coefs[k], cispi(-freq * (k - (n + 1) / 2)), c) + end + return abs(c) +end 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 + out = (coefs .*= proto.window) proto.scale ? rmul!(out, 1/scalefactor(out, ftype, fs)) : out end diff --git a/test/filter_design.jl b/test/filter_design.jl index 9918adbca..e07d697c6 100644 --- a/test/filter_design.jl +++ b/test/filter_design.jl @@ -1012,6 +1012,36 @@ end winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt"),'\t') @test winfirtaps_jl ≈ winfirtaps_scipy + @test_throws ArgumentError ComplexBandpass(2, 1) + + winfirtaps_jl = digitalfilter(ComplexBandpass(0.1, 0.2), FIRWindow(hamming(128))) + winfirfreq_db = 20log10.(abs.(fft([winfirtaps_jl; zeros(400)]))) + f = range(0; step=2/length(winfirfreq_db), length=length(winfirfreq_db)) + # above -6dB in the whole passband + @test minimum(winfirfreq_db[0.1 .< f .< 0.2]) > -6.02 + # close to 0dB in the passband a bit away from the edges + @test all(-0.05 .< winfirfreq_db[0.125 .< f .< 0.175] .< 0.05) + # exactly unity in the middle of the passband + @test abs(sum(winfirtaps_jl .* cispi.(-0.15 * axes(winfirtaps_jl,1)))) ≈ 1 + # below -6dB outside passband + @test maximum(winfirfreq_db[.!(0.1 .< f .< 0.2)]) < -6.02 + # reasonable attenuation outside passband a bit away from the edges + @test maximum(winfirfreq_db[.!(0.07 .< f .< 0.23)]) < -50 + + winfirtaps_jl = digitalfilter(ComplexBandpass(30, 40), FIRWindow(hamming(511)); fs=100) + winfirfreq_db = 20log10.(abs.(fft([winfirtaps_jl; zeros(400)]))) + f = range(0; step=100/length(winfirfreq_db), length=length(winfirfreq_db)) + # above -6dB in the whole passband + @test minimum(winfirfreq_db[30 .< f .< 40]) > -6.02 + # close to 0dB in the passband a bit away from the edges + @test all(-0.01 .< winfirfreq_db[31 .< f .< 39] .< 0.01) + # exactly unity in the middle of the passband + @test abs(sum(winfirtaps_jl .* cispi.(-2*35/100 * axes(winfirtaps_jl,1)))) ≈ 1 + # below -6dB outside passband + @test maximum(winfirfreq_db[.!(30 .< f .< 40)]) < -6.02 + # reasonable attenuation outside passband a bit away from the edges + @test maximum(winfirfreq_db[.!(28 .< f .< 42)]) < -60 + @test_throws ArgumentError digitalfilter(Bandstop(0.1, 0.2),FIRWindow(hamming(128), scale=false); fs=1) winfirtaps_jl = digitalfilter(Bandstop(0.1, 0.2),FIRWindow(hamming(129), scale=false); fs=1)