I have time series data in which the recorded values capture a signal (i.e., peak) that emerges from an elevated noise floor (i.e., baseline; can be modeled as a straight line). The shape of the peak in the various time series is identical. What distinguishes the time series are the parameters of the baselines (i.e., slope and intercept), the scale of each peak, and the detection times. However, in this specific use case, both the difference in detection times between the time series (δt) and the peaks scaling factors are a priori known. The example code below simulates two such time series (Figure column 1).
I would like to deduce the slope and the intercept of each time series’ baseline (Figure column 2) and approximate the shape of the peak with a spline (Figure column 3) using the data from all time series (in the example below, two) simultaneously.
I was able to approximate the signal in noisy data with a spline in the absence of a baseline. I was also able to fit baselines to flat noise floors (i.e., without a peak) while imposing certain constraints (e.g., the inferred baseline must not predict negative values of the noise floor within the time range). However, I am struggling to piece these components together and solve this optimization problem with the available packages in Julia. Any assistance in this regard would be greatly appreciated.
As indicated above, the optimization problem should ideally impose certain constraints, namely:
- Signals predicted by the deduced baselines cannot be negative within the analyzed time range.
- Signals of the predicted peaks cannot be negative in any of the time series.
- The weight of each sample point should be adjusted based on its associated error (optional).
using CairoMakie, Distributions, Random
################
baseline(intercept, slope, δt) = t -> intercept .+ slope * (t + δt)
signal(baseline, peak, δt, scale) = t -> baseline(t) + scale * peak(t + δt)
noisy!(signal, noise_factor) = t -> signal(t) + rand.(Normal.(0, σ_error(signal(t), noise_factor)))
σ_error(signal, noise_factor) = @. noise_factor * sqrt(signal) / sqrt(2/π)
################
ts = 5:5:100
noise_factor = 1.7
# Some signal distribution
peak(t) = 6 * 10^4 * pdf.(Normal(50, 5), t) # peak sigmal (to be inferred from the empirical data)
# Time series data 1
δt₁ = -0.5 # ion detection time shift relativ to scantime (known)
bl₁ = baseline(2000, -5.1, δt₁) # baseline signal (to be inferred from the empirical data)
scale₁ = 1 # peak signal scaling factor (known)
# Time series data 2
δt₂ = -4.2 # ion detection time shift relativ to scantime (known)
bl₂ = baseline(1000, 2, δt₂) # baseline signal (to be inferred from the empirical data)
scale₂ = 0.2 # peak signal scaling factor (known)
# Generate the signals of each ion by combining baseline signal and peak signal
signal₁ = signal(bl₁, peak, δt₁, scale₁) # ion 1
signal₂ = signal(bl₂, peak, δt₂, scale₂) # ion 2
fig = Figure(resolution = (800, 500))
# Ion 1
ax₁ = Axis(fig[1, 1])
ylims!(ax₁, 500, 7000)
lines!(ax₁, 1..100, signal₁, color=:red) # true signal
scatter!(ax₁, ts, noisy!(signal₁, noise_factor), color=:blue) # detected signal
# Ion 2
ax₂ = Axis(fig[2, 1])
ylims!(ax₂, 500, 2200)
lines!(ax₂, 1..100, signal₂, color=:red) # true signal
scatter!(ax₂, ts, noisy!(signal₂, noise_factor), color=:blue) # detected signal
# To be inferred baseline in time series data of ion 1
ax₃ = Axis(fig[1, 2])
ylims!(ax₃, 500, 7000)
lines!(ax₃, 5..100, bl₁)
# To be inferred baseline in time series data of ion 2
ax₄ = Axis(fig[2, 2])
ylims!(ax₄, 500, 2200)
lines!(ax₄, 5..100, bl₂)
# To be inferred empirical peak shape from signals in time series data of ion 1 and ion 2
ax₅ = Axis(fig[1:2, 3])
lines!(ax₅, 5..100, peak)
display(fig)