From e445b1dd3adb0a72e18baac9c8da21c5bcedadb4 Mon Sep 17 00:00:00 2001
From: Carl Philipp Klemm <philipp@uvos.xyz>
Date: Fri, 3 Jun 2022 13:31:10 +0200
Subject: [PATCH] add the ability to add random noise to output spectra

---
 basicmath.cpp            | 16 ++++++++++++++++
 eisgenerator/basicmath.h |  1 +
 main.cpp                 |  9 +++++++--
 options.h                |  5 +++++
 4 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/basicmath.cpp b/basicmath.cpp
index fa68fad..3c52f4c 100644
--- a/basicmath.cpp
+++ b/basicmath.cpp
@@ -1,5 +1,6 @@
 #include "basicmath.h"
 #include <algorithm>
+#include <random>
 
 #include "eistype.h"
 
@@ -115,3 +116,18 @@ std::vector<eis::DataPoint> eis::rescale(const std::vector<eis::DataPoint>& data
 	}
 	return output;
 }
+
+void eis::noise(std::vector<eis::DataPoint>& data, double amplitude, bool relative)
+{
+	std::random_device rd;
+	std::default_random_engine rande(rd());
+
+	for(eis::DataPoint& point : data)
+	{
+		double realNoise = (static_cast<double>(rande() >> 1) / (rande.max() >> 1))*amplitude;
+		point.im.real(relative ? point.im.real()+realNoise/point.im.real() : point.im.real()+realNoise);
+
+		double imgNoise = (static_cast<double>(rande() >> 1) / (rande.max() >> 1))*amplitude;
+		point.im.imag(relative ? point.im.imag()+imgNoise/point.im.imag() : point.im.imag()+imgNoise);
+	}
+}
diff --git a/eisgenerator/basicmath.h b/eisgenerator/basicmath.h
index 0cc0d01..b54a726 100644
--- a/eisgenerator/basicmath.h
+++ b/eisgenerator/basicmath.h
@@ -12,4 +12,5 @@ namespace eis
 	fvalue median(std::vector<fvalue> data);
 	std::complex<fvalue> median(const std::vector<eis::DataPoint>& data);
 	std::vector<eis::DataPoint> rescale(const std::vector<eis::DataPoint>& data, size_t outputSize);
+	void noise(std::vector<eis::DataPoint>& data, double amplitude, bool relative);
 }
diff --git a/main.cpp b/main.cpp
index ffe507a..28fad67 100644
--- a/main.cpp
+++ b/main.cpp
@@ -3,6 +3,7 @@
 #include <chrono>
 #include <cmath>
 
+#include "basicmath.h"
 #include "model.h"
 #include "log.h"
 #include "options.h"
@@ -36,7 +37,8 @@ static void paramSweepCb(std::vector<eis::DataPoint>& data, const std::vector<fv
 		std::cout<<'.'<<std::flush;
 }
 
-static void runSweep(const std::string& modelString, eis::Range omega, bool normalize = false, bool reduce = false, bool hertz = false, bool invert = false)
+static void runSweep(const std::string& modelString, eis::Range omega, bool normalize = false,
+					 bool reduce = false, bool hertz = false, bool invert = false, double noise = 0)
 {
 	std::vector<eis::DataPoint> results;
 
@@ -64,6 +66,9 @@ static void runSweep(const std::string& modelString, eis::Range omega, bool norm
 		eis::Log(eis::Log::INFO)<<"results:";
 	}
 
+	if(noise > 0)
+		eis::noise(results, noise, false);
+
 	std::cout<<(hertz ? "freqency" : "omega")<<",real,im\n";
 
 	for(const eis::DataPoint& res : results)
@@ -107,7 +112,7 @@ int main(int argc, char** argv)
 	switch(config.mode)
 	{
 		case MODE_SWEEP:
-			runSweep(config.modelStr, config.omegaRange, config.normalize, config.reduce, config.hertz, config.invert);
+			runSweep(config.modelStr, config.omegaRange, config.normalize, config.reduce, config.hertz, config.invert, config.noise);
 			break;
 		case MODE_PARAM_SWEEP:
 			runParamSweep();
diff --git a/options.h b/options.h
index ac04ebb..b5f5cc2 100644
--- a/options.h
+++ b/options.h
@@ -24,6 +24,7 @@ static struct argp_option options[] =
   {"reduce",    'r', 0,      0,  "reduce values to \"interesting\" range" },
   {"hz",        'h', 0,      0,  "freqency values as temporal frequency instead of angular frequency"},
   {"invert",        'i', 0,      0,  "inverts the imaginary axis"},
+  {"noise",        'x', "[AMPLITUDE]",      0,  "add noise to output"},
   { 0 }
 };
 
@@ -42,6 +43,7 @@ struct Config
 	bool reduce = false;
 	bool hertz = false;
 	bool invert = false;
+	double noise = 0;
 
 	Config(): omegaRange(1, 1001, 1, true)
 	{}
@@ -123,6 +125,9 @@ parse_opt (int key, char *arg, struct argp_state *state)
 		}
 		break;
 	}
+	case 'x':
+		config->noise = std::stod(std::string(arg));
+		break;
 	default:
 		return ARGP_ERR_UNKNOWN;
 	}
-- 
GitLab