Implement suggestions from #39

src/LegendDSP.jl

@@ -40,6 +40,7 @@ include("saturation.jl")
 include("interpolation.jl")
 include("intersect_maximum.jl")
 include("dsp_sipm.jl")
+include("dsp_sipm_optimization.jl")
 include("dsp_routines.jl")
 include("dsp_puls.jl")
 include("haar_filter.jl")

src/moving_window_multi.jl

@@ -3,10 +3,21 @@ import RadiationDetectorDSP: AbstractConvFilterInstance,
     AbstractRadSigFilterInstance, LinearFiltering, _filterlen, _floattype
 
 """
-    struct MovingWindowMulti <: AbstractRadFIRFilter
+    struct MovingWindowMultiFilter <: AbstractRadFIRFilter
 
 apply left and right moving average windows to signal.
-
+Working example:
+```julia
+using RadiationDetectorSignals
+using Unitful
+
+signal = rand(100)
+t = range(0u"ms", 20u"ms", 100)
+
+wf = RDWaveform(t, signal)
+flt = MovingWindowMultiFilter(1u"ms")
+wf_new = flt(wf)
+```
 Constructors:
 
 * ```$(FUNCTIONNAME)(; fields...)```

src/ms_filter.jl

@@ -2,7 +2,7 @@
 """
     struct ModifiedSincFilter{T<:RealQuantity} <: AbstractRadFIRFilter{}
 
-A [Modified-Sinc filter](https://doi.org/10.1021/acsmeasuresciau.3c00017)
+A [Modified-Sinc filter](https://doi.org/10.1021/acsmeasuresciau.1c00054)
 Working example
 ```julia
 using RadiationDetectorSignals
@@ -35,7 +35,7 @@ Base.@kwdef struct ModifiedSincFilter{T<:RealQuantity} <: AbstractRadFIRFilter
     "half-width of the kernel"
     m::T = 3
     function ModifiedSincFilter(d::Int, m::T) where T
-        s1 = "degree of ModifiedSincFilter needs to be even and be in the interval [2, 10]"
+        s1 = "degree of ModifiedSincFilter needs to be even and be in the interval [2, $(D_MAX)]"
         @assert valid_degree(d) s1
 
         new{T}(d, m)
@@ -59,7 +59,7 @@ Adapt.adapt_structure(to, flt::ModifiedSincFilter) = flt
 function fltinstance(flt::ModifiedSincFilter, fi::SamplingInfo{T}) where {T}
     _m = round(Int, ustrip(NoUnits, flt.m/step(fi.axis)))
     m_min = flt.d/2 + 2
-    s2 = "size of kernel to small for given degree"
+    s2 = "size of kernel too small for given degree"
     @assert _m >= m_min s2
 
     ModifiedSincFilterInstance{T}(flt.d, _m, length(fi.axis))
@@ -130,12 +130,12 @@ const MS_COEFFS = Dict(
 end
 
 """
-    makeFitWeights(d::Int, m::Int)
+    makeFitWeights(d::Int, m::Int) where {U}
 
 return the weights for the linear fit user for linear extrapolation at 
 at the right boundary.
 """
-function makeFitWeights(d::U, m::Int) where U
+function makeFitWeights(d::U, m::Int) where {U <: Integer}
     firstZero = (m+1)/(1.5 + d/2)
     β = 0.7 + 0.14*exp(-0.6*(d-4))
     l = ceil(Int, firstZero*β)
@@ -174,7 +174,7 @@ end
     weighted_linear_reg(w::AbstractVector, y::AbstractVector)
 
 Do a weighted linear regression of the data `y` with weights `w` and return
-the offset and slope. its assumed that the 
+the offset and slope. its assumed that x is a range from 1 to `length(y)`
 """
 function weighted_linear_reg(w::AbstractVector, y::AbstractVector)
     T = RadiationDetectorDSP._floattype(promote_type(eltype(w), eltype(y)))

src/sgw_filter.jl

@@ -1,7 +1,7 @@
 """
     struct WeightedSavitzkyGolayFilter{T<:RealQuantity} <: AbstractRadFIRFilter
 
-A [Weighted-Savitzky-Golay filter](https://doi.org/10.1021/acsmeasuresciau.3c00017).
+A [Weighted-Savitzky-Golay filter](https://doi.org/10.1021/acsmeasuresciau.1c00054).
 Working example:
 ```julia
 using RadiationDetectorSignals
@@ -65,7 +65,7 @@ const SGW_COEFFS = Dict(
     4 => [0.62303, 0.25310, -7.07317]    # HANNCUBE
 )
 
-const WEIGHTTYPES = (0, 1, 2, 3, 4)
+const WEIGHTTYPES = keys(SGW_COEFFS)
 
 function fltinstance(flt::WeightedSavitzkyGolayFilter, fi::SamplingInfo{T}
     ) where {T}
@@ -79,7 +79,7 @@ function fltinstance(flt::WeightedSavitzkyGolayFilter, fi::SamplingInfo{T}
     @assert length(fi.axis) >= 2*m + 1 msg1
 
     d_max = 2*m
-    msg2 = "degree to big for given kernel size; max degree: $(d_max)"
+    msg2 = "degree too big for given kernel size; max degree: $(d_max)"
     @assert flt.degree <= 2*m msg2
 
     WeightedSavitzkyGolayFilterInstance{T}(m, flt.degree, flt.weightType, l)
@@ -127,7 +127,7 @@ Adapt.adapt_structure(to, flt::WeightedSavitzkyGolayFilter) = flt
             sum = fma(kernel[j], x[idx], sum)
         end
         y[i] = sum
-        sum = 0.
+        sum = zero(T)
     end
     # near boundary points at the right
     sum = zero(T)
@@ -145,7 +145,7 @@ Adapt.adapt_structure(to, flt::WeightedSavitzkyGolayFilter) = flt
             sum = fma(kernel[j], x[idx], sum)
         end
         y[L-m+i] = sum
-        sum = 0.
+        sum = zero(T)
     end
     y
 end
@@ -205,27 +205,27 @@ function _unsafe_dot_w(x::AbstractVector, y::AbstractVector, w::AbstractVector)
     dotprod
 end
 
-function weightFunctionScale(missingFrac, weightType::Int)
+function weightFunctionScale(missingFrac::Real, weightType::Int)
     coeffs = SGW_COEFFS[weightType]
     missingFrac <= 0 ? 1 : _w(missingFrac, coeffs)
 end
 
-function weight_function(weight_type::Int, x::T) where T
-    decay = 2.0 # for GAUSS only
+function weight_function(weight_type::Int, x::T) where {T <: Real}
     if x <= -0.999999999999 || x >= 0.999999999999
-        return 0.0
+        return zero(T)
     elseif weight_type == 0
-        return 1.0
+        return one(T)
     elseif weight_type == 1
+        decay = T(2.0)
         return exp(-x^2 * decay) + exp(-(x - 2)^2 * decay) + exp(-(x + 2)^2 * decay) -
                2 * exp(-decay) - exp(-9 * decay) # Gaussian-like alpha=2
     elseif weight_type == 2
-        return cos(0.5 * π * x)^2 # Hann
+        return cos(π/2 * x)^2 # Hann
     elseif weight_type == 3
-        return (cos(0.5 * π * x))^4 # Hann-squared
+        return (cos(π/2 * x))^4 # Hann-squared
     else weight_type == 4
-        return (cos(0.5 * π * x))^6 # Hann-cube
+        return (cos(π/2 * x))^6 # Hann-cube
     end
 end
 
-@inline _w(x, v::AbstractVector) = 1 - v[1]/(1 + v[2]*x^v[3])
+@inline _w(x::Real, v::AbstractVector) = 1 - v[1]/(1 + v[2]*x^v[3])

src/wh_filter.jl

@@ -1,8 +1,8 @@
 """
 
-    WhittakerHendersonFilter <: AbstractRadSigFilter{LinearFiltering}
+    struct WhittakerHendersonFilter <: AbstractRadSigFilter{LinearFiltering}
 
-A [Whittaker-Henderson filter](https://doi.org/10.1021/acsmeasuresciau.3c00017).
+A [Whittaker-Henderson filter](https://doi.org/10.1021/acsmeasuresciau.1c00054).
 Working example:
 ```julia
 using RadiationDetectorSignals
@@ -45,7 +45,7 @@ end
 function rdfilt!(output::AbstractVector, fi::WhittakerHendersonFilterInstance, input::AbstractVector)
     A = makeIPlusLambdaDprimeD(fi.λ, fi.p, length(output))
     choleskyL!(A)
-    solve!(output, A, input)
+    banded_cholesky_solve!(output, A, input)
 end
 
 RadiationDetectorDSP.flt_output_smpltype(fi::WhittakerHendersonFilterInstance{T}) where {T} = RadiationDetectorDSP._floattype(flt_input_smpltype(fi))
@@ -115,13 +115,13 @@ function choleskyL!(b::AbstractMatrix{T}) where {T}
 end
 
 """
-    solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
+    banded_cholesky_solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
 
-solve the lienar equation AA'x = y and store the result in x, where
+solve the linear equation `AA'x = y` and store the result in `x`, where
 `A` is the cholesky decomposition of a banded centro symmetric matrix.
-A[1, i] is the diagonal, A[2, i] is the first subdiagonal and so on.
+`A[1, i]` is the diagonal, `A[2, i]` is the first subdiagonal and so on.
 """
-function solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
+function banded_cholesky_solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     n = size(A, 2)
     dmax = size(A, 1) - 1
     T = RadiationDetectorDSP._floattype(eltype(x))
@@ -137,7 +137,7 @@ function solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
 
     # Backward substitution
     @inbounds for i in reverse(axes(A, 2))
-        sum = 0.0
+        sum = zero(T)
         @simd for j in (i+1):min(i+dmax, n)
             sum = fma(A[j-i+1, i], x[j], sum)
         end