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/KatydidPSNoiseCavityConfig.yaml b/Examples/ConfigFiles/KatydidPSNoiseCavityConfig.yaml new file mode 100644 index 000000000..fd2e7bfd8 --- /dev/null +++ b/Examples/ConfigFiles/KatydidPSNoiseCavityConfig.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/Examples/ConfigFiles/KatydidPSNoiseConfig.yaml b/Examples/ConfigFiles/KatydidPSNoiseConfig.yaml new file mode 100644 index 000000000..78f3d56e9 --- /dev/null +++ b/Examples/ConfigFiles/KatydidPSNoiseConfig.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 + noise-floor-psd: 2.2e-13 # Noise power in W/Hz + seed: 12345 + +writer: + output-file: "ps_output_noise.root" + file-flag: recreate diff --git a/Source/Simulation/CMakeLists.txt b/Source/Simulation/CMakeLists.txt index be8f37b0a..e99c78d09 100644 --- a/Source/Simulation/CMakeLists.txt +++ b/Source/Simulation/CMakeLists.txt @@ -4,6 +4,8 @@ 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 ) @@ -11,20 +13,22 @@ 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}/KTMeasuredNoiseGenerator.cc ${CMAKE_CURRENT_SOURCE_DIR}/KTSinusoidGenerator.cc ${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/KTCavityNoiseGenerator.cc b/Source/Simulation/KTCavityNoiseGenerator.cc new file mode 100644 index 000000000..5c63c343a --- /dev/null +++ b/Source/Simulation/KTCavityNoiseGenerator.cc @@ -0,0 +1,183 @@ +/* + * 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::KTCavityNoiseGenerator(const string& name) : + KTGaussianNoiseGenerator(name), + fF0(25.904e9), + fQL(625.0), + fQ0(1.e4), + fA(0.90), + fTLineStart(80.0), + fTLineEnd(5.2), + fTCav(80.0), + fTIsol(5.2), + fEpsilon(0.5), + fFLo(25.9702e9), + fTransformFlag("ESTIMATE"), + fNoiseScaling(1.0) + { + } + + KTCavityNoiseGenerator::~KTCavityNoiseGenerator() + { + } + + bool KTCavityNoiseGenerator::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)); // Cavity noise should have fRNG() with default (mean, sigma), not inherited from KTGaussianNoiseGenerator + + if (node->has("cavity")) + { + const scarab::param_node& m = (*node)["cavity"].as_node(); + fF0 = m.get_value("f0", fF0); + fQL = m.get_value("Q_L", fQL); + fQ0 = m.get_value("Q0", fQ0); + fA = m.get_value("A", fA); + fTLineStart = m.get_value("T_line_start", fTLineStart); + fTLineEnd = m.get_value("T_line_end", fTLineEnd); + fTCav = m.get_value("T_cav", fTCav); + fTIsol = m.get_value("T_isol", fTIsol); + fEpsilon = m.get_value("epsilon", fEpsilon); + fFLo = m.get_value("f_lo", fFLo); + } + + 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; + + bool isComplex = dynamic_cast< KTTimeSeriesFFTW* >(data.GetTimeSeries(0)) != NULL; + + KTFrequencySpectrumFFTW spec(sliceSize, -fs*0.5, fs*0.5, false); + spec.SetNTimeBins(sliceSize); + + if (isComplex) + { + for (unsigned k = 0; k < sliceSize; ++k) + { + double fIf = (k <= N2) ? k * df : (static_cast(k) - static_cast(sliceSize)) * df; + double fRf = fIf + fFLo; // Down-converted + double pBin = NoisePSD(fRf) * df; // PSD -> power in one FFT bin + 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()); + } + } + else + { + for (unsigned k = 0; k <= N2; ++k) + { + double fIf = k * df; + double fRf = fIf + fFLo; // Down-converted + double pBin = NoisePSD(fRf) * df; // PSD -> power in one FFT bin + 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() + + 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; + } + + 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; + } + + double KTCavityNoiseGenerator::Eta(double x) const + { + return x/(std::exp(x)-1.0) + 0.5; + } + + double KTCavityNoiseGenerator::NoisePSD(double f) const + { + const double lor = 1.0 / ( 1.0 + ( 2.0 * fQL * (f - fF0) / fF0)*( 2.0 * fQL * (f - fF0) / fF0) ); // Lorentzian + + const double kB = KTMath::BoltzmannConstant(); + const double hbar= KTMath::ReducedPlanckConstant(); + const double omega = KTMath::TwoPi() * f; + + const double Tcav = fTCav * Eta(hbar*omega/(kB*fTCav)); + const double Tisol = fTIsol * Eta(hbar*omega/(kB*fTIsol)); + + const double P_cav = kB*Tcav * (4.*fG/((1.+fG)*(1.+fG))) * lor; + const double P_loss = kB*( fA*fTLineStart + (1.-fA)*fTLineEnd ); + const double P_isol = kB*Tisol * (1. - (4.*fG/((1.+fG)*(1.+fG)))*lor); + const double P_amp = hbar*omega / fEpsilon; + + 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..dbc61ce34 --- /dev/null +++ b/Source/Simulation/KTCavityNoiseGenerator.hh @@ -0,0 +1,124 @@ +/* + * KTCavityNoiseGenerator.hh + * + * Created on: May 28, 2025 + * Author: ehtkarim +*/ + +#ifndef KTCAVITYNOISEGENERATOR_HH_ +#define KTCAVITYNOISEGENERATOR_HH_ + +#include "KTGaussianNoiseGenerator.hh" + +#include + +namespace Katydid +{ + + /*! + @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: + KTCavityNoiseGenerator(const std::string& name = "cavity-noise-generator"); + virtual ~KTCavityNoiseGenerator(); + + virtual bool ConfigureDerivedGenerator(const scarab::param_node* node); + + protected: + double fF0; + double fA; + double fQL; + double fQ0; + double fTLineStart; + double fTLineEnd; + double fTCav; + double fTIsol; + double fEpsilon; + double fFLo; + + double GetQL() const; + double GetQ0() const; + + void SetQL(double ql); + void SetQ0(double q0); + + std::string fTransformFlag; + double fNoiseScaling; + + double NoisePSD(double f) const; + double Eta(double x) const; + + private: + double fG; + + public: + virtual bool GenerateTS(KTTimeSeriesData& data); + }; + + inline double KTCavityNoiseGenerator::GetQL() const + { + return fQL; + } + + inline void KTCavityNoiseGenerator::SetQL(double ql) + { + fQL = ql; + fG = fQ0 / fQL; + return; + } + + inline double KTCavityNoiseGenerator::GetQ0() const + { + return fQ0; + } + + inline void KTCavityNoiseGenerator::SetQ0(double q0) + { + fQ0 = q0; + fG = fQ0 / fQL; + return; + } + +} /* namespace Katydid */ + +#endif /* KTCAVITYNOISEGENERATOR_HH_ */ diff --git a/Source/Simulation/KTDCOffsetGenerator.cc b/Source/Simulation/KTDCOffsetGenerator.cc index d648451aa..b8d0f85fd 100644 --- a/Source/Simulation/KTDCOffsetGenerator.cc +++ b/Source/Simulation/KTDCOffsetGenerator.cc @@ -35,17 +35,17 @@ namespace Katydid { if (node == NULL) return false; - const KTParamArray* offsetPairs = node->ArrayAt("offsets"); + const scarab::param_array* offsetPairs = node-> array_at("offsets"); //change 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)); + 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); fOffsets[pair.first] = pair.second; } diff --git a/Source/Simulation/KTGaussianNoiseGenerator.cc b/Source/Simulation/KTGaussianNoiseGenerator.cc index 541c96892..bdc46e288 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" @@ -11,8 +13,10 @@ #include "KTMath.hh" #include "KTTimeSeriesData.hh" #include "KTTimeSeries.hh" +#include "KTTimeSeriesFFTW.hh" #include +#include using std::string; @@ -24,7 +28,10 @@ namespace Katydid KTGaussianNoiseGenerator::KTGaussianNoiseGenerator(const string& name) : KTTSGenerator(name), - fRNG() + fRNG(), + fSigmaPSD(0.0), + fGain(1.0), + fResistance(50.0) // (Ω) { } @@ -37,16 +44,63 @@ 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")) + { + sigma = std::sqrt(KTMath::BoltzmannConstant() * 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)); + 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(); @@ -61,12 +115,20 @@ namespace Katydid continue; } - //double binCenter = 0.5 * binWidth; - for (unsigned iBin = 0; iBin < sliceSize; iBin++) + double binCenter = 0.5 * binWidth; + if (auto* tsFFTW = dynamic_cast(timeSeries)) // Handling complex FFTW time series correctly + { + for (unsigned iBin = 0; iBin < sliceSize; ++iBin) + { + tsFFTW->SetRect(iBin, tsFFTW->GetReal(iBin) + fRNG() * kInvSqrt2, tsFFTW->GetImag(iBin) + fRNG() * kInvSqrt2); // Complex white-Gaussian noise + } + } + else { - timeSeries->SetValue(iBin, fRNG() + timeSeries->GetValue(iBin)); - //binCenter += binWidth; - //KTDEBUG(genlog, iBin << " " << timeSeries->GetValue(iBin)); + for (unsigned iBin = 0; iBin < sliceSize; ++iBin) + { + timeSeries->SetValue(iBin, timeSeries->GetValue(iBin) + fRNG()); + } } } diff --git a/Source/Simulation/KTGaussianNoiseGenerator.hh b/Source/Simulation/KTGaussianNoiseGenerator.hh index fde677024..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. @@ -61,8 +64,11 @@ namespace Katydid double GetSigma() const; void SetSigma(double sigma); - private: + protected: KTRNGGaussian<> fRNG; + double fSigmaPSD; + double fGain; + double fResistance; public: virtual bool GenerateTS(KTTimeSeriesData& data); diff --git a/Source/Simulation/KTMeasuredNoiseGenerator.cc b/Source/Simulation/KTMeasuredNoiseGenerator.cc new file mode 100644 index 000000000..c87bede5d --- /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) * sliceSize; + + if (isComplex) + { + for (unsigned k = 0; k < sliceSize; ++k) + { + const double fIf = (k <= N2) ? k * df : (static_cast(k) - static_cast(sliceSize)) * df; + const double fAbs = std::fabs(fIf); + + const double psd = DrawPSD(fAbs); + 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 fIf = k * df; + const double psd = DrawPSD(fIf); + 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(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_ */ + diff --git a/Source/Simulation/KTTSGenerator.cc b/Source/Simulation/KTTSGenerator.cc index ae0f92bef..26339db23 100644 --- a/Source/Simulation/KTTSGenerator.cc +++ b/Source/Simulation/KTTSGenerator.cc @@ -12,12 +12,13 @@ #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" :obsolete #include @@ -165,7 +166,7 @@ namespace Katydid } - newHeader->SetTimestamp(get_absolute_time_string()); + newHeader->SetTimestamp(scarab::get_absolute_time_string()); return newHeader; } @@ -232,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 ed2678d36..9b6e63a06 100644 --- a/Source/Simulation/KTTSGenerator.hh +++ b/Source/Simulation/KTTSGenerator.hh @@ -125,7 +125,7 @@ namespace Katydid private: Nymph::KTSignalOneArg< KTEggHeader* > fHeaderSignal; Nymph::KTSignalData fDataSignal; - Nymph::KTSignalOneArg< void > fDoneSignal; + Nymph::KTSignalOneArg< void > fDoneSignal; Nymph::KTSignalOneArg< const KTProcSummary* > fSummarySignal; }; diff --git a/Source/Utility/KTMath.hh b/Source/Utility/KTMath.hh index ba276bf65..4083329b1 100644 --- a/Source/Utility/KTMath.hh +++ b/Source/Utility/KTMath.hh @@ -21,15 +21,21 @@ namespace Katydid /* * From ROOT's TMath * */ /* ************************* */ - inline double Pi() { return 3.14159265358979323846; } - inline double TwoPi() { return 6.28318530717958623; } - inline double PiOver2() { return 1.57079632679489656; } - inline double PiOver4() { return 0.785398163397448279; } - inline double InvPi() { return 0.318309886183790691; } - inline double RadToDeg() { return 57.2957795130823229; } - inline double DegToRad() { return 1.74532925199432955e-02; } - inline double Sqrt2() { return 1.4142135623730950488016887242097; } - + constexpr double Pi() { return 3.14159265358979323846; } + constexpr double TwoPi() { return 6.28318530717958623; } + constexpr double PiOver2() { return 1.57079632679489656; } + constexpr double PiOver4() { return 0.785398163397448279; } + constexpr double InvPi() { return 0.318309886183790691; } + constexpr double RadToDeg() { return 57.2957795130823229; } + constexpr double DegToRad() { return 1.74532925199432955e-02; } + constexpr double Sqrt2() { return 1.4142135623730950488016887242097; } + + /* ************************ */ + /* * Physical constants * */ + /* ************************ */ + + constexpr double BoltzmannConstant() { return 1.380649e-23; } + constexpr double ReducedPlanckConstant() { return 1.054571817e-34; } /// Round to nearest integer. Rounds half integers to the nearest even integer. /// From ROOT's TMath diff --git a/Source/Utility/KTRandom.hh b/Source/Utility/KTRandom.hh index c82e2780d..a3e7689dc 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) {