From 6953c9bf9e35e96d2168d12f09cb5f5e52bb86b5 Mon Sep 17 00:00:00 2001
From: Joakim Skogholt <onlinereg@skogholt.xyz>
Date: Sat, 20 May 2023 08:55:37 +0200
Subject: [PATCH] Added basic simulation of spectroscopic data

---
 src/MinPakke.jl                  |  3 +++
 src/simulateSpectroscopicData.jl | 43 ++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)
 create mode 100644 src/simulateSpectroscopicData.jl

diff --git a/src/MinPakke.jl b/src/MinPakke.jl
index 6163d27..172b261 100644
--- a/src/MinPakke.jl
+++ b/src/MinPakke.jl
@@ -48,11 +48,14 @@ export TRSegCV
 export TRVirCV
 export TRBidiagDecomp
 
+export simulateSpectrum
+
 
 include("preprocessing.jl")
 include("convenience.jl")
 include("conveniencePlots.jl")
 include("variousRegressionFunctions.jl")
 include("TR.jl")
+include("simulateSpectroscopicData.jl")
 
 end
\ No newline at end of file
diff --git a/src/simulateSpectroscopicData.jl b/src/simulateSpectroscopicData.jl
new file mode 100644
index 0000000..6e0b1cf
--- /dev/null
+++ b/src/simulateSpectroscopicData.jl
@@ -0,0 +1,43 @@
+Lk(waves, c, gamma) = @. gamma / (pi*(waves-c)^2+gamma^2);
+Gk(waves, c, gamma) = @. 1/(sqrt(2*pi)*gamma) * exp(-(waves-c)^2 / (2*gamma^2));
+
+function pseudoVoigtPeak(waves, c, gamma, eta, alpha)
+    
+fk = zeros(length(waves));
+
+for i in 1:length(c)
+        fk += alpha[i] .* pseudoVoigtPeak(waves, c[i], gamma[i], eta[i]);
+end
+    
+return fk
+end
+    
+function pseudoVoigtPeak(waves, c, gamma::Float64, eta::Float64)
+
+fk = eta * Lk(waves, c, gamma) + (1-eta) * Gk(waves, c, gamma);
+
+return fk
+end
+
+function simulateBaseline(waves, a)
+
+vm = zeros(length(waves), length(a));
+
+for i in 0:(length(a)-1)
+    vm[:, i+1] = waves.^i ./ norm(waves.^i);
+end
+
+b = vm * a;
+    
+return b
+end
+
+function simulateSpectrum(waves, c, gamma, eta, alpha, a, sigma=0.0)
+    
+pure_spec = pseudoVoigtPeak(waves, c, gamma, eta, alpha);
+b         = simulateBaseline(waves, a);
+noise     = randn(length(waves)) .* sigma;
+spec      = pure_spec + b + noise;
+return spec, pure_spec, b, noise
+end
+