Merge pull request #47 from legend-exp/waveform-flt-time

Add `TimeAxis` filter

src/LegendDSP.jl

@@ -52,5 +52,6 @@ include("dsp_ml_routines.jl")
 include("multi_intersect.jl")
 include("moving_window_multi.jl")
 include("alternative_filters.jl")
+include("timeaxis.jl")
 
 end # module

src/timeaxis.jl

@@ -0,0 +1,57 @@
+# This file is a part of LegendDSP.jl, licensed under the MIT License (MIT).
+
+"""
+    struct TimeAxisFilter <: AbstractRadIIRFilter
+
+Updates the time range step to `period` and shifts the time axis by `offset`.
+
+Working example:
+```julia
+using LegendDSP
+using RadiationDetectorSignals
+using Unitful
+
+signal = rand(100)
+t = range(0u"ms", 20u"ms", 100)
+
+wf = RDWaveform(t, signal)
+flt = TimeAxisFilter(4u"ns")
+wf_new = flt(wf)
+```
+Constructors:
+
+* ```$(FUNCTIONNAME)(; fields...)```
+
+Fields:
+
+$(TYPEDFIELDS)
+"""
+Base.@kwdef struct TimeAxisFilter{T <: RadiationDetectorDSP.RealQuantity} <: AbstractRadIIRFilter
+    "Filter period"
+    period::T
+    "Filter offset"
+    offset::T = false * period
+    function TimeAxisFilter(period::T, offset::T=false * period) where T
+        new{T}(period, offset)
+    end
+end
+export TimeAxisFilter
+
+struct TimeAxisFilterInstance{T <: RadiationDetectorDSP.RealQuantity, TT <: RadiationDetectorDSP.RealQuantity} <: AbstractRadSigFilterInstance{LinearFiltering}
+    period::TT
+    offset::TT
+    n_input::Int
+end
+
+function RadiationDetectorDSP.fltinstance(flt::TimeAxisFilter{TT}, si::SamplingInfo{T}) where {T <: RadiationDetectorDSP.RealQuantity, TT <: RadiationDetectorDSP.RealQuantity}
+    TimeAxisFilterInstance{T,TT}(flt.period, flt.offset, _smpllen(si))
+end
+
+RadiationDetectorDSP.flt_output_smpltype(fi::TimeAxisFilterInstance{T})  where {T} = flt_input_smpltype(fi)
+RadiationDetectorDSP.flt_input_smpltype(::TimeAxisFilterInstance{T}) where {T} = T
+RadiationDetectorDSP.flt_output_length(fi::TimeAxisFilterInstance) = fi.n_input
+RadiationDetectorDSP.flt_input_length(fi::TimeAxisFilterInstance) = fi.n_input
+RadiationDetectorDSP.flt_output_time_axis(fi::TimeAxisFilterInstance, time::AbstractVector) = throw(ArgumentError("Non-range time axis not supported"))
+RadiationDetectorDSP.flt_output_time_axis(fi::TimeAxisFilterInstance, time::AbstractRange) = range(first(time) + fi.offset, step=fi.period, length=length(time))
+
+rdfilt!(y::AbstractVector, fi::TimeAxisFilterInstance, x::AbstractVector) = y .= x

test/runtests.jl

@@ -5,8 +5,8 @@ import Test
 Test.@testset "Package LegendDSP" begin
     include("test_aqua.jl")
     include("test_haar_filter.jl")
-    # include("test_SOMETHING.jl")
     include("test_multiintersect.jl")
+    include("test_timeaxis.jl")
     include("test_docs.jl")
     include("test_alternative_filters.jl")
     isempty(Test.detect_ambiguities(LegendDSP))

test/test_timeaxis.jl

@@ -0,0 +1,56 @@
+# This file is a part of LegendDSP.jl, licensed under the MIT License (MIT).
+
+using Test
+using LegendDSP
+using RadiationDetectorSignals
+using Unitful
+
+@testset "TimeAxisFilter" begin
+    @testset "Single waveform" begin
+        # create random waveform with random start point
+        old_offset = rand()u"ns"
+        old_step   = rand()u"ns"
+        t = range(old_offset, step = old_step, length = 100)
+        signal = rand(100)
+        wf = RDWaveform(t, signal)
+        @test first(wf.time) == old_offset
+        @test step(wf.time)  == old_step
+
+        # overwrite step and add offset (using random values to test flexibility)
+        new_offset = rand()u"ns"
+        new_step   = rand()u"ns"
+        flt = TimeAxisFilter(new_step, new_offset)
+        wf_new = flt(wf)
+        @test first(wf_new.time) == old_offset + new_offset
+        @test step(wf_new.time)  == new_step
+
+        # explicitly set the start point to zero and use a step of 4ns
+        fixed_step = 4.0u"ns"
+        flt = TimeAxisFilter(fixed_step, -first(wf.time))
+        wf_to0 = flt(wf)
+        @test iszero(first(wf_to0.time))
+        @test step(wf_to0.time) == fixed_step
+    end
+
+    @testset "Multiple waveforms" begin
+        # create 10 random waveforms with a step of 16ns
+        old_offset = 0.0u"ns"
+        old_step   = 16.0u"ns"
+        t = range(old_offset, step = old_step, length = 100)
+        wfs = ArrayOfRDWaveforms(RDWaveform.(Ref(t), rand.(fill(100,10))))
+        @test length(wfs) == 10
+        @test all(wfs.time .== Ref(t))
+
+        new_offset = 100u"ns"
+        new_step   =   4u"ns"
+        flt = TimeAxisFilter(new_step, new_offset)
+        wfs_new = flt.(wfs)
+        # the broadcasting result and individual result should be identical
+        @test first(wfs_new) == flt(first(wfs))
+        # all time axes should be the same
+        new_t = flt(first(wfs)).time
+        @test all(wfs_new.time .== Ref(new_t))
+        @test first(new_t) == old_offset + new_offset
+        @test step(new_t) == new_step
+    end
+end