Add complex valued bandpass filter

docs/src/filters.md

@@ -119,6 +119,7 @@ bilinear
 Lowpass
 Highpass
 Bandpass
+ComplexBandpass
 Bandstop
 ```
 

src/Filters/Filters.jl

@@ -38,6 +38,7 @@ export FilterType,
         Lowpass,
         Highpass,
         Bandpass,
+        ComplexBandpass,
         Bandstop,
         analogfilter,
         digitalfilter,

src/Filters/design.jl

@@ -234,61 +234,84 @@
 # 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 Nyquist frequency $(fs)"))
+    f
 end
 
 """
     Lowpass(Wn::Real)
 
 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
 
 """
     Highpass(Wn::Real)
 
 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
 
 """
     Bandpass(Wn1::Real, Wn2::Real)
 
 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
 
 """
     Bandstop(Wn1::Real, Wn2::Real)
 
 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 @@
 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,19 +594,31 @@
     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)
 
     [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 @@
 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 @@
     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
 

test/data/digitalfilter_hamming_128_complexbandpass_fc1.txt

@@ -0,0 +1,128 @@
+0.0003244322612916426 + 0.0im
+0.00024117018074241317 + 0.0003319422765306228im
+-0.00010642680276501818 + 0.0003275480187840939im
+-0.00013152109701620738 + 4.2733794890724646e-5im
+0.00011900249938979866 + 8.646037674562085e-5im
+-2.886564280178518e-19 + 0.00041380061803169714im
+-0.00044838738988835745 + 0.00032577250777521207im
+-0.0004651663326831893 - 0.0001511417035131978im
+-6.329013745427158e-5 - 0.00019478701410857128im
+-0.0001324490327176748 + 0.0001823004540165023im
+-0.0006505695595044714 - 9.07640428943599e-19im
+-0.0005219277061327729 - 0.0007183718585141435im
+0.0002452429823247377 - 0.0007547802893086026im
+0.0003184480892798616 - 0.00010347005643476167im
+-0.0002990175633422033 - 0.00021724897638864901im
+2.9773099685853895e-18 - 0.001067023445080161im
+0.001175088353257653 - 0.0008537516628054679im
+0.0012289381776472005 + 0.00039930621937015213im
+0.00016745130899444005 + 0.0005153621369706367im
+0.00034910596341747123 - 0.0004805031363764749im
+0.0017014070539441253 + 4.747427252593109e-18im
+0.0013502875329514353 + 0.001858511347791988im
+-0.0006262938009802742 + 0.0019275341207118755im
+-0.0008015851671894735 + 0.00026045080902800373im
+0.0007412096457101393 + 0.0005385203297762062im
+-7.264329664649744e-18 + 0.0026034273040960584im
+-0.002821707494657874 + 0.002050090496460395im
+-0.0029046951541786647 - 0.000943792667144941im
+-0.00038970635403492494 - 0.0011993928301459067im
+-0.0008003944204631822 + 0.001101648409571519im
+-0.00384526997501156 - 2.1458873619490014e-17im
+-0.0030104645068493494 - 0.004143548919447631im
+0.0013785580633875695 - 0.004242765456728057im
+0.001743474025759286 - 0.0005664890508396567im
+-0.0015945161910667057 - 0.0011584838244030914im
+3.0942851386844633e-17 - 0.005544727998724577im
+0.005955788483457553 - 0.004327133621036445im
+0.006082566786490974 + 0.0019763457525830428im
+0.000810535204434714 + 0.0024945708549897545im
+0.0016553907350729034 - 0.0022784498790698144im
+0.007918228053812622 + 4.4188381102305984e-17im
+0.006180384007964175 + 0.008506568810131199im
+-0.002825553047025619 + 0.00869615809624599im
+-0.0035731668502424025 + 0.0011609922875702257im
+0.003272999669444053 + 0.002377973453998638im
+-6.37301723008596e-17 + 0.011419971168861046im
+-0.012332860590669896 + 0.008960347711082848im
+-0.012691866302459233 - 0.004123837343623776im
+-0.001708538077954521 - 0.005258339515157798im
+-0.0035352813046994826 + 0.004865897271568019im
+-0.0171899576750999 - 9.593009897133222e-17im
+-0.01369264317092155 - 0.018846306503919812im
+0.006418337674197185 - 0.01975361219590808im
+0.008368767275319924 - 0.002719177320940775im
+-0.0079587605978641 - 0.0057823780445617505im
+3.245897315016622e-16 - 0.029082028508867944im
+0.033260866159613885 - 0.024165433783253146im
+0.036790610403702906 + 0.011953993956594534im
+0.005432989869401969 + 0.01672102347869875im
+0.012702228363241048 - 0.017483117468860793im
+0.07307046227651212 + 8.15552523229359e-16im
+0.07457615512753078 + 0.10264527161578987im
+-0.05298447973043177 + 0.16306946099214328im
+-0.18707152579089845 + 0.06078322333380365im
+-0.15913254463367585 - 0.11561656126607897im
+1.9137111664316116e-15 - 0.17146137816013457im
+0.10264527161578987 - 0.07457615512753078im
+0.06949413929677596 + 0.022580014630276066im
+0.0066779566437977 + 0.020552637224586754im
+0.010334160836520934 - 0.014223752138628407im
+0.03868393704616105 + 4.3175835328721346e-16im
+0.024165433783253146 + 0.033260866159613885im
+-0.00898684104013703 + 0.027658652720440248im
+-0.009356084211901972 + 0.00303997604006061im
+0.0071188986476609295 + 0.005172182620085615im
+-2.318196675617288e-16 + 0.0207701770162497im
+-0.018846306503919892 + 0.013692643170921445im
+-0.016348621261741567 - 0.005311989054192185im
+-0.0018586073719735797 - 0.005720205310795779im
+-0.0032498325447541956 + 0.004473010759158416im
+-0.013345017972118517 - 1.4894613692899587e-16im
+-0.008960347711082778 - 0.012332860590669946im
+0.0035289651664501646 - 0.010861037996048024im
+0.003847641872914742 - 0.0012501746285604402im
+-0.003039517269535386 - 0.0022083385609242326im
+1.020542703249778e-16 - 0.009143681734206438im
+0.008506568810131234 - 0.006180384007964127im
+0.007530682388089548 + 0.0024468670339646836im
+0.0008702904121415359 + 0.0026784784740096334im
+0.0015417295757285691 - 0.0021220087142884555im
+0.006395589202401618 + 7.138231713683442e-17im
+0.004327133621036377 + 0.0059557884834576im
+-0.0017134151807925684 + 0.005273349694271168im
+-0.0018744662033011688 + 0.0006090509893754281im
+0.001483087589352909 + 0.0010775262064218595im
+-4.97912213708856e-17 + 0.004461107604052578im
+-0.004143548919447665 + 0.003010464506849303im
+-0.003657069066648828 - 0.0011882537702383586im
+-0.0004207922488040332 - 0.0012950653767152048im
+-0.0007412655348931074 + 0.0010202644804952476im
+-0.00305417722754682 - 3.408821932574686e-17im
+-0.0020500904964603633 - 0.0028217074946578964im
+0.0008045032805854749 - 0.0024760065022612684im
+0.0008713441972107048 - 0.0002831168918719907im
+-0.0006818690704327253 - 0.0004954068782008672im
+2.2620689048956738e-17 - 0.002026729313858864im
+0.0018585113477920033 - 0.0013502875329514143im
+0.0016181342655240977 + 0.0005257636940181548im
+0.00018353586639487782 + 0.0005648653144847465im
+0.000318511317162636 - 0.00043839321840810375im
+0.001292182069720249 + 2.884455126188064e-17im
+0.0008537516628054544 + 0.0011750883532576624im
+-0.0003297283779263003 + 0.0010147996004831823im
+-0.00035151622781795787 + 0.00011421454596354103im
+0.000270888124564075 + 0.00019681174282741785im
+-1.7715536177824818e-17 + 0.0007936229618076265im
+-0.0007183718585141494 + 0.0005219277061327647im
+-0.0006187284188699898 - 0.00020103704990992184im
+-6.96325772697628e-5 - 0.00021430703671424575im
+-0.0001203849952862038 + 0.00016569573100793522im
+-0.0004891048268038243 - 1.091797323285137e-17im
+-0.00032577250777520643 - 0.00044838738988836163im
+0.00012787142325466127 - 0.00039354777422600395im
+0.00013989582825453626 - 4.545491002071199e-5im
+-0.00011187852749218314 - 8.128450819368954e-5im
+7.687917741015551e-18 - 0.00034440436837556115im
+0.00033194227653062604 - 0.00024117018074240875im
+0.0003085534161977859 + 0.00010025508226261925im

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)