From 81000d6542783c320a58dae8a7878086d5677c59 Mon Sep 17 00:00:00 2001 From: adil Date: Wed, 23 Apr 2025 10:46:48 -0400 Subject: [PATCH 01/12] Initial commmit for implementing working aWGN in Katydid, need to set a working seed for random noise generation --- CMakeLists.txt | 4 +- Examples/ConfigFiles/meg_noise_config.yaml | 57 +++++++++++++++++++ Source/Simulation/CMakeLists.txt | 20 +++---- Source/Simulation/KTDCOffsetGenerator.cc | 16 +++--- Source/Simulation/KTGaussianNoiseGenerator.cc | 11 ++-- Source/Simulation/KTTSGenerator.cc | 9 ++- Source/Simulation/KTTSGenerator.hh | 6 +- Source/Time/KTEggProcessor.hh | 2 +- 8 files changed, 94 insertions(+), 31 deletions(-) create mode 100644 Examples/ConfigFiles/meg_noise_config.yaml diff --git a/CMakeLists.txt b/CMakeLists.txt index 33dfd2146..f619526ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -285,7 +285,7 @@ include_directories (BEFORE ${PROJECT_SOURCE_DIR}/Source/IO/ROOTTreeWriter ${PROJECT_SOURCE_DIR}/Source/IO/TerminalWriter ${PROJECT_SOURCE_DIR}/Source/Time - #${PROJECT_SOURCE_DIR}/Source/Simulation + ${PROJECT_SOURCE_DIR}/Source/Simulation ${PROJECT_SOURCE_DIR}/Source/Evaluation ${PROJECT_SOURCE_DIR}/Source/Transform ${PROJECT_SOURCE_DIR}/Source/SpectrumAnalysis @@ -297,7 +297,7 @@ add_subdirectory (Source/Utility) add_subdirectory (Source/Data) add_subdirectory (Source/IO) add_subdirectory (Source/Time) -#add_subdirectory (Source/Simulation) +add_subdirectory (Source/Simulation) #add_subdirectory (Source/Evaluation) add_subdirectory (Source/Transform) add_subdirectory (Source/SpectrumAnalysis) diff --git a/Examples/ConfigFiles/meg_noise_config.yaml b/Examples/ConfigFiles/meg_noise_config.yaml new file mode 100644 index 000000000..f8140805b --- /dev/null +++ b/Examples/ConfigFiles/meg_noise_config.yaml @@ -0,0 +1,57 @@ +processor-toolbox: + processors: + - type: egg-processor + name: egg + - type: forward-fftw + name: fft + - type: root-spectrogram-writer + name: writer + - type: basic-root-writer + name: one_d_hist + - type: realistic-noise-generator + name: noise_gen + + + connections: + - signal: "egg:header" + slot: "fft:header" + - signal: "egg:ts" + slot: "noise_gen:slice" + - signal: "noise_gen:slice" + slot: "fft:ts-fftw" + - signal: "fft:fft" + slot: "one_d_hist:fs-fftw" + - signal: "fft:fft" + slot: "one_d_hist:fs-fftw-power-dist" + - signal: "egg:egg-done" + slot: "writer:write-file" + + run-queue: + - egg + +egg: + filename: "/work/locust_mc_Seed640_Angle90.0000000_Pos0.0040000_Energy18600.000000.egg" + egg-reader: egg3 + slice-size: 4096 + number-of-slices: 1 #0 will run thru everything + +fft: + transform-flag: ESTIMATE + +noise_gen: + mean: 0.0 #to get a meaningful difference as of oct 1 2024 these values need to be pretty large + sigma: 1.0 + + +writer: + output-file: "/work/ps_output_noise.root" + file-flag: recreate + min-time: 0.0 + max-time: 0.1 + min-freq: 0.0 + max-freq: 100.0e6 + +one_d_hist: + output-file: "/work/test_output_hist.root" + file-flag: recreate + diff --git a/Source/Simulation/CMakeLists.txt b/Source/Simulation/CMakeLists.txt index be8f37b0a..bcbca15e5 100644 --- a/Source/Simulation/CMakeLists.txt +++ b/Source/Simulation/CMakeLists.txt @@ -15,16 +15,16 @@ set (SIMULATION_SOURCEFILES ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.cc ) -#if (ROOT_FOUND) -# set (SIMULATION_HEADERFILES -# ${SIMULATION_HEADERFILES} -# ${CMAKE_CURRENT_SOURCE_DIR}/KT.hh -# ) -# set (SIMULATION_SOURCEFILES -# ${SIMULATION_SOURCEFILES} -# ${CMAKE_CURRENT_SOURCE_DIR}/KT.cc -# ) -#endif (ROOT_FOUND) +if (ROOT_FOUND) + set (SIMULATION_HEADERFILES + ${SIMULATION_HEADERFILES} + ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.hh + ) + set (SIMULATION_SOURCEFILES + ${SIMULATION_SOURCEFILES} + ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.cc + ) +endif (ROOT_FOUND) set (KATYDID_LIBS KatydidUtility diff --git a/Source/Simulation/KTDCOffsetGenerator.cc b/Source/Simulation/KTDCOffsetGenerator.cc index d648451aa..90a32f6dd 100644 --- a/Source/Simulation/KTDCOffsetGenerator.cc +++ b/Source/Simulation/KTDCOffsetGenerator.cc @@ -31,22 +31,22 @@ namespace Katydid { } - bool KTDCOffsetGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) + bool KTDCOffsetGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) //changed scarab::param_node to scarab::param { if (node == NULL) return false; - const KTParamArray* offsetPairs = node->ArrayAt("offsets"); + const scarab::param_array* offsetPairs = node-> array_at("offsets"); //changing this based on DevGuide &KTCorrelator if (offsetPairs != NULL) { - for (KTParamArray::const_iterator pairIt = offsetPairs->Begin(); pairIt != offsetPairs->End(); ++pairIt) + for (scarab::param_array::const_iterator pairIt = offsetPairs->begin(); pairIt != offsetPairs->end(); ++pairIt) { - if (! ((*pairIt)->IsArray() && (*pairIt)->AsArray().Size() == 2)) + if (! ((*pairIt)->is_array() && (*pairIt)->as_array().size() == 2)) { - KTERROR(genlog, "Invalid pair: " << (*pairIt)->ToString()); + KTERROR(genlog, "Invalid pair: " << (*pairIt)->to_string()); return false; } - UIntDoublePair pair((*pairIt)->AsArray().GetValue< unsigned >(0), (*pairIt)->AsArray().GetValue< double >(1)); - if (fOffsets.size() <= pair.first) fOffsets.resize(pair.first + 1); + UIntDoublePair pair((*pairIt)->as_array().get_value< unsigned >(0), (*pairIt)->as_array().get_value< double >(1)); + if (fOffsets.size() <= pair.first) fOffsets.resize(pair.first + 1); //keeping fOffsets for now bc it is defined in the .hh fOffsets[pair.first] = pair.second; } } @@ -73,7 +73,7 @@ namespace Katydid for (unsigned iBin = 0; iBin < sliceSize; ++iBin) { - timeSeries->SetValue(iBin, fOffsets[iComponent] + timeSeries->GetValue(iBin)); + timeSeries->SetValue(iBin, fOffsets[iComponent] + timeSeries->GetValue(iBin)); //is this GetValue diff } } diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 541c96892..4067afb9f 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -32,7 +32,7 @@ namespace Katydid { } - bool KTGaussianNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) + bool KTGaussianNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) //I changed this to scarab::param somewhere else { if (node == NULL) return false; @@ -40,13 +40,14 @@ namespace Katydid input_type mean = node->get_value< input_type >("mean", fRNG.mean()); input_type sigma = node->get_value< input_type >("sigma", fRNG.sigma()); fRNG.param(KTRNGGaussian<>::param_type(mean, sigma)); + // fRNG.SetSeed(1) we still need to add a set seed function here but this one does not work return true; } bool KTGaussianNoiseGenerator::GenerateTS(KTTimeSeriesData& data) { - //const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); + const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); const unsigned sliceSize = data.GetTimeSeries(0)->GetNTimeBins(); unsigned nComponents = data.GetNComponents(); @@ -61,12 +62,12 @@ namespace Katydid continue; } - //double binCenter = 0.5 * binWidth; + double binCenter = 0.5 * binWidth; for (unsigned iBin = 0; iBin < sliceSize; iBin++) { timeSeries->SetValue(iBin, fRNG() + timeSeries->GetValue(iBin)); - //binCenter += binWidth; - //KTDEBUG(genlog, iBin << " " << timeSeries->GetValue(iBin)); + binCenter += binWidth; + KTDEBUG(genlog, iBin << " " << timeSeries->GetValue(iBin)); } } diff --git a/Source/Simulation/KTTSGenerator.cc b/Source/Simulation/KTTSGenerator.cc index ae0f92bef..224331ea6 100644 --- a/Source/Simulation/KTTSGenerator.cc +++ b/Source/Simulation/KTTSGenerator.cc @@ -12,12 +12,14 @@ #include "KTLogger.hh" #include "KTProcSummary.hh" #include "param.hh" +#include "time.hh" #include "KTSliceHeader.hh" #include "KTTimeSeriesData.hh" #include "KTTimeSeriesFFTW.hh" #include "KTTimeSeriesReal.hh" -#include "thorax.hh" +//#include "thorax.hh" it appears this isn't defined anywhere? + #include @@ -165,7 +167,7 @@ namespace Katydid } - newHeader->SetTimestamp(get_absolute_time_string()); + newHeader->SetTimestamp(scarab::get_absolute_time_string()); //this is a membervariableref... whatever that means return newHeader; } @@ -194,7 +196,8 @@ namespace Katydid for (unsigned iComponent = 0; iComponent < fNChannels; ++iComponent) { - sliceHeader.SetTimeStamp((uint64_t)(sliceHeader.GetTimeInRun() * (double)NSEC_PER_SEC), iComponent); // TODO: change this to 1e3 when switch to usec is made + sliceHeader.SetTimeStamp((uint64_t)(sliceHeader.GetTimeInRun() * (double)CLOCKS_PER_SEC), iComponent); // TODO: change this to 1e3 when switch to usec is made + // changed above to clocks_per_sec bc that was what the builder suggested... how much do we trust it sliceHeader.SetAcquisitionID(0); sliceHeader.SetRecordID(0); } diff --git a/Source/Simulation/KTTSGenerator.hh b/Source/Simulation/KTTSGenerator.hh index ed2678d36..babca06c3 100644 --- a/Source/Simulation/KTTSGenerator.hh +++ b/Source/Simulation/KTTSGenerator.hh @@ -123,9 +123,11 @@ namespace Katydid // Signals //*************** private: - Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; + //Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; this was the old one + Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; //this is the new one according to KTEggProcessor.hh in Source/Time could still be wrong tho + //compare in soruce/time directory in KTeggprocessor make the header signal look the same also check done signal in line129 Nymph::KTSignalData fDataSignal; - Nymph::KTSignalOneArg< void > fDoneSignal; + Nymph::KTSignalOneArg< void > fDoneSignal; //this looks the same as last Nymph::KTSignalOneArg< const KTProcSummary* > fSummarySignal; }; diff --git a/Source/Time/KTEggProcessor.hh b/Source/Time/KTEggProcessor.hh index 2af98b7a0..17a11eb10 100644 --- a/Source/Time/KTEggProcessor.hh +++ b/Source/Time/KTEggProcessor.hh @@ -1,5 +1,5 @@ /** - @file KTEggProcessor.hh + @file KTEggprocessor.hh @brief Contains KTEggProcessor @details Iterates over slices in an Egg file. @author: N. S. Oblath From 22c4938d23794a9e13dfea834478b056ab20f4f6 Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 17:29:50 -0400 Subject: [PATCH 02/12] Added a fix for handling noise addition properly (resolves mirrored spectrum issue in PS) --- Source/Simulation/KTGaussianNoiseGenerator.cc | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 4067afb9f..0766b9666 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -11,6 +11,7 @@ #include "KTMath.hh" #include "KTTimeSeriesData.hh" #include "KTTimeSeries.hh" +#include "KTTimeSeriesFFTW.hh" #include @@ -63,13 +64,21 @@ namespace Katydid } double binCenter = 0.5 * binWidth; - for (unsigned iBin = 0; iBin < sliceSize; iBin++) + if (auto* tsFFTW = dynamic_cast(timeSeries)) // Handling complex FFTW time series correctly { - timeSeries->SetValue(iBin, fRNG() + timeSeries->GetValue(iBin)); - binCenter += binWidth; - KTDEBUG(genlog, iBin << " " << timeSeries->GetValue(iBin)); + for (unsigned iBin = 0; iBin < sliceSize; ++iBin) + { + tsFFTW->SetRect(iBin, tsFFTW->GetReal(iBin) + fRNG(), tsFFTW->GetImag(iBin) + fRNG()); // Complex white-Gaussian noise + } + } + else + { + for (unsigned iBin = 0; iBin < sliceSize; ++iBin) + { + timeSeries->SetValue(iBin, timeSeries->GetValue(iBin) + fRNG()); + } + } } - } return true; } From a12bb5585d9ecbaba6c50572c72da4dcf4d3a194 Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 18:07:06 -0400 Subject: [PATCH 03/12] Added a way to incorporate seed for the RNG() --- Source/Simulation/KTGaussianNoiseGenerator.cc | 4 ++-- Source/Utility/KTRandom.hh | 8 ++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 0766b9666..d461b753e 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -41,8 +41,8 @@ namespace Katydid input_type mean = node->get_value< input_type >("mean", fRNG.mean()); input_type sigma = node->get_value< input_type >("sigma", fRNG.sigma()); fRNG.param(KTRNGGaussian<>::param_type(mean, sigma)); - // fRNG.SetSeed(1) we still need to add a set seed function here but this one does not work - + fRNG.SetSeed(node->get_value("seed")); + return true; } diff --git a/Source/Utility/KTRandom.hh b/Source/Utility/KTRandom.hh index c82e2780d..ca6a4b654 100644 --- a/Source/Utility/KTRandom.hh +++ b/Source/Utility/KTRandom.hh @@ -109,6 +109,7 @@ namespace Katydid Engine* GetEngine() const; void SetEngine(Engine* rng); + void SetSeed(unsigned seed); protected: Engine* fEngine; @@ -132,6 +133,13 @@ namespace Katydid return; } + template< class Engine > + inline void KTRNGDistribution< Engine >::SetSeed(unsigned seed) + { + fEngine->SetSeed(seed); + return; + } + template< class Engine > inline bool KTRNGDistribution< Engine >::Configure(const scarab::param_node* node) { From de31c1e914464eb85c7c30084a8934f395a2aa0a Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 18:20:11 -0400 Subject: [PATCH 04/12] Added an example config file for applying aWGN --- .../ConfigFiles/KatydidPSConfig_noise.yaml | 47 +++++++++++++++ Examples/ConfigFiles/meg_noise_config.yaml | 57 ------------------- 2 files changed, 47 insertions(+), 57 deletions(-) create mode 100644 Examples/ConfigFiles/KatydidPSConfig_noise.yaml delete mode 100644 Examples/ConfigFiles/meg_noise_config.yaml diff --git a/Examples/ConfigFiles/KatydidPSConfig_noise.yaml b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml new file mode 100644 index 000000000..3a9f422ac --- /dev/null +++ b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml @@ -0,0 +1,47 @@ +processor-toolbox: + processors: + - type: egg-processor + name: egg + - type: forward-fftw + name: fft + - type: gaussian-noise-generator + name: noise_gen + - type: convert-to-power + name: to-ps + - type: basic-root-writer + name: writer + + connections: + - signal: "egg:header" + slot: "fft:header" + - signal: "egg:ts" + slot: "noise_gen:slice" + - signal: "noise_gen:slice" + slot: "fft:ts-fftw" + - signal: "fft:fft" + slot: "to-ps:fs-fftw-to-psd" + - signal: "to-ps:psd" + slot: "writer:ps" + - signal: "egg:ts" + slot: "writer:ts" + + run-queue: + - egg + +egg: + filename: "locust_mc_Seed600_Angle90.000000000_Pos0.0040000_Energy18563.251000000.egg" + egg-reader: egg3 + slice-size: 4096 + number-of-slices: 1 + +fft: + transform-flag: ESTIMATE + +noise_gen: + mean: 0.0 + sigma: 0.02 + seed: 12345 + +writer: + output-file: "ps_output_noise.root" + file-flag: recreate diff --git a/Examples/ConfigFiles/meg_noise_config.yaml b/Examples/ConfigFiles/meg_noise_config.yaml deleted file mode 100644 index f8140805b..000000000 --- a/Examples/ConfigFiles/meg_noise_config.yaml +++ /dev/null @@ -1,57 +0,0 @@ -processor-toolbox: - processors: - - type: egg-processor - name: egg - - type: forward-fftw - name: fft - - type: root-spectrogram-writer - name: writer - - type: basic-root-writer - name: one_d_hist - - type: realistic-noise-generator - name: noise_gen - - - connections: - - signal: "egg:header" - slot: "fft:header" - - signal: "egg:ts" - slot: "noise_gen:slice" - - signal: "noise_gen:slice" - slot: "fft:ts-fftw" - - signal: "fft:fft" - slot: "one_d_hist:fs-fftw" - - signal: "fft:fft" - slot: "one_d_hist:fs-fftw-power-dist" - - signal: "egg:egg-done" - slot: "writer:write-file" - - run-queue: - - egg - -egg: - filename: "/work/locust_mc_Seed640_Angle90.0000000_Pos0.0040000_Energy18600.000000.egg" - egg-reader: egg3 - slice-size: 4096 - number-of-slices: 1 #0 will run thru everything - -fft: - transform-flag: ESTIMATE - -noise_gen: - mean: 0.0 #to get a meaningful difference as of oct 1 2024 these values need to be pretty large - sigma: 1.0 - - -writer: - output-file: "/work/ps_output_noise.root" - file-flag: recreate - min-time: 0.0 - max-time: 0.1 - min-freq: 0.0 - max-freq: 100.0e6 - -one_d_hist: - output-file: "/work/test_output_hist.root" - file-flag: recreate - From 8b0b99d561e7a85e8c68ca03a080b09dd9021b42 Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 18:28:11 -0400 Subject: [PATCH 05/12] Minor fix in KTRandom.hh --- Source/Utility/KTRandom.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utility/KTRandom.hh b/Source/Utility/KTRandom.hh index ca6a4b654..a3e7689dc 100644 --- a/Source/Utility/KTRandom.hh +++ b/Source/Utility/KTRandom.hh @@ -109,7 +109,7 @@ namespace Katydid Engine* GetEngine() const; void SetEngine(Engine* rng); - void SetSeed(unsigned seed); + void SetSeed(unsigned seed); protected: Engine* fEngine; From 8265bdd48566d12da2b3828774f3aee4d7b4c93e Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 18:48:45 -0400 Subject: [PATCH 06/12] Removed commits and made minor changes before merging to /develop --- Source/Simulation/KTDCOffsetGenerator.cc | 8 ++++---- Source/Simulation/KTGaussianNoiseGenerator.cc | 6 +++--- Source/Simulation/KTTSGenerator.cc | 10 ++++------ Source/Simulation/KTTSGenerator.hh | 6 ++---- Source/Time/KTEggProcessor.hh | 2 +- 5 files changed, 14 insertions(+), 18 deletions(-) diff --git a/Source/Simulation/KTDCOffsetGenerator.cc b/Source/Simulation/KTDCOffsetGenerator.cc index 90a32f6dd..b8d0f85fd 100644 --- a/Source/Simulation/KTDCOffsetGenerator.cc +++ b/Source/Simulation/KTDCOffsetGenerator.cc @@ -31,11 +31,11 @@ namespace Katydid { } - bool KTDCOffsetGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) //changed scarab::param_node to scarab::param + bool KTDCOffsetGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) { if (node == NULL) return false; - const scarab::param_array* offsetPairs = node-> array_at("offsets"); //changing this based on DevGuide &KTCorrelator + const scarab::param_array* offsetPairs = node-> array_at("offsets"); //change based on DevGuide &KTCorrelator if (offsetPairs != NULL) { for (scarab::param_array::const_iterator pairIt = offsetPairs->begin(); pairIt != offsetPairs->end(); ++pairIt) @@ -46,7 +46,7 @@ namespace Katydid return false; } UIntDoublePair pair((*pairIt)->as_array().get_value< unsigned >(0), (*pairIt)->as_array().get_value< double >(1)); - if (fOffsets.size() <= pair.first) fOffsets.resize(pair.first + 1); //keeping fOffsets for now bc it is defined in the .hh + if (fOffsets.size() <= pair.first) fOffsets.resize(pair.first + 1); fOffsets[pair.first] = pair.second; } } @@ -73,7 +73,7 @@ namespace Katydid for (unsigned iBin = 0; iBin < sliceSize; ++iBin) { - timeSeries->SetValue(iBin, fOffsets[iComponent] + timeSeries->GetValue(iBin)); //is this GetValue diff + timeSeries->SetValue(iBin, fOffsets[iComponent] + timeSeries->GetValue(iBin)); } } diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index d461b753e..84f7a0ead 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -33,7 +33,7 @@ namespace Katydid { } - bool KTGaussianNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) //I changed this to scarab::param somewhere else + bool KTGaussianNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) { if (node == NULL) return false; @@ -41,8 +41,8 @@ namespace Katydid input_type mean = node->get_value< input_type >("mean", fRNG.mean()); input_type sigma = node->get_value< input_type >("sigma", fRNG.sigma()); fRNG.param(KTRNGGaussian<>::param_type(mean, sigma)); - fRNG.SetSeed(node->get_value("seed")); - + fRNG.SetSeed(node->get_value< unsigned >("seed")); + return true; } diff --git a/Source/Simulation/KTTSGenerator.cc b/Source/Simulation/KTTSGenerator.cc index 224331ea6..26339db23 100644 --- a/Source/Simulation/KTTSGenerator.cc +++ b/Source/Simulation/KTTSGenerator.cc @@ -18,8 +18,7 @@ #include "KTTimeSeriesFFTW.hh" #include "KTTimeSeriesReal.hh" -//#include "thorax.hh" it appears this isn't defined anywhere? - +//#include "thorax.hh" :obsolete #include @@ -167,7 +166,7 @@ namespace Katydid } - newHeader->SetTimestamp(scarab::get_absolute_time_string()); //this is a membervariableref... whatever that means + newHeader->SetTimestamp(scarab::get_absolute_time_string()); return newHeader; } @@ -196,8 +195,7 @@ namespace Katydid for (unsigned iComponent = 0; iComponent < fNChannels; ++iComponent) { - sliceHeader.SetTimeStamp((uint64_t)(sliceHeader.GetTimeInRun() * (double)CLOCKS_PER_SEC), iComponent); // TODO: change this to 1e3 when switch to usec is made - // changed above to clocks_per_sec bc that was what the builder suggested... how much do we trust it + sliceHeader.SetTimeStamp((uint64_t)(sliceHeader.GetTimeInRun() * (double)NSEC_PER_SEC), iComponent); // TODO: change this to 1e3 when switch to usec is made sliceHeader.SetAcquisitionID(0); sliceHeader.SetRecordID(0); } @@ -235,4 +233,4 @@ namespace Katydid } -} /* namespace Katydid */ +} /* namespace Katydid */ \ No newline at end of file diff --git a/Source/Simulation/KTTSGenerator.hh b/Source/Simulation/KTTSGenerator.hh index babca06c3..9b6e63a06 100644 --- a/Source/Simulation/KTTSGenerator.hh +++ b/Source/Simulation/KTTSGenerator.hh @@ -123,11 +123,9 @@ namespace Katydid // Signals //*************** private: - //Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; this was the old one - Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; //this is the new one according to KTEggProcessor.hh in Source/Time could still be wrong tho - //compare in soruce/time directory in KTeggprocessor make the header signal look the same also check done signal in line129 + Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; Nymph::KTSignalData fDataSignal; - Nymph::KTSignalOneArg< void > fDoneSignal; //this looks the same as last + Nymph::KTSignalOneArg< void > fDoneSignal; Nymph::KTSignalOneArg< const KTProcSummary* > fSummarySignal; }; diff --git a/Source/Time/KTEggProcessor.hh b/Source/Time/KTEggProcessor.hh index 17a11eb10..2af98b7a0 100644 --- a/Source/Time/KTEggProcessor.hh +++ b/Source/Time/KTEggProcessor.hh @@ -1,5 +1,5 @@ /** - @file KTEggprocessor.hh + @file KTEggProcessor.hh @brief Contains KTEggProcessor @details Iterates over slices in an Egg file. @author: N. S. Oblath From 33a521512e90822044efc3bbf5bbd2925085596d Mon Sep 17 00:00:00 2001 From: adil Date: Tue, 27 May 2025 18:53:40 -0400 Subject: [PATCH 07/12] Minor indentation fix --- Source/Simulation/KTGaussianNoiseGenerator.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 84f7a0ead..8bfef84cd 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -78,7 +78,7 @@ namespace Katydid timeSeries->SetValue(iBin, timeSeries->GetValue(iBin) + fRNG()); } } - } + } return true; } From ac7f48445b0b6b6ed0e64c4687a6f182aebc5443 Mon Sep 17 00:00:00 2001 From: adil Date: Thu, 29 May 2025 12:57:27 -0400 Subject: [PATCH 08/12] Initial implementation of cavity noise model under KTCavityNoiseGenerator --- .../ConfigFiles/KatydidPSConfig_noise.yaml | 10 +- .../KatydidPSConfig_noise_cavity.yaml | 58 ++++++ Source/Simulation/CMakeLists.txt | 2 + Source/Simulation/KTCavityNoiseGenerator.cc | 167 ++++++++++++++++++ Source/Simulation/KTCavityNoiseGenerator.hh | 58 ++++++ Source/Simulation/KTGaussianNoiseGenerator.hh | 2 +- 6 files changed, 291 insertions(+), 6 deletions(-) create mode 100644 Examples/ConfigFiles/KatydidPSConfig_noise_cavity.yaml create mode 100644 Source/Simulation/KTCavityNoiseGenerator.cc create mode 100644 Source/Simulation/KTCavityNoiseGenerator.hh diff --git a/Examples/ConfigFiles/KatydidPSConfig_noise.yaml b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml index 3a9f422ac..039e0658c 100644 --- a/Examples/ConfigFiles/KatydidPSConfig_noise.yaml +++ b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml @@ -5,7 +5,7 @@ processor-toolbox: - type: forward-fftw name: fft - type: gaussian-noise-generator - name: noise_gen + name: noise-gen - type: convert-to-power name: to-ps - type: basic-root-writer @@ -15,8 +15,8 @@ processor-toolbox: - signal: "egg:header" slot: "fft:header" - signal: "egg:ts" - slot: "noise_gen:slice" - - signal: "noise_gen:slice" + slot: "noise-gen:slice" + - signal: "noise-gen:slice" slot: "fft:ts-fftw" - signal: "fft:fft" slot: "to-ps:fs-fftw-to-psd" @@ -37,9 +37,9 @@ egg: fft: transform-flag: ESTIMATE -noise_gen: +noise-gen: mean: 0.0 - sigma: 0.02 + sigma: 0.02 # varying this shifts the noise floor seed: 12345 writer: diff --git a/Examples/ConfigFiles/KatydidPSConfig_noise_cavity.yaml b/Examples/ConfigFiles/KatydidPSConfig_noise_cavity.yaml new file mode 100644 index 000000000..fd2e7bfd8 --- /dev/null +++ b/Examples/ConfigFiles/KatydidPSConfig_noise_cavity.yaml @@ -0,0 +1,58 @@ +processor-toolbox: + processors: + - type: egg-processor + name: egg + - type: forward-fftw + name: fft + - type: cavity-noise-generator + name: noise-gen + - type: convert-to-power + name: to-ps + - type: basic-root-writer + name: writer + + connections: + - signal: "egg:header" + slot: "fft:header" + - signal: "egg:ts" + slot: "noise-gen:slice" + - signal: "noise-gen:slice" + slot: "fft:ts-fftw" + - signal: "fft:fft" + slot: "to-ps:fs-fftw-to-psd" + - signal: "to-ps:psd" + slot: "writer:ps" + - signal: "egg:ts" + slot: "writer:ts" + + run-queue: + - egg + +egg: + filename: "locust_mc_Seed600_Angle90.000000000_Pos0.0040000_Energy18563.251000000.egg" + egg-reader: egg3 + slice-size: 4096 + number-of-slices: 1 + +fft: + transform-flag: ESTIMATE + +noise-gen: + seed: 1234 + transform-flag: ESTIMATE + noise-scaling: 5000000.0 # noise amplitude scaling + cavity: # optional + f0 : 25.904e9 + Q_L : 625 + Q0 : 10000 + A : 0.90 + T_line_start : 80.0 + T_line_end : 5.2 + T_cav : 80.0 + T_isol : 5.2 + epsilon : 0.5 + f_lo : 25.9702e9 + +writer: + output-file: "ps_output_noise_cavity.root" + file-flag: recreate diff --git a/Source/Simulation/CMakeLists.txt b/Source/Simulation/CMakeLists.txt index bcbca15e5..62b164fe6 100644 --- a/Source/Simulation/CMakeLists.txt +++ b/Source/Simulation/CMakeLists.txt @@ -4,6 +4,7 @@ set (SIMULATION_HEADERFILES ${CMAKE_CURRENT_SOURCE_DIR}/KTDCOffsetGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTGaussianNoiseGenerator.hh + ${CMAKE_CURRENT_SOURCE_DIR}/KTCavityNoiseGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTSinusoidGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.hh ) @@ -11,6 +12,7 @@ set (SIMULATION_HEADERFILES set (SIMULATION_SOURCEFILES ${CMAKE_CURRENT_SOURCE_DIR}/KTDCOffsetGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTGaussianNoiseGenerator.cc + ${CMAKE_CURRENT_SOURCE_DIR}/KTCavityNoiseGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTSinusoidGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.cc ) diff --git a/Source/Simulation/KTCavityNoiseGenerator.cc b/Source/Simulation/KTCavityNoiseGenerator.cc new file mode 100644 index 000000000..588b84b17 --- /dev/null +++ b/Source/Simulation/KTCavityNoiseGenerator.cc @@ -0,0 +1,167 @@ +/* + * KTCavityNoiseGenerator.cc + * + * Created on: May 28, 2025 + * Author: ehtkarim + */ + +#include "KTCavityNoiseGenerator.hh" + +#include "param.hh" +#include "KTMath.hh" +#include "KTTimeSeriesData.hh" +#include "KTTimeSeries.hh" +#include "KTTimeSeriesFFTW.hh" +#include "KTFrequencySpectrumFFTW.hh" +#include "KTReverseFFTW.hh" + +#include +#include + +using std::string; + +namespace Katydid +{ + KTLOGGER(genlog, "KTCavityNoiseGenerator"); + + KT_REGISTER_PROCESSOR(KTCavityNoiseGenerator, "cavity-noise-generator"); + + KTCavityNoiseGenerator::ModelPars::ModelPars() : + f0(25.904e9), + Q_L(625.0), + Q0(1.e4), + A(0.90), + T_line_start(80.0), + T_line_end(5.2), + T_cav(80.0), + T_isol(5.2), + epsilon(0.5), + f_lo(25.9702e9) + { + } + + KTCavityNoiseGenerator::KTCavityNoiseGenerator(const string& name) : + KTGaussianNoiseGenerator(name), + fPars(), + fTransformFlag("ESTIMATE"), + fNoiseScaling(1.0) + { + } + + KTCavityNoiseGenerator::~KTCavityNoiseGenerator() + { + } + + bool KTCavityNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) + { + if (! KTGaussianNoiseGenerator::ConfigureDerivedGenerator(node)) return false; + if (node == NULL) return false; + + fRNG.param(KTRNGGaussian<>::param_type(0.0, 1.0)); // Cavity noise should have fRNG() with default (mean, sigma), not derived from KTGaussianNoiseGenerator + + if (node->as_node().has("cavity")) + { + const scarab::param_node& m = (*node)["cavity"].as_node(); + #define GET(v) v = m.get_value(#v, v) + GET(fPars.f0); GET(fPars.Q_L); GET(fPars.Q0); GET(fPars.A); + GET(fPars.T_line_start); GET(fPars.T_line_end); + GET(fPars.T_cav); GET(fPars.T_isol); + GET(fPars.epsilon); GET(fPars.f_lo); + #undef GET + } + + fNoiseScaling = node->get_value("noise-scaling", fNoiseScaling); + if (fNoiseScaling <= 0.0) + { + KTWARN(genlog, "\"noise-scaling\" must be > 0; using 1.0"); + fNoiseScaling = 1.0; + } + + fTransformFlag = node->get_value("transform-flag", fTransformFlag); + + return true; + } + + bool KTCavityNoiseGenerator::GenerateTS(KTTimeSeriesData& data) + { + const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); + const unsigned sliceSize = data.GetTimeSeries(0)->GetNTimeBins(); + const unsigned nComponents = data.GetNComponents(); + + const double fs = 1.0 / binWidth; + const double df = fs / sliceSize; + const unsigned N2 = sliceSize / 2; + + KTFrequencySpectrumFFTW spec(sliceSize, -fs*0.5, fs*0.5, false); + spec.SetNTimeBins(sliceSize); + + for (unsigned k = 0; k <= N2; ++k) + { + const double f_if = k * df; + const double f_rf = -f_if + fPars.f_lo; // Down-converted + const double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin + const double amp = fNoiseScaling * std::sqrt(pBin) * N2; + + const double re = amp * fRNG(); + const double im = (k==0 || (sliceSize%2==0 && k==N2)) ? 0.0 : amp * fRNG(); // Set imag component 0 for the DC bin (k = 0) and for the Nyquist bin (k = N/2) (even); otherwise amp * fRNG() + + spec.SetRect(k, re, im); + if (k>0 && k noiseTS( rfft.TransformToComplex(&spec) ); + if (! noiseTS) + { + KTERROR(genlog, "Inverse FFT failed while producing cavity noise"); + return false; + } + + const double norm = 1.0 / sliceSize; // 1/N normalization + + for (unsigned iComponent = 0; iComponent < nComponents; ++iComponent) + { + KTTimeSeries* ts = data.GetTimeSeries(iComponent); + + if (auto* tsFFTW = dynamic_cast< KTTimeSeriesFFTW* >(ts)) + { + for (unsigned i = 0; i < sliceSize; ++i) + tsFFTW->SetRect(i, tsFFTW->GetReal(i) + noiseTS->GetReal(i)*norm, tsFFTW->GetImag(i) + noiseTS->GetImag(i)*norm); + } + else + { + for (unsigned i = 0; i < sliceSize; ++i) + ts->SetValue(i, ts->GetValue(i) + noiseTS->GetReal(i)*norm); + } + } + + return true; + } + + double KTCavityNoiseGenerator::NoisePSD(double f) const + { + const double g = fPars.Q0 / fPars.Q_L; + const double lor = 1.0 / ( 1.0 + std::pow( 2.0 * fPars.Q_L * (f - fPars.f0) / fPars.f0, 2.0 ) ); // Lorentzian + + const double kB = 1.380649e-23; + const double hbar= 1.054571817e-34; + const double omega = 2.0 * M_PI * f; + + auto eta = [](double x){ return x/(std::exp(x)-1.0) + 0.5; }; // Bose-Einstein factor + + const double Tcav = fPars.T_cav * eta(hbar*omega/(kB*fPars.T_cav)); + const double Tisol = fPars.T_isol * eta(hbar*omega/(kB*fPars.T_isol)); + + const double P_cav = kB*Tcav * (4.*g/std::pow(1.+g,2)) * lor; + const double P_loss = kB*( fPars.A*fPars.T_line_start + (1.-fPars.A)*fPars.T_line_end ); + const double P_isol = kB*Tisol * (1. - (4.*g/std::pow(1.+g,2))*lor); + const double P_amp = hbar*omega / fPars.epsilon; + + return P_cav + P_loss + P_isol + P_amp; + } + +} /* namespace Katydid */ diff --git a/Source/Simulation/KTCavityNoiseGenerator.hh b/Source/Simulation/KTCavityNoiseGenerator.hh new file mode 100644 index 000000000..2f4ba64bd --- /dev/null +++ b/Source/Simulation/KTCavityNoiseGenerator.hh @@ -0,0 +1,58 @@ +/* + * KTCavityNoiseGenerator.hh + * + * Created on: May 28, 2025 + * Author: ehtkarim +*/ + +#ifndef KTCAVITYNOISEGENERATOR_HH_ +#define KTCAVITYNOISEGENERATOR_HH_ + +#include "KTGaussianNoiseGenerator.hh" + +#include + +namespace Katydid +{ + + /*! @class KTCavityNoiseGenerator + * @brief Generates cavity noise in a time series + */ + class KTCavityNoiseGenerator : public KTGaussianNoiseGenerator + { + public: + KTCavityNoiseGenerator(const std::string& name = "cavity-noise-generator"); + virtual ~KTCavityNoiseGenerator(); + + virtual bool ConfigureDerivedGenerator(const scarab::param_node* node); + + protected: + struct ModelPars + { + double f0; + double Q_L; + double Q0; + double A; + double T_line_start; + double T_line_end; + double T_cav; + double T_isol; + double epsilon; + double f_lo; + + ModelPars(); + }; + + ModelPars fPars; + std::string fTransformFlag; + double fNoiseScaling; + + double NoisePSD(double f) const; + + public: + virtual bool GenerateTS(KTTimeSeriesData& data); + }; + +} /* namespace Katydid */ + +#endif /* KTCAVITYNOISEGENERATOR_HH_ */ diff --git a/Source/Simulation/KTGaussianNoiseGenerator.hh b/Source/Simulation/KTGaussianNoiseGenerator.hh index fde677024..41e159a36 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.hh +++ b/Source/Simulation/KTGaussianNoiseGenerator.hh @@ -61,7 +61,7 @@ namespace Katydid double GetSigma() const; void SetSigma(double sigma); - private: + protected: KTRNGGaussian<> fRNG; public: From 8b0f5a3c6fcdc6eda6ad244d661fe3946080694a Mon Sep 17 00:00:00 2001 From: adil Date: Mon, 16 Jun 2025 23:28:49 -0400 Subject: [PATCH 09/12] Fixed cavity parameters reading-in, f_lo implementation & complex signal noise handling --- Source/Simulation/KTCavityNoiseGenerator.cc | 57 ++++++++++++++------- 1 file changed, 39 insertions(+), 18 deletions(-) diff --git a/Source/Simulation/KTCavityNoiseGenerator.cc b/Source/Simulation/KTCavityNoiseGenerator.cc index 588b84b17..0b1718191 100644 --- a/Source/Simulation/KTCavityNoiseGenerator.cc +++ b/Source/Simulation/KTCavityNoiseGenerator.cc @@ -54,20 +54,24 @@ namespace Katydid bool KTCavityNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) { - if (! KTGaussianNoiseGenerator::ConfigureDerivedGenerator(node)) return false; if (node == NULL) return false; + if (! KTGaussianNoiseGenerator::ConfigureDerivedGenerator(node)) return false; fRNG.param(KTRNGGaussian<>::param_type(0.0, 1.0)); // Cavity noise should have fRNG() with default (mean, sigma), not derived from KTGaussianNoiseGenerator - if (node->as_node().has("cavity")) + if (node->has("cavity")) { const scarab::param_node& m = (*node)["cavity"].as_node(); - #define GET(v) v = m.get_value(#v, v) - GET(fPars.f0); GET(fPars.Q_L); GET(fPars.Q0); GET(fPars.A); - GET(fPars.T_line_start); GET(fPars.T_line_end); - GET(fPars.T_cav); GET(fPars.T_isol); - GET(fPars.epsilon); GET(fPars.f_lo); - #undef GET + fPars.f0 = m.get_value("f0", fPars.f0); + fPars.Q_L = m.get_value("Q_L", fPars.Q_L); + fPars.Q0 = m.get_value("Q0", fPars.Q0); + fPars.A = m.get_value("A", fPars.A); + fPars.T_line_start = m.get_value("T_line_start", fPars.T_line_start); + fPars.T_line_end = m.get_value("T_line_end", fPars.T_line_end); + fPars.T_cav = m.get_value("T_cav", fPars.T_cav); + fPars.T_isol = m.get_value("T_isol", fPars.T_isol); + fPars.epsilon = m.get_value("epsilon", fPars.epsilon); + fPars.f_lo = m.get_value("f_lo", fPars.f_lo); } fNoiseScaling = node->get_value("noise-scaling", fNoiseScaling); @@ -92,22 +96,39 @@ namespace Katydid const double df = fs / sliceSize; const unsigned N2 = sliceSize / 2; + bool isComplex = dynamic_cast< KTTimeSeriesFFTW* >(data.GetTimeSeries(0)) != NULL; + KTFrequencySpectrumFFTW spec(sliceSize, -fs*0.5, fs*0.5, false); spec.SetNTimeBins(sliceSize); - for (unsigned k = 0; k <= N2; ++k) + if (isComplex) + { + for (unsigned k = 0; k < sliceSize; ++k) + { + double f_if = (k <= N2) ? k * df : (static_cast(k) - static_cast(sliceSize)) * df; + double f_rf = f_if + fPars.f_lo; // Down-converted + double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin + double amp = fNoiseScaling * std::sqrt(pBin) * N2; + + spec.SetRect(k, amp * fRNG(), amp * fRNG()); + } + } + else { - const double f_if = k * df; - const double f_rf = -f_if + fPars.f_lo; // Down-converted - const double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin - const double amp = fNoiseScaling * std::sqrt(pBin) * N2; + for (unsigned k = 0; k <= N2; ++k) + { + double f_if = k * df; + double f_rf = f_if + fPars.f_lo; // Down-converted + double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin + double amp = fNoiseScaling * std::sqrt(pBin) * N2; - const double re = amp * fRNG(); - const double im = (k==0 || (sliceSize%2==0 && k==N2)) ? 0.0 : amp * fRNG(); // Set imag component 0 for the DC bin (k = 0) and for the Nyquist bin (k = N/2) (even); otherwise amp * fRNG() + double re = amp * fRNG(); + double im = (k==0 || (sliceSize%2==0 && k==N2)) ? 0.0 : amp * fRNG(); // Set imag component 0 for the DC bin (k = 0) and for the Nyquist bin (k = N/2) (even); otherwise amp * fRNG() - spec.SetRect(k, re, im); - if (k>0 && k0 && k Date: Wed, 30 Jul 2025 23:54:10 -0400 Subject: [PATCH 10/12] Added feature to set aWGN by noise-floor-psd or noise-temperature, and corrected cavity noise model power from the PSD --- .../ConfigFiles/KatydidPSConfig_noise.yaml | 2 +- Source/Simulation/KTCavityNoiseGenerator.cc | 10 ++- Source/Simulation/KTGaussianNoiseGenerator.cc | 65 +++++++++++++++++-- Source/Simulation/KTGaussianNoiseGenerator.hh | 6 ++ 4 files changed, 70 insertions(+), 13 deletions(-) diff --git a/Examples/ConfigFiles/KatydidPSConfig_noise.yaml b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml index 039e0658c..78f3d56e9 100644 --- a/Examples/ConfigFiles/KatydidPSConfig_noise.yaml +++ b/Examples/ConfigFiles/KatydidPSConfig_noise.yaml @@ -39,7 +39,7 @@ fft: noise-gen: mean: 0.0 - sigma: 0.02 # varying this shifts the noise floor + noise-floor-psd: 2.2e-13 # Noise power in W/Hz seed: 12345 writer: diff --git a/Source/Simulation/KTCavityNoiseGenerator.cc b/Source/Simulation/KTCavityNoiseGenerator.cc index 0b1718191..338e4e83a 100644 --- a/Source/Simulation/KTCavityNoiseGenerator.cc +++ b/Source/Simulation/KTCavityNoiseGenerator.cc @@ -108,7 +108,7 @@ namespace Katydid double f_if = (k <= N2) ? k * df : (static_cast(k) - static_cast(sliceSize)) * df; double f_rf = f_if + fPars.f_lo; // Down-converted double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin - double amp = fNoiseScaling * std::sqrt(pBin) * N2; + double amp = fNoiseScaling*fGain*std::sqrt(fResistance)*std::pow(sliceSize, 1.5)*std::sqrt(pBin / 2.0); // N^{3/2}*sqrt(P_bin/2) - N^{3/2} since ReverseFFTW does a sqrt(N) normalization spec.SetRect(k, amp * fRNG(), amp * fRNG()); } @@ -120,7 +120,7 @@ namespace Katydid double f_if = k * df; double f_rf = f_if + fPars.f_lo; // Down-converted double pBin = NoisePSD(f_rf) * df; // PSD -> power in one FFT bin - double amp = fNoiseScaling * std::sqrt(pBin) * N2; + double amp = fNoiseScaling*fGain*std::sqrt(fResistance)*std::pow(sliceSize, 1.5)*std::sqrt(pBin); // N^{3/2}*sqrt(P_bin) - for real signal bins P_bin/2 -> P_bin double re = amp * fRNG(); double im = (k==0 || (sliceSize%2==0 && k==N2)) ? 0.0 : amp * fRNG(); // Set imag component 0 for the DC bin (k = 0) and for the Nyquist bin (k = N/2) (even); otherwise amp * fRNG() @@ -142,8 +142,6 @@ namespace Katydid return false; } - const double norm = 1.0 / sliceSize; // 1/N normalization - for (unsigned iComponent = 0; iComponent < nComponents; ++iComponent) { KTTimeSeries* ts = data.GetTimeSeries(iComponent); @@ -151,12 +149,12 @@ namespace Katydid if (auto* tsFFTW = dynamic_cast< KTTimeSeriesFFTW* >(ts)) { for (unsigned i = 0; i < sliceSize; ++i) - tsFFTW->SetRect(i, tsFFTW->GetReal(i) + noiseTS->GetReal(i)*norm, tsFFTW->GetImag(i) + noiseTS->GetImag(i)*norm); + tsFFTW->SetRect(i, tsFFTW->GetReal(i) + noiseTS->GetReal(i), tsFFTW->GetImag(i) + noiseTS->GetImag(i)); } else { for (unsigned i = 0; i < sliceSize; ++i) - ts->SetValue(i, ts->GetValue(i) + noiseTS->GetReal(i)*norm); + ts->SetValue(i, ts->GetValue(i) + noiseTS->GetReal(i)); } } diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 8bfef84cd..5f6e67a75 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.cc +++ b/Source/Simulation/KTGaussianNoiseGenerator.cc @@ -3,6 +3,8 @@ * * Created on: May 3, 2013 * Author: nsoblath + * Edited on: Jul 25, 2025 + * Author: ehtkarim */ #include "KTGaussianNoiseGenerator.hh" @@ -14,6 +16,7 @@ #include "KTTimeSeriesFFTW.hh" #include +#include using std::string; @@ -25,7 +28,10 @@ namespace Katydid KTGaussianNoiseGenerator::KTGaussianNoiseGenerator(const string& name) : KTTSGenerator(name), - fRNG() + fRNG(), + fSigmaPSD(0.0), + fGain(1.0), + fResistance(50.0) // (Ω) { } @@ -38,17 +44,64 @@ namespace Katydid if (node == NULL) return false; typedef KTRNGGaussian<>::input_type input_type; - input_type mean = node->get_value< input_type >("mean", fRNG.mean()); - input_type sigma = node->get_value< input_type >("sigma", fRNG.sigma()); + + if (node->has("noise-floor-psd") && node->has("noise-temperature")) + { + KTERROR(genlog, "Both noise-floor-psd and noise-temperature are defined. Only one can be used!"); + return false; + } + + input_type sigma = 0.0; + + if (node->has("noise-floor-psd")) + { + sigma = std::sqrt(node->get_value("noise-floor-psd")); + } + else if (node->has("noise-temperature")) + { + static constexpr double kBoltzmann = 1.38064852e-23; // J/K - Boltzmann constant; there's none defined in Source/Utility, so keeping a local constexpr + sigma = std::sqrt(kBoltzmann * node->get_value("noise-temperature")); + } + else // falling back to the old "sigma" parameter + { + sigma = node->get_value("sigma", fRNG.sigma()); + } + + input_type mean = node->get_value("mean", fRNG.mean()); + + fSigmaPSD = sigma; // to avoid sigma growing in each time-slice fRNG.param(KTRNGGaussian<>::param_type(mean, sigma)); - fRNG.SetSeed(node->get_value< unsigned >("seed")); + + unsigned seed; + if (node-> has("seed")) + { + seed = node->get_value< unsigned >("seed"); + } + else + { + std::random_device rd; + seed = rd(); + } + fRNG.SetSeed(seed); + + fGain = node->get_value("gain", fGain); + fResistance = node->get_value("resistance", fResistance); return true; } bool KTGaussianNoiseGenerator::GenerateTS(KTTimeSeriesData& data) { - const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); + const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); + const double acquisitionRate = 1.0 / binWidth; // (Hz) + static constexpr double kInvSqrt2 = M_SQRT1_2; // sqrt(0.5) - sigma scaling for complex signal + + // Converting the stored PSD‑sigma (V/sqrt(Hz)) into the per‑sample sigma (V) + double sampledSigma = fSigmaPSD * std::sqrt(acquisitionRate); + + sampledSigma *= fGain * std::sqrt(fResistance); // V = sigma x Gain x sqrt(R) + fRNG.param(KTRNGGaussian<>::param_type(fRNG.mean(), sampledSigma)); // Updating the RNG + const unsigned sliceSize = data.GetTimeSeries(0)->GetNTimeBins(); unsigned nComponents = data.GetNComponents(); @@ -68,7 +121,7 @@ namespace Katydid { for (unsigned iBin = 0; iBin < sliceSize; ++iBin) { - tsFFTW->SetRect(iBin, tsFFTW->GetReal(iBin) + fRNG(), tsFFTW->GetImag(iBin) + fRNG()); // Complex white-Gaussian noise + tsFFTW->SetRect(iBin, tsFFTW->GetReal(iBin) + fRNG() * kInvSqrt2, tsFFTW->GetImag(iBin) + fRNG() * kInvSqrt2); // Complex white-Gaussian noise } } else diff --git a/Source/Simulation/KTGaussianNoiseGenerator.hh b/Source/Simulation/KTGaussianNoiseGenerator.hh index 41e159a36..4b5658306 100644 --- a/Source/Simulation/KTGaussianNoiseGenerator.hh +++ b/Source/Simulation/KTGaussianNoiseGenerator.hh @@ -36,8 +36,11 @@ namespace Katydid - "time-series-type": string -- Type of time series to produce (options: real [default], fftw) - "record-size": unsigned -- Size of the imaginary record that this slice came from (only used to fill in the egg header; does not affect the simulation at all) - From KTGaussianNoiseGenerator + - "seed": int -- Seed used to generate random noise. If this is omitted then the noise spectrum is irreproducible - "mean": double -- Mean for the randomly-chosen time-series values - "sigma": double -- Standard deviation for the randomly-chosen time-series values + - "noise-floor-psd": double -- Noise power in W/Hz + - "noise-temperature" - double -- Noise temperature in K Slots: (inherited from KTTSGenerator) - "slice": void (Nymph::KTDataPtr) -- Add a signal to an existing time series; Requires KTTimeSeriesData; Emits signal "slice" when done. @@ -63,6 +66,9 @@ namespace Katydid protected: KTRNGGaussian<> fRNG; + double fSigmaPSD; + double fGain; + double fResistance; public: virtual bool GenerateTS(KTTimeSeriesData& data); From 62dcc257d2cb2341ec0b5fcd002b26d979c1c923 Mon Sep 17 00:00:00 2001 From: adil Date: Wed, 1 Oct 2025 13:41:39 -0400 Subject: [PATCH 11/12] Added documentation for new processors, and created 'KTMeasuredNoiseGenerator' processor. Updated 'CMakeLists.txt' --- Source/Simulation/CMakeLists.txt | 2 + Source/Simulation/KTCavityNoiseGenerator.cc | 2 +- Source/Simulation/KTCavityNoiseGenerator.hh | 44 +- Source/Simulation/KTMeasuredNoiseGenerator.cc | 592 ++++++++++++++++++ Source/Simulation/KTMeasuredNoiseGenerator.hh | 125 ++++ 5 files changed, 761 insertions(+), 4 deletions(-) create mode 100644 Source/Simulation/KTMeasuredNoiseGenerator.cc create mode 100644 Source/Simulation/KTMeasuredNoiseGenerator.hh diff --git a/Source/Simulation/CMakeLists.txt b/Source/Simulation/CMakeLists.txt index 62b164fe6..e99c78d09 100644 --- a/Source/Simulation/CMakeLists.txt +++ b/Source/Simulation/CMakeLists.txt @@ -5,6 +5,7 @@ set (SIMULATION_HEADERFILES ${CMAKE_CURRENT_SOURCE_DIR}/KTDCOffsetGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTGaussianNoiseGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTCavityNoiseGenerator.hh + ${CMAKE_CURRENT_SOURCE_DIR}/KTMeasuredNoiseGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTSinusoidGenerator.hh ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.hh ) @@ -13,6 +14,7 @@ set (SIMULATION_SOURCEFILES ${CMAKE_CURRENT_SOURCE_DIR}/KTDCOffsetGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTGaussianNoiseGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTCavityNoiseGenerator.cc + ${CMAKE_CURRENT_SOURCE_DIR}/KTMeasuredNoiseGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTSinusoidGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTTSGenerator.cc ) diff --git a/Source/Simulation/KTCavityNoiseGenerator.cc b/Source/Simulation/KTCavityNoiseGenerator.cc index 338e4e83a..ab452ccc7 100644 --- a/Source/Simulation/KTCavityNoiseGenerator.cc +++ b/Source/Simulation/KTCavityNoiseGenerator.cc @@ -57,7 +57,7 @@ namespace Katydid if (node == NULL) return false; if (! KTGaussianNoiseGenerator::ConfigureDerivedGenerator(node)) return false; - fRNG.param(KTRNGGaussian<>::param_type(0.0, 1.0)); // Cavity noise should have fRNG() with default (mean, sigma), not derived from KTGaussianNoiseGenerator + fRNG.param(KTRNGGaussian<>::param_type(0.0, 1.0)); // Cavity noise should have fRNG() with default (mean, sigma), not inherited from KTGaussianNoiseGenerator if (node->has("cavity")) { diff --git a/Source/Simulation/KTCavityNoiseGenerator.hh b/Source/Simulation/KTCavityNoiseGenerator.hh index 2f4ba64bd..a34b75d37 100644 --- a/Source/Simulation/KTCavityNoiseGenerator.hh +++ b/Source/Simulation/KTCavityNoiseGenerator.hh @@ -15,9 +15,47 @@ namespace Katydid { - /*! @class KTCavityNoiseGenerator - * @brief Generates cavity noise in a time series - */ + /*! + @class KTCavityNoiseGenerator + @author E. Karim (Adil) + + @brief Generates a time series with Gaussian noise using the frequency-dependent cavity noise model + + @details + Can create a new time series and drive processing, or can add Gaussian noise to an existing time series. + + Available configuration options: + - Inherited from KTTSGenerator + - "number-of-slices": unsigned -- Number of slices to create (used only if creating new slices) + - "n-channels": unsigned -- Number of channels per slice to create (used only if creating new slices) + - "slice-size": unsigned -- Specify the size of the time series (used only if creating new slices) + - "bin-width": double -- Specify the bin width + - "time-series-type": string -- Type of time series to produce (options: real [default], fftw) + - "record-size": unsigned -- Size of the imaginary record that this slice came from (only used to fill in the egg header; does not affect the simulation at all) + - Inherited from KTGaussianNoiseGenerator + - "seed": int -- Seed used to generate random noise. If this is omitted then the noise spectrum is irreproducible + - From KTCavityNoiseGenerator + - "noise-scaling": double -- Scaling factor of noise amplitude, must be positive-definite. Defaults to unity + - "cavity": node -- Node for defining parameters of the cavity noise model + - "f0": -- Cavity Resonance Frequency (in Hz) + - "Q_L": -- Loaded Cavity Q-factor + - "Q0": -- Unloaded Cavity Q-factor + - "A": -- Transmission line attenuation factor. Can take values from 0 to 1.0 + - "T_line_start": -- Temperature at the cavity end of the transmission line (in K) + - "T_line_end": -- Temperature at the isolator end of the transmission line (in K) + - "T_cav": -- Physical temperature of the cavity (in K) + - "T_isol": -- Physical Temperature of the isolator (in K) + - "epsilon": -- Efficiency factor for the Quantum-Amplifier. Can take values from 0 to 1.0 + - "f_lo": -- Local Oscillator Frequency in the DAQ chain (in Hz) + + Slots: (inherited from KTTSGenerator) + - "slice": void (Nymph::KTDataPtr) -- Add a signal to an existing time series; Requires KTTimeSeriesData; Emits signal "slice" when done. + + Signals: (inherited from KTTSGenerator) + - "header": void (KTEggHeader*) -- emitted when the egg header is created. + - "slice": void (Nymph::KTDataPtr) -- emitted when the new time series is produced or processed. + - "done": void () -- emitted when the job is complete. + */ class KTCavityNoiseGenerator : public KTGaussianNoiseGenerator { public: diff --git a/Source/Simulation/KTMeasuredNoiseGenerator.cc b/Source/Simulation/KTMeasuredNoiseGenerator.cc new file mode 100644 index 000000000..43d8103e9 --- /dev/null +++ b/Source/Simulation/KTMeasuredNoiseGenerator.cc @@ -0,0 +1,592 @@ +/* + * KTMeasuredNoiseGenerator.cc + * + * Created on: Sep 17, 2025 + * Author: ehtkarim + */ + +#include "KTMeasuredNoiseGenerator.hh" + +#include "param.hh" +#include "param_yaml.hh" +#include "KTMath.hh" +#include "KTTimeSeriesData.hh" +#include "KTTimeSeries.hh" +#include "KTTimeSeriesFFTW.hh" +#include "KTFrequencySpectrumFFTW.hh" +#include "KTReverseFFTW.hh" + +#include +#include +#include + +using std::string; + +namespace Katydid +{ + KTLOGGER(genlog, "KTMeasuredNoiseGenerator"); + + KT_REGISTER_PROCESSOR(KTMeasuredNoiseGenerator, "measured-noise-generator"); + + // ---------- KCubicSpline ---------- + + KTMeasuredNoiseGenerator::KCubicSpline::KCubicSpline() : + fX(), fA(), fB(), fC(), fD() + {} + + KTMeasuredNoiseGenerator::KCubicSpline::~KCubicSpline() = default; + + void KTMeasuredNoiseGenerator::KCubicSpline::Clear() + { + fX.clear(); fA.clear(); fB.clear(); fC.clear(); fD.clear(); + } + + bool KTMeasuredNoiseGenerator::KCubicSpline::IsBuilt() const + { + return ! fX.empty(); + } + + bool KTMeasuredNoiseGenerator::KCubicSpline::Build(const std::vector< double >& x, const std::vector< double >& y) + { + Clear(); + + if (x.size() < 2 || x.size() != y.size()) return false; + + for (size_t i = 1; i < x.size(); ++i) + { + if (x[i] <= x[i-1]) return false; + } + + size_t n = x.size(); + fX = x; + fA = y; + fB.assign(n - 1, 0.0); + fC.assign(n, 0.0); + fD.assign(n - 1, 0.0); + + std::vector< double > h(n - 1, 0.0); + for (size_t i = 0; i < n - 1; ++i) h[i] = fX[i+1] - fX[i]; + + std::vector< double > alpha(n, 0.0); + for (size_t i = 1; i < n - 1; ++i) + { + alpha[i] = 3.0/h[i]*(fA[i+1]-fA[i]) - 3.0/h[i-1]*(fA[i]-fA[i-1]); + } + + std::vector< double > l(n, 0.0), mu(n, 0.0), z(n, 0.0); + l[0] = 1.0; + mu[0] = 0.0; + z[0] = 0.0; + + for (size_t i = 1; i < n - 1; ++i) + { + l[i] = 2.0*(fX[i+1]-fX[i-1]) - h[i-1]*mu[i-1]; + if (l[i] == 0.0) return false; + mu[i] = h[i]/l[i]; + z[i] = (alpha[i] - h[i-1]*z[i-1]) / l[i]; + } + + l[n-1] = 1.0; + z[n-1] = 0.0; + fC[n-1] = 0.0; + + for (int j = static_cast(n) - 2; j >= 0; --j) + { + fC[j] = z[j] - mu[j]*fC[j+1]; + fB[j] = (fA[j+1]-fA[j])/h[j] - h[j]*(fC[j+1] + 2.0*fC[j])/3.0; + fD[j] = (fC[j+1]-fC[j]) / (3.0*h[j]); + } + + return true; + } + + double KTMeasuredNoiseGenerator::KCubicSpline::Evaluate(double xval) const + { + if (fX.empty()) return 0.0; + + if (xval <= fX.front()) return fA.front(); + if (xval >= fX.back()) return fA.back(); + + std::vector::const_iterator it = std::upper_bound(fX.begin(), fX.end(), xval); + size_t j = static_cast(std::distance(fX.begin(), it) - 1); + + double dx = xval - fX[j]; + return fA[j] + fB[j]*dx + fC[j]*dx*dx + fD[j]*dx*dx*dx; + } + + // ---------- KTMeasuredNoiseGenerator ---------- + + KTMeasuredNoiseGenerator::KTMeasuredNoiseGenerator(const string& name) : + KTGaussianNoiseGenerator(name), + fTransformFlag("ESTIMATE"), + fNoiseScaling(1.0), + fUseVariance(true), + fFixedRadius(true), + fFreqScale(1.0e6), + fFreqMinHz(0.0), + fFreqMaxHz(400.0e6), + fMeanXHz(), + fMeanY(), + fVarY(), + fMeanSpline(), + fVarSpline(), + fHaveMean(false), + fHaveVar(false) + {} + + KTMeasuredNoiseGenerator::~KTMeasuredNoiseGenerator() = default; + + bool KTMeasuredNoiseGenerator::LoadPointsFromArray(const scarab::param_array& arr, std::vector< double >& xs, std::vector< double >& ys, double freqScale) const + { + xs.clear(); + ys.clear(); + + for (unsigned i = 0; i < arr.size(); ++i) + { + double f_in = 0.0; + double v = 0.0; + + if (arr[i].is_array()) + { + const scarab::param_array& pair = arr[i].as_array(); + if (pair.size() < 2 || ! pair[0].is_value() || ! pair[1].is_value()) + { + KTERROR(genlog, "Each point must be a two-element array [frequency, value]"); + return false; + } + f_in = pair[0].as_value().as_double(); + v = pair[1].as_value().as_double(); + } + else if (arr[i].is_node()) + { + const scarab::param_node& pn = arr[i].as_node(); + f_in = pn.get_value("f", 0.0); + if (pn.has("value")) v = pn["value"].as_value().as_double(); + else if (pn.has("psd")) v = pn["psd"].as_value().as_double(); + else if (pn.has("var")) v = pn["var"].as_value().as_double(); + else + { + KTERROR(genlog, "Point node must contain \"value\", \"psd\", or \"var\""); + return false; + } + } + else + { + KTERROR(genlog, "Unexpected YAML element type in measured spectrum array"); + return false; + } + + const double f_hz = f_in * freqScale; + + if (f_hz < fFreqMinHz || f_hz > fFreqMaxHz) + { + KTERROR(genlog, "Frequency value " << f_in << " (scaled to " << f_hz << " Hz) is outside the allowed range [0, 400 MHz]"); + return false; + } + + xs.push_back(f_hz); + ys.push_back(v); + } + + if (xs.size() < 2) + { + KTERROR(genlog, "At least two points are required to build a spline"); + return false; + } + + std::vector< size_t > order(xs.size()); + for (size_t i = 0; i < order.size(); ++i) order[i] = i; + std::sort(order.begin(), order.end(), [&](size_t a, size_t b){ return xs[a] < xs[b]; }); + + std::vector< double > xs_sorted(xs.size()), ys_sorted(xs.size()); + for (size_t i = 0; i < order.size(); ++i) + { + xs_sorted[i] = xs[order[i]]; + ys_sorted[i] = ys[order[i]]; + } + for (size_t i = 1; i < xs_sorted.size(); ++i) + { + if (xs_sorted[i] <= xs_sorted[i-1]) + { + KTERROR(genlog, "Input frequencies must be strictly increasing"); + return false; + } + } + + xs.swap(xs_sorted); + ys.swap(ys_sorted); + return true; + } + + bool KTMeasuredNoiseGenerator::BuildSplines() + { + if (! fHaveMean) + { + KTERROR(genlog, "Mean PSD points were not provided"); + return false; + } + if (! fMeanSpline.Build(fMeanXHz, fMeanY)) + { + KTERROR(genlog, "Failed to build mean PSD spline"); + return false; + } + + if (fHaveVar) + { + if (! fVarSpline.Build(fMeanXHz, fVarY)) + { + KTERROR(genlog, "Failed to build variance PSD spline"); + return false; + } + } + + return true; + } + + bool KTMeasuredNoiseGenerator::ConfigureDerivedGenerator(const scarab::param_node* node) + { + if (node == NULL) return false; + if (! KTGaussianNoiseGenerator::ConfigureDerivedGenerator(node)) return false; + + fRNG.param(KTRNGGaussian<>::param_type(0.0, 1.0)); + + fNoiseScaling = node->get_value("noise-scaling", fNoiseScaling); + if (fNoiseScaling <= 0.0) + { + KTWARN(genlog, "\"noise-scaling\" must be > 0; using 1.0"); + fNoiseScaling = 1.0; + } + + fTransformFlag = node->get_value("transform-flag", fTransformFlag); + fUseVariance = node->get_value("use-variance", fUseVariance); + fFixedRadius = node->get_value("fixed-radius", fFixedRadius); + + // Default units for values provided in *this* config + string units = node->get_value("frequency-units", string("MHz")); + if (units == "Hz") fFreqScale = 1.0; + else if (units == "kHz") fFreqScale = 1.0e3; + else if (units == "MHz") fFreqScale = 1.0e6; + else if (units == "GHz") fFreqScale = 1.0e9; + else + { + KTWARN(genlog, "Unrecognized frequency-units \"" << units << "\"; assuming MHz"); + fFreqScale = 1.0e6; + } + + if (! node->has("measured-noise")) + { + KTERROR(genlog, "Missing required \"measured-noise\" configuration node"); + return false; + } + + const scarab::param_node& mn = (*node)["measured-noise"].as_node(); + + // ---- loading arrays from an external YAML file if provided ---- + string yamlFile; + if (mn.has("psd-file")) yamlFile = mn.get_value("psd-file", string()); + else if (mn.has("psd-values-file")) yamlFile = mn.get_value("psd-values-file", string()); + else if (mn.has("psd-file-dir")) + { + string dir = mn.get_value("psd-file-dir", string()); + string name = mn.get_value("psd-file-name", string()); + if (name.empty()) + { + KTERROR(genlog, "\"psd-file-dir\" provided but \"psd-file-name\" is missing"); + return false; + } + yamlFile = dir; + if (! yamlFile.empty() && yamlFile.back() != '/' ) yamlFile += "/"; + yamlFile += name; + } + + if (! yamlFile.empty()) + { + scarab::param_input_yaml reader; + auto doc = reader.read_file(yamlFile); + if (!doc || !doc->is_node()) + { + KTERROR(genlog, "Failed to read YAML file \"" << yamlFile << "\""); + return false; + } + const scarab::param_node& root = doc->as_node(); + const scarab::param_node* src = &root; + if (root.has("measured-noise") && root["measured-noise"].is_node()) + src = &root["measured-noise"].as_node(); + + // Allowing the file to override frequency units (for values in the file) + double fileFreqScale = fFreqScale; + if (src->has("frequency-units")) + { + string fu = src->get_value("frequency-units", string("MHz")); + if (fu == "Hz") fileFreqScale = 1.0; + else if (fu == "kHz") fileFreqScale = 1.0e3; + else if (fu == "MHz") fileFreqScale = 1.0e6; + else if (fu == "GHz") fileFreqScale = 1.0e9; + else + { + KTERROR(genlog, "Unrecognized frequency-units in PSD file: \"" << fu << "\""); + return false; + } + } + + // PSD units must be W/Hz (default) + string psdUnits = src->get_value("psd-units", string("W/Hz")); + if (psdUnits != "W/Hz") + { + KTERROR(genlog, "Unsupported psd-units \"" << psdUnits + << "\"; only W/Hz is accepted"); + return false; + } + + if (! src->has("psd-mean") || ! (*src)["psd-mean"].is_array()) + { + KTERROR(genlog, "PSD file is missing required array \"psd-mean\""); + return false; + } + + const scarab::param_array& meanArr = (*src)["psd-mean"].as_array(); + fHaveMean = LoadPointsFromArray(meanArr, fMeanXHz, fMeanY, fileFreqScale); + if (! fHaveMean) return false; + + fHaveVar = false; + if (src->has("psd-variance")) + { + if (! (*src)["psd-variance"].is_array()) + { + KTERROR(genlog, "\"psd-variance\" must be an array of [frequency, value] pairs"); + return false; + } + const scarab::param_array& varArr = (*src)["psd-variance"].as_array(); + std::vector< double > varXHz, varVals; + if (! LoadPointsFromArray(varArr, varXHz, varVals, fileFreqScale)) return false; + + if (varXHz.size() != fMeanXHz.size()) + { + KTERROR(genlog, "Variance and mean arrays must use the same frequency knots"); + return false; + } + for (size_t i = 0; i < varXHz.size(); ++i) + { + if (std::abs(varXHz[i] - fMeanXHz[i]) > 0.5) + { + KTERROR(genlog, "Variance and mean arrays must share identical frequency knots"); + return false; + } + } + + for (double v : varVals) + { + if (v < 0.0) + { + KTERROR(genlog, "Variance values must be non-negative"); + return false; + } + } + + fVarY = varVals; + fHaveVar= true; + } + else + { + KTWARN(genlog, "PSD file has no \"psd-variance\"; will inject using mean only"); + } + } + else + { + // ---- Backward-compatible path: arrays inline in main config ---- + if (! mn.has("psd-mean")) + { + KTERROR(genlog, "Missing required array \"measured-noise.psd-mean\""); + return false; + } + if (! mn["psd-mean"].is_array()) + { + KTERROR(genlog, "\"measured-noise.psd-mean\" must be an array of [frequency, value] pairs"); + return false; + } + const scarab::param_array& meanArr = mn["psd-mean"].as_array(); + fHaveMean = LoadPointsFromArray(meanArr, fMeanXHz, fMeanY, fFreqScale); + if (! fHaveMean) return false; + + fHaveVar = false; + if (mn.has("psd-variance")) + { + if (! mn["psd-variance"].is_array()) + { + KTERROR(genlog, "\"measured-noise.psd-variance\" must be an array of [frequency, value] pairs"); + return false; + } + const scarab::param_array& varArr = mn["psd-variance"].as_array(); + std::vector< double > varXHz, varVals; + if (! LoadPointsFromArray(varArr, varXHz, varVals, fFreqScale)) return false; + + if (varXHz.size() != fMeanXHz.size()) + { + KTERROR(genlog, "Variance and mean arrays must use the same frequency knots"); + return false; + } + for (size_t i = 0; i < varXHz.size(); ++i) + { + if (std::abs(varXHz[i] - fMeanXHz[i]) > 0.5) + { + KTERROR(genlog, "Variance and mean arrays must share identical frequency knots"); + return false; + } + } + + for (double v : varVals) + { + if (v < 0.0) + { + KTERROR(genlog, "Variance values must be non-negative"); + return false; + } + } + + fVarY = varVals; + fHaveVar = true; + } + else + { + KTWARN(genlog, "No \"psd-variance\" provided; will inject using mean only"); + } + } + + if (! BuildSplines()) return false; + + return true; + } + + double KTMeasuredNoiseGenerator::DrawPSD(double f_abs_hz) + { + const double mean_psd = std::max(0.0, fMeanSpline.Evaluate(f_abs_hz)); + + if (! fUseVariance || ! fVarSpline.IsBuilt()) return mean_psd; + + const double var_psd = std::max(0.0, fVarSpline.Evaluate(f_abs_hz)); + const double stdev_psd = std::sqrt(var_psd); + + const double draw = mean_psd + stdev_psd * fRNG(); + return (draw > 0.0) ? draw : 0.0; + } + + void KTMeasuredNoiseGenerator::RandomUnitComplex(double& c, double& s) + { + const double x = fRNG(); + const double y = fRNG(); + const double r = std::sqrt(x*x + y*y); + if (r > 0.0) { c = x / r; s = y / r; } + else { c = 1.0; s = 0.0; } + } + + bool KTMeasuredNoiseGenerator::GenerateTS(KTTimeSeriesData& data) + { + const double binWidth = data.GetTimeSeries(0)->GetTimeBinWidth(); + const unsigned sliceSize = data.GetTimeSeries(0)->GetNTimeBins(); + const unsigned nComponents = data.GetNComponents(); + + const double fs = 1.0 / binWidth; + const double df = fs / sliceSize; + const unsigned N2 = sliceSize / 2; + + bool isComplex = dynamic_cast< KTTimeSeriesFFTW* >(data.GetTimeSeries(0)) != NULL; + + KTFrequencySpectrumFFTW spec(sliceSize, -fs*0.5, fs*0.5, false); + spec.SetNTimeBins(sliceSize); + + const double scale = fNoiseScaling * fGain * std::sqrt(fResistance) * std::pow(sliceSize, 1.5); + + if (isComplex) + { + for (unsigned k = 0; k < sliceSize; ++k) + { + const double f_if = (k <= N2) ? k * df : (static_cast(k) - static_cast(sliceSize)) * df; + const double f_abs = std::fabs(f_if); + + const double psd = DrawPSD(f_abs); + const double pBin = psd * df; + + if (fFixedRadius) + { + const double R = scale * std::sqrt(pBin); + double c = 1.0, s = 0.0; RandomUnitComplex(c, s); + spec.SetRect(k, R * c, R * s); + } + else + { + const double amp = scale * std::sqrt(pBin / 2.0); + spec.SetRect(k, amp * fRNG(), amp * fRNG()); + } + } + } + else + { + for (unsigned k = 0; k <= N2; ++k) + { + const double f_if = k * df; + const double psd = DrawPSD(f_if); + const double pBin = psd * df; + + if (fFixedRadius) + { + if (k==0 || (sliceSize%2==0 && k==N2)) + { + const double R0 = scale * std::sqrt(pBin); + const double sign = (fRNG() >= 0.0) ? 1.0 : -1.0; + spec.SetRect(k, sign * R0, 0.0); + } + else + { + const double R = scale * std::sqrt(2.0 * pBin); + double c = 1.0, s = 0.0; RandomUnitComplex(c, s); + const double re = R * c; + const double im = R * s; + + spec.SetRect(k, re, im); + spec.SetRect(sliceSize - k, re, -im); + } + } + else + { + const double amp = scale * std::sqrt(pBin); + const double re = amp * fRNG(); + const double im = (k==0 || (sliceSize%2==0 && k==N2)) ? 0.0 : amp * fRNG(); + + spec.SetRect(k, re, im); + if (k>0 && k noiseTS( rfft.TransformToComplex(&spec) ); + if (! noiseTS) + { + KTERROR(genlog, "Inverse FFT failed while producing measured noise"); + return false; + } + + for (unsigned iComponent = 0; iComponent < nComponents; ++iComponent) + { + KTTimeSeries* ts = data.GetTimeSeries(iComponent); + + if (auto* tsFFTW = dynamic_cast< KTTimeSeriesFFTW* >(ts)) + { + for (unsigned i = 0; i < sliceSize; ++i) + tsFFTW->SetRect(i, tsFFTW->GetReal(i) + noiseTS->GetReal(i), tsFFTW->GetImag(i) + noiseTS->GetImag(i)); + } + else + { + for (unsigned i = 0; i < sliceSize; ++i) + ts->SetValue(i, ts->GetValue(i) + noiseTS->GetReal(i)); + } + } + + return true; + } + +} /* namespace Katydid */ diff --git a/Source/Simulation/KTMeasuredNoiseGenerator.hh b/Source/Simulation/KTMeasuredNoiseGenerator.hh new file mode 100644 index 000000000..f98fa261a --- /dev/null +++ b/Source/Simulation/KTMeasuredNoiseGenerator.hh @@ -0,0 +1,125 @@ +/* + * KTMeasuredNoiseGenerator.hh + * + * Created on: Sep 17, 2025 + * Author: ehtkarim + */ + +#ifndef KTMEASUREDNOISEGENERATOR_HH_ +#define KTMEASUREDNOISEGENERATOR_HH_ + +#include "KTGaussianNoiseGenerator.hh" + +#include +#include + +namespace scarab { + class param_node; + class param_array; +} + +namespace Katydid +{ + /*! + @class KTMeasuredNoiseGenerator + @author E. Karim (Adil) + + @brief Generates additive noise using user-provided (measured) PSD mean and variance + + @details + Can create a new time series and drive processing, or can add a user-provided (measured) noise spectrum to an existing time series. + + Available configuration options: + - Inherited from KTTSGenerator + - "number-of-slices": unsigned -- Number of slices to create (used only if creating new slices) + - "n-channels": unsigned -- Number of channels per slice to create (used only if creating new slices) + - "slice-size": unsigned -- Specify the size of the time series (used only if creating new slices) + - "bin-width": double -- Specify the bin width + - "time-series-type": string -- Type of time series to produce (options: real [default], fftw) + - "record-size": unsigned -- Size of the imaginary record that this slice came from (only used to fill in the egg header; does not affect the simulation at all) + - Inherited from KTCavityNoiseGenerator + - "noise-scaling": double -- Scaling factor of noise amplitude, must be positive-definite. Defaults to unity + - From KTMeasuredNoiseGenerator + - "frequency-units": string -- Unit of the input psd-array frequency. Can be either "Hz", "kHz", "MHz" or "GHz". If the psd-array input file has units defined, this is overriden + - "fixed-radius": bool -- Mode for which power in each bin is set to exactly match the per-bin power in the input file. If false, complex Gaussian drawing causes additional variance + - "use-variance": bool -- Mode for which the variance in each bin is set to exactly match the per-bin variance in the input file; phase is random and uniform. If false, the resulting spectrum is stochastic with unknown variance per slice (default) + - "measured-noise": node -- Node for defining the file directory of the input psd-arrays + - "psd-file": path -- Directs to the path where the .yaml file is saved with both psd-mean and psd-variance arrays + + Slots: (inherited from KTTSGenerator) + - "slice": void (Nymph::KTDataPtr) -- Add a signal to an existing time series; Requires KTTimeSeriesData; Emits signal "slice" when done. + + Signals: (inherited from KTTSGenerator) + - "header": void (KTEggHeader*) -- emitted when the egg header is created. + - "slice": void (Nymph::KTDataPtr) -- emitted when the new time series is produced or processed. + - "done": void () -- emitted when the job is complete. + */ + class KTMeasuredNoiseGenerator : public KTGaussianNoiseGenerator + { + public: + KTMeasuredNoiseGenerator(const std::string& name = "measured-noise-generator"); + virtual ~KTMeasuredNoiseGenerator(); + + virtual bool ConfigureDerivedGenerator(const scarab::param_node* node); + + public: + virtual bool GenerateTS(KTTimeSeriesData& data); + + protected: + // Simple natural cubic spline for doubles + class KCubicSpline + { + public: + KCubicSpline(); + ~KCubicSpline(); + + void Clear(); + bool Build(const std::vector< double >& x, const std::vector< double >& y); + double Evaluate(double xval) const; + bool IsBuilt() const; + + private: + std::vector< double > fX; + std::vector< double > fA; + std::vector< double > fB; + std::vector< double > fC; + std::vector< double > fD; + }; + + // Helpers for reading arrays and building splines + bool LoadPointsFromArray(const scarab::param_array& arr, std::vector< double >& xs, std::vector< double >& ys, double freqScale) const; + + bool BuildSplines(); + + // Drawing a PSD (W/Hz) sample at |f| using mean and variance splines + double DrawPSD(double f_abs_hz); + + // Random unit complex (cos, sin) using 2D Gaussians; avoids a new RNG + void RandomUnitComplex(double& c, double& s); + + protected: + // Config + std::string fTransformFlag; // FFTW estimate/measure flag + double fNoiseScaling; // global scaling > 0 + bool fUseVariance; // if true, draw bin powers using provided variance + bool fFixedRadius; // if true, use fixed-radius, random-phase (default true) + double fFreqScale; // units scale for input frequencies (Hz multiplier) + double fFreqMinHz; // required input min frequency (Hz) for validation + double fFreqMaxHz; // required input max frequency (Hz) for validation + + // Spline data + std::vector< double > fMeanXHz; + std::vector< double > fMeanY; // W/Hz + std::vector< double > fVarY; // (W/Hz)^2 + + KCubicSpline fMeanSpline; + KCubicSpline fVarSpline; + + bool fHaveMean; + bool fHaveVar; + }; + +} /* namespace Katydid */ + +#endif /* KTMEASUREDNOISEGENERATOR_HH_ */ + From 4a93ba505cc23c4d29d67750f3a47064705b949a Mon Sep 17 00:00:00 2001 From: adil Date: Fri, 10 Oct 2025 23:19:01 -0400 Subject: [PATCH 12/12] Correcting power normalization KTMeasuredNoiseGenerator.cc --- Source/Simulation/KTMeasuredNoiseGenerator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Simulation/KTMeasuredNoiseGenerator.cc b/Source/Simulation/KTMeasuredNoiseGenerator.cc index 43d8103e9..2e7ad8ec9 100644 --- a/Source/Simulation/KTMeasuredNoiseGenerator.cc +++ b/Source/Simulation/KTMeasuredNoiseGenerator.cc @@ -494,7 +494,7 @@ namespace Katydid KTFrequencySpectrumFFTW spec(sliceSize, -fs*0.5, fs*0.5, false); spec.SetNTimeBins(sliceSize); - const double scale = fNoiseScaling * fGain * std::sqrt(fResistance) * std::pow(sliceSize, 1.5); + const double scale = fNoiseScaling * fGain * std::sqrt(fResistance) * sliceSize; if (isComplex) { @@ -537,7 +537,7 @@ namespace Katydid } else { - const double R = scale * std::sqrt(2.0 * pBin); + const double R = scale * std::sqrt(pBin); double c = 1.0, s = 0.0; RandomUnitComplex(c, s); const double re = R * c; const double im = R * s;