Added baseline correction

TODO: Test

Project.toml

@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

docs/src/reference.md

@@ -75,6 +75,12 @@ adjust_phase
 ernst
 ```
 
+### Baseline correction
+
+```@docs
+baseline_als
+```
+
 ## Plotting
 
 ```@docs

src/MagneticResonanceSignals.jl

@@ -10,6 +10,7 @@ using FFTW
 
 using LinearAlgebra
 using Optim
+using SparseArrays
 
 # Plotting
 using RecipesBase
@@ -91,6 +92,7 @@ export
     simple_averaging,
     # Fourier transform
     spectrum,
+    baseline_als
 
 # Plotting
 export

src/processing.jl

@@ -145,4 +145,34 @@ function single_spectrum_version(spectrum::AbstractArray)
     result.minimizer
 end
 
+"""
+    function baseline_als(spectrum, lambda::Float64, p::Float64; niter::Int64=10)
+
+Baseline correction based on ALS (Asymmetric Least Square).
+
+Lambda is 2nd derivative constraint, which contributes to the smoothness, while p is
+weighting of positive residuals, which affect to the asymmetry.
+
+This implementation is inspired by https://stackoverflow.com/a/29185844/5023889
+"""
+function baseline_als(spectrum::AbstractArray, lambda::Float64, p::Float64; niter::Int64=10)
+    L = size(spectrum)[1]
+    x = diff(Matrix{Float64}(I, L, L), dims=2)
+    D = sparse(diff(x, dims=2))
+    w = ones(L)
+    z = nothing
+    for i in 1:niter
+        W = sparse(diagm(ones(L)))
+        Z = W + lambda * (D * transpose(D))
+        z = Z \ (w .* spectrum)
+        w = (p .* (spectrum > z)) + ((1-p) .* (spectrum < z))
+    end
+    z
+end
 
+function baseline_als(spectrum::AxisArray, lambda::Float64, p::Float64; niter::Int64=10)
+    AxisArray(
+        baseline_als(spectrum, lambda, p; niter=niter),
+        AxisArrays.axes(y)
+    )
+end