Merge pull request #45 from fhagemann/derivative

Add `DerivativeFilter`

src/LegendDSP.jl

@@ -27,13 +27,15 @@ using LegendDataTypes
 
 import Adapt
 import RadiationDetectorDSP: fltinstance, rdfilt!, flt_output_length, 
-    flt_input_smpltype, flt_output_smpltype
+    flt_input_smpltype, flt_output_smpltype, _smpllen, smplinfo, SamplingInfo,
+    AbstractRadSigFilterInstance, LinearFiltering, _filterlen, _floattype
 
 using Unitful: RealOrRealQuantity as RealQuantity
 
 include("tailstats.jl")
 include("types.jl")
 include("utils.jl")
+include("derivative.jl")
 include("dsp_decaytime.jl")
 include("dsp_filter_optimization.jl")
 include("dsp_icpc.jl")

src/derivative.jl

@@ -0,0 +1,55 @@
+# This file is a part of LegendDSP.jl, licensed under the MIT License (MIT).
+
+"""
+    struct DerivativeFilter <: AbstractRadIIRFilter
+
+Working example:
+```julia
+using LegendDSP
+using Unitful
+
+signal = rand(100)
+t = range(0u"ms", 20u"ms", 100)
+
+wf = RDWaveform(t, signal)
+flt = DerivativeFilter()
+wf_new = flt(wf)
+```
+Constructors:
+
+* ```$(FUNCTIONNAME)(; fields...)```
+
+Fields:
+
+$(TYPEDFIELDS)
+"""
+Base.@kwdef struct DerivativeFilter{T <: RadiationDetectorDSP.RealQuantity} <: AbstractRadIIRFilter
+    "Filter gain"
+    gain::T = 1
+end
+export DerivativeFilter
+
+struct DerivativeFilterInstance{T <: RadiationDetectorDSP.RealQuantity, TG <: RadiationDetectorDSP.RealQuantity} <: AbstractRadSigFilterInstance{LinearFiltering}
+    gain::TG
+    n_input::Int
+end
+
+function RadiationDetectorDSP.fltinstance(flt::DerivativeFilter{TG}, si::SamplingInfo{T}) where {T, TG}
+    DerivativeFilterInstance{T, TG}(flt.gain, _smpllen(si))
+end
+
+RadiationDetectorDSP.flt_output_smpltype(::DerivativeFilterInstance{T, TG})  where {T,TG} = promote_type(typeof(zero(T) * unit(TG)), TG)
+RadiationDetectorDSP.flt_input_smpltype(::DerivativeFilterInstance{T}) where {T} = T
+RadiationDetectorDSP.flt_output_length(fi::DerivativeFilterInstance) = fi.n_input
+RadiationDetectorDSP.flt_input_length(fi::DerivativeFilterInstance) = fi.n_input
+RadiationDetectorDSP.flt_output_time_axis(::DerivativeFilterInstance, time::AbstractVector) = time
+
+function rdfilt!(y::AbstractVector, fi::DerivativeFilterInstance, x::AbstractVector)
+
+    @assert firstindex(y) == firstindex(x)
+    @assert lastindex(y) == lastindex(x)
+    @inbounds @simd for i in eachindex(y)
+        y[i] = fi.gain * (x[max(i,begin+1)] - x[max(i-1,begin)])
+    end
+    y
+end

test/test_derivative.jl

@@ -0,0 +1,32 @@
+using Test
+using LegendDSP
+using RadiationDetectorSignals: RDWaveform
+using Unitful
+
+@testset "Test DerivativeFilter" begin
+    T = Float64
+    sig = rand(100)
+    wvf = RDWaveform(range(0u"ns", step = 16u"ns", length = length(sig)), sig)
+
+    flt = DerivativeFilter()
+    wvf_out = flt(wvf)
+    @test isapprox(wvf_out.signal, vcat(sig[2] - sig[1], diff(sig)...))
+
+    gain = rand()
+    flt = DerivativeFilter(gain)
+    wvf_out = flt(wvf)
+    @test isapprox(wvf_out.signal, gain * vcat(sig[2] - sig[1], diff(sig)...))
+
+    @testset "Test units and precision types" begin
+        for T1 in (Float16, Float32, Float64)
+            for T2 in (Float16, Float32, Float64)
+                sig = rand(T1, 100)
+                wvf = RDWaveform(range(0u"ns", step = 16u"ns", length = length(sig)), sig)
+                gain = rand(T1) * u"keV"
+                flt = DerivativeFilter(gain)
+                wvf_out = flt(wvf)
+                @test isapprox(wvf_out.signal, gain * vcat(sig[2] - sig[1], diff(sig)...))
+            end
+        end
+    end
+end