diff --git a/src/LegendDSP.jl b/src/LegendDSP.jl index a359e1d..26f6a56 100644 --- a/src/LegendDSP.jl +++ b/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") diff --git a/src/moving_window_multi.jl b/src/moving_window_multi.jl index cf80bf3..2647a08 100644 --- a/src/moving_window_multi.jl +++ b/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...)``` diff --git a/src/ms_filter.jl b/src/ms_filter.jl index 5af1ea8..9ad8536 100644 --- a/src/ms_filter.jl +++ b/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))) diff --git a/src/sgw_filter.jl b/src/sgw_filter.jl index ccfd4a6..403da4e 100644 --- a/src/sgw_filter.jl +++ b/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]) \ No newline at end of file +@inline _w(x::Real, v::AbstractVector) = 1 - v[1]/(1 + v[2]*x^v[3]) \ No newline at end of file diff --git a/src/wh_filter.jl b/src/wh_filter.jl index 22696a9..2313f19 100644 --- a/src/wh_filter.jl +++ b/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