From 8348ca363b119570f1fa7104902b98e6c1d5c212 Mon Sep 17 00:00:00 2001 From: Penny Slocum Date: Thu, 9 May 2024 16:15:25 -0400 Subject: [PATCH] Progress merging single-mode without indices, and multimode. --- Source/Core/LMCHFSSResponseFileHandler.cc | 294 ++++++++++++++-------- Source/Core/LMCHFSSResponseFileHandler.hh | 57 ++--- 2 files changed, 224 insertions(+), 127 deletions(-) diff --git a/Source/Core/LMCHFSSResponseFileHandler.cc b/Source/Core/LMCHFSSResponseFileHandler.cc index 205e5650..399bb29e 100644 --- a/Source/Core/LMCHFSSResponseFileHandler.cc +++ b/Source/Core/LMCHFSSResponseFileHandler.cc @@ -19,16 +19,16 @@ namespace locust LOGGER( lmclog, "HFSSResponseFileHandlerCore" ); HFSSResponseFileHandlerCore::HFSSResponseFileHandlerCore(): - //fTFNBins(1000), - //fFIRNBins(2000), + fTFNBins(1000), + fFIRNBins(2000), fCropIndex(2100), - //fResolution(1e-12), + fResolution(1e-12), fCharacteristicImpedance(1.0), fNSkips(1), fNModes(2), fComplexFFT(), fHFSSFiletype(""), - //fIsFIRCreated(false), + fIsFIRCreated(false), fWindowName("tukey"), fWindowParam(0.5), fPrintFIR ( false ), @@ -47,39 +47,39 @@ namespace locust { fPrintFIR = aParam["print-fir-debug"]().as_bool(); } - if( aParam.has( "n-modes" ) ) - { - fNModes = aParam["n-modes"]().as_int(); - } + if( aParam.has( "n-modes" ) ) + { + fNModes = aParam["n-modes"]().as_int(); + } - fFilterComplexArray.resize(2); - fTFNBinsArray.resize(2); + fFilterComplexArray.resize(2); // TE or TM + fTFNBinsArray.resize(2); fFIRNBinsArray.resize(2); fResolutionArray.resize(2); - fIsFIRCreatedArray.resize(2); - for ( unsigned bTE=0; bTE<2; bTE++) - { - fFilterComplexArray[bTE].resize(fNModes); - fTFNBinsArray[bTE].resize(fNModes); - fFIRNBinsArray[bTE].resize(fNModes); - fResolutionArray[bTE].resize(fNModes); - fIsFIRCreatedArray[bTE].resize(fNModes); - for (unsigned l=0; l inputBuffer) + { + double convolution=0.0; + if(fFIRNBins<=0) + { + LERROR(lmclog,"Number of bins in the filter should be positive"); + } + int firBinNumber=0; + for (auto it = inputBuffer.begin();it!=inputBuffer.end(); ++it) + { + convolution+=*(it)*fFilter[firBinNumber]; + firBinNumber++; + } + + return convolution; + } + + double HFSSResponseFileHandlerCore::ConvolveWithFIRFilterArray(int bTE, int l, int m, int n, std::deque inputBuffer) { double convolution=0.0; @@ -134,6 +152,43 @@ namespace locust return convolution; } + std::pair HFSSResponseFileHandlerCore::ConvolveWithComplexFIRFilter(std::deque inputBuffer) + { + double convolutionMag = 0.0; + double convolutionValueReal = 0.0; + double convolutionValueImag = 0.0; + + if(fFIRNBins<=0) + { + LERROR(lmclog,"Number of bins in the filter should be positive"); + } + int firBinNumber=0; + + for (auto it = inputBuffer.begin();it!=inputBuffer.end(); ++it) + { + convolutionValueReal += *(it)*fFilterComplex[firBinNumber][0]; + convolutionValueImag += *(it)*fFilterComplex[firBinNumber][1]; + firBinNumber++; + } + std::pair complexConvolution; + double complexPhase = atan(convolutionValueImag/convolutionValueReal); + double complexMag = pow(convolutionValueReal*convolutionValueReal + convolutionValueImag*convolutionValueImag, 0.5); + + /* Note that the convolution above uses a real signal to drive both the real and + * imaginary filter components, without a phase shift. For sinusoidal signals, this + * seems like a reasonable approach. It also increases the effective gain of + * the complex filter, compared to a real filter. To align the complex gain with + * that of a real filter, a gain factor of 0.5x for sinusoidal drive signals is + * applied below. For non-sinusoidal drive signals, this scaling should not necessarily + * be applicable. + */ + + double gainFactor = 0.5; + + return std::make_pair(gainFactor*complexMag, complexPhase); + } + + std::pair HFSSResponseFileHandlerCore::ConvolveWithComplexFIRFilterArray(int bTE, int l, int m, int n, std::deque inputBuffer) { double convolutionMag = 0.0; @@ -186,17 +241,18 @@ namespace locust } - bool TFFileHandlerCore::ConvertTFtoFIR(int bTE, int l, int m, int n, std::vector> &tfArray, bool GeneratedTF) + + bool TFFileHandlerCore::ConvertTFtoFIR(std::vector> &tfArray, bool GeneratedTF) { - if(fTFNBinsArray[bTE][l][m][n]<=0) + if(fTFNBins<=0) { LERROR(lmclog,"The size of transfer function has to be positive integer"); return false; } //Might need to be moved to a different function - fTFComplex=(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fTFNBinsArray[bTE][l][m][n]); + fTFComplex=(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fTFNBins); - for (int i = 0; i < fTFNBinsArray[bTE][l][m][n]; ++i) + for (int i = 0; i < fTFNBins; ++i) { fTFComplex[i][0]=tfArray.at(i).real(); fTFComplex[i][1]=tfArray.at(i).imag(); @@ -204,48 +260,89 @@ namespace locust // this length is somewhat arbitrary, but works for now // future improvement: make it adjust depending on length of FIR - fFIRNBinsArray[bTE][l][m][n]=fTFNBinsArray[bTE][l][m][n]+ 2*fComplexFFT.GetShiftNBins(); + fFIRNBins = fTFNBins + 2*fComplexFFT.GetShiftNBins(); + if(GeneratedTF) + { + // If TF generated based on config file (frequency ranges from 0.9 - 1.1 times the center given in .json config file), calculate the TF bin width given number of bins (also set in .json config file)). + double AnalyticTFBinWidth = 2./9.*fInitialTFIndex/(1.0*fTFNBins); + fComplexFFT.SetupIFFTWindow(fTFNBins,fInitialTFIndex,AnalyticTFBinWidth, fWindowName, fWindowParam);//Uses binwidth as calculated in the previous line based on internally generated TF + } + else + { + // If TF read from externally generated TF file, use TF bin width as given in .json config file + fComplexFFT.SetupIFFTWindow(fTFNBins,fInitialTFIndex,fTFBinWidth, fWindowName, fWindowParam); + } + fFIRComplex = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fFIRNBins); + fComplexFFT.ApplyWindowFunction(fTFNBins, fTFComplex); + fComplexFFT.GenerateFIR(fTFNBins,fTFComplex,fFIRComplex); + fResolution = fComplexFFT.GetTimeResolution(); - if(GeneratedTF) - { - // If TF generated based on config file (frequency ranges from 0.9 - 1.1 times the center given in .json config file), calculate the TF bin width given number of bins (also set in .json config file)). - double AnalyticTFBinWidth = 2./9.*fInitialTFIndex/(1.0*fTFNBinsArray[bTE][l][m][n]); - fComplexFFT.SetupIFFTWindow(fTFNBinsArray[bTE][l][m][n],fInitialTFIndex,AnalyticTFBinWidth, fWindowName, fWindowParam);//Uses binwidth as calculated in the previous line based on internally generated TF - } - else + for (int i = 0; i < fFIRNBins; i++) + { + fFilter.push_back(fFIRComplex[i][0]); + } + fFilterComplex = fFIRComplex; + + LDEBUG( lmclog, "Finished IFFT to convert transfer function to FIR"); + return true; + } + + + bool TFFileHandlerCore::ConvertTFtoFIR(int bTE, int l, int m, int n, std::vector> &tfArray, bool GeneratedTF) { - // If TF read from externally generated TF file, use TF bin width as given in .json config file - fComplexFFT.SetupIFFTWindow(fTFNBinsArray[bTE][l][m][n],fInitialTFIndex,fTFBinWidth, fWindowName, fWindowParam); - } + if(fTFNBinsArray[bTE][l][m][n]<=0) + { + LERROR(lmclog,"The size of transfer function has to be positive integer"); + return false; + } + //Might need to be moved to a different function + fTFComplex=(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fTFNBinsArray[bTE][l][m][n]); - fFIRComplex=(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fFIRNBinsArray[bTE][l][m][n]); + for (int i = 0; i < fTFNBinsArray[bTE][l][m][n]; ++i) + { + fTFComplex[i][0]=tfArray.at(i).real(); + fTFComplex[i][1]=tfArray.at(i).imag(); + } + // this length is somewhat arbitrary, but works for now + // future improvement: make it adjust depending on length of FIR - fComplexFFT.ApplyWindowFunction(fTFNBinsArray[bTE][l][m][n], fTFComplex); + fFIRNBinsArray[bTE][l][m][n]=fTFNBinsArray[bTE][l][m][n]+ 2*fComplexFFT.GetShiftNBins(); + if(GeneratedTF) + { + // If TF generated based on config file (frequency ranges from 0.9 - 1.1 times the center given in .json config file), calculate the TF bin width given number of bins (also set in .json config file)). + double AnalyticTFBinWidth = 2./9.*fInitialTFIndex/(1.0*fTFNBinsArray[bTE][l][m][n]); + fComplexFFT.SetupIFFTWindow(fTFNBinsArray[bTE][l][m][n],fInitialTFIndex,AnalyticTFBinWidth, fWindowName, fWindowParam);//Uses binwidth as calculated in the previous line based on internally generated TF + } + else + { + // If TF read from externally generated TF file, use TF bin width as given in .json config file + fComplexFFT.SetupIFFTWindow(fTFNBinsArray[bTE][l][m][n],fInitialTFIndex,fTFBinWidth, fWindowName, fWindowParam); + } - fComplexFFT.GenerateFIR(fTFNBinsArray[bTE][l][m][n],fTFComplex,fFIRComplex); + fFIRComplex=(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fFIRNBinsArray[bTE][l][m][n]); + fComplexFFT.ApplyWindowFunction(fTFNBinsArray[bTE][l][m][n], fTFComplex); + fComplexFFT.GenerateFIR(fTFNBinsArray[bTE][l][m][n],fTFComplex,fFIRComplex); + fResolutionArray[bTE][l][m][n]=fComplexFFT.GetTimeResolution(); - fResolutionArray[bTE][l][m][n]=fComplexFFT.GetTimeResolution(); - for (int i = 0; i < fFIRNBinsArray[bTE][l][m][n]; i++) - { - fFilter.push_back(fFIRComplex[i][0]); - } - //fFilterComplex = fFIRComplex; + for (int i = 0; i < fFIRNBinsArray[bTE][l][m][n]; i++) + { + fFilter.push_back(fFIRComplex[i][0]); + } + //fFilterComplex = fFIRComplex; - LDEBUG( lmclog, "Finished IFFT to convert transfer function to FIR"); - return true; + LDEBUG( lmclog, "Finished IFFT to convert transfer function to FIR"); + return true; } - bool TFFileHandlerCore::ReadHFSSFile(int bTE, int l, int m, int n) //read HFSS file for a given mode + bool TFFileHandlerCore::ReadHFSSFile() //read HFSS file for a given mode { - - - if(fIsFIRCreatedArray[bTE][l][m][n]) + if(fIsFIRCreated) { return true; } - fTFNBinsArray[bTE][l][m][n]=0; + fTFNBins = 0; if(!ends_with(fHFSSFilename,".txt")) { LERROR(lmclog,"The TF file " << fHFSSFilename.c_str() <<"doesn't end in .txt"); @@ -256,11 +353,11 @@ namespace locust double tfImaginaryValue; std::vector> tfArray; std::fstream tfFile(fHFSSFilename.c_str(),std::ios::in); - if (tfFile.fail()) - { + if (tfFile.fail()) + { LERROR(lmclog,"The TF file " << fHFSSFilename.c_str() <<" doesn't exist"); return false; - } + } //logic copied from /LMCPatchSignalGenerator.cc int totalcount=0; @@ -290,10 +387,10 @@ namespace locust ++wordCount; } // The TF values from HFSS are in GHz, so need to convert to Hz - if(fTFNBinsArray[bTE][l][m][n]==0)fInitialTFIndex=tfIndex*pow(10.0,9); + if(fTFNBins == 0)fInitialTFIndex=tfIndex*pow(10.0,9); const std::complex temp(tfRealValue,tfImaginaryValue); tfArray.push_back(temp); - ++fTFNBinsArray[bTE][l][m][n]; + ++fTFNBins; } } } @@ -305,22 +402,22 @@ namespace locust return false; } - if (!ConvertTFtoFIR(bTE, l, m, n, tfArray, false ) ) //bool determines if TF was generated dynamically + if (!ConvertTFtoFIR(tfArray, false ) ) //bool determines if TF was generated dynamically { return false; } - if (!CropFIR(bTE, l, m, n, fFIRComplex, fConvertStoZ ) ) + if (!CropFIR(fFIRComplex, fConvertStoZ ) ) { return false; } - fIsFIRCreatedArray[bTE][l][m][n]=true; + fIsFIRCreated = true; if (fPrintFIR) { - PrintFIR( fFIRComplex, fFIRNBinsArray[bTE][l][m][n], fOutputPath+"/FIRhisto.root"); - PrintFIR( fTFComplex, fTFNBinsArray[bTE][l][m][n], fOutputPath+"/TFhisto.root"); + PrintFIR( fFIRComplex, fFIRNBins, fOutputPath+"/FIRhisto.root"); + PrintFIR( fTFComplex, fTFNBins, fOutputPath+"/TFhisto.root"); LPROG( lmclog, "Finished writing histos to output/FIRhisto.root and output/TFhisto.root"); LPROG( lmclog, "Press Return to continue, or Cntrl-C to quit."); @@ -330,7 +427,7 @@ namespace locust return true; } - bool TFFileHandlerCore::CropFIR(int bTE, int l, int m, int n, fftw_complex* anArray, bool bConvert) + bool TFFileHandlerCore::CropFIR(fftw_complex* anArray, bool bConvert) { if (bConvert) { @@ -345,7 +442,7 @@ namespace locust aSubArray[i][1] = anArray[i][1]; } anArray = aSubArray; - fFIRNBinsArray[bTE][l][m][n] = cropRange; + fFIRNBins = cropRange; } return true; } @@ -412,16 +509,14 @@ namespace locust } if (fPrintFIR) { - - scarab::path dataDir = TOSTRING(PB_DATA_INSTALL_DIR); - std::string modeIndexStr = std::to_string(bTE) + std::to_string(l) + std::to_string(m) + std::to_string(n); - std::string fileName = (dataDir / "../output/FIRhisto").string() + modeIndexStr + ".root"; + std::string modeIndexStr = std::to_string(bTE) + std::to_string(l) + std::to_string(m) + std::to_string(n); + std::string fileName = fOutputPath + modeIndexStr + ".root"; PrintFIR( fFilterComplexArray[bTE][l][m][n], fFIRNBinsArray[bTE][l][m][n], fileName ); LPROG( lmclog, "Finished writing histos to output/FIRhisto"+modeIndexStr+".root"); LPROG( lmclog, "Press Return to continue, or Cntrl-C to quit."); getchar(); } - fIsFIRCreatedArray[bTE][l][m][n]=true; + fIsFIRCreatedArray[bTE][l][m][n]=true; LDEBUG( lmclog, "Finished populating FIR filter with Green's function."); return true; @@ -431,7 +526,7 @@ namespace locust bool TFFileHandlerCore::ConvertAnalyticTFtoFIR(int bTE, int l, int m, int n, double initialFreq, std::vector> tfArray) { - //Replaces ReadHFSSFile() in the case where a Transfer Funciton has been generated analytically + //Replaces ReadHFSSFile() in the case where a Transfer Funciton has been generated analytically if(fIsFIRCreatedArray[bTE][l][m][n]) { @@ -441,7 +536,8 @@ namespace locust fTFNBinsArray[bTE][l][m][n]=tfArray.size(); fInitialTFIndex = initialFreq; - if(!ConvertTFtoFIR(bTE, l, m, n, tfArray, true)){ //bool determines if TF was generated dynamically + if(!ConvertTFtoFIR(bTE, l, m, n, tfArray, true)) + { //bool determines if TF was generated dynamically return false; } fIsFIRCreatedArray[bTE][l][m][n]=true; @@ -535,13 +631,13 @@ namespace locust return true; } - bool FIRFileHandlerCore::ReadHFSSFile(int bTE, int l, int m, int n) + bool FIRFileHandlerCore::ReadHFSSFile() { - if(fIsFIRCreatedArray[bTE][l][m][n]) - { - return true; - } - fFIRNBinsArray[bTE][l][m][n]=0; + if(fIsFIRCreated) + { + return true; + } + fFIRNBins = 0; if(!ends_with(fHFSSFilename,".txt")) { LERROR(lmclog,"The FIR file " << fHFSSFilename.c_str() <<" doesn't end in .txt"); @@ -551,27 +647,27 @@ namespace locust double filterMagnitude; FILE *firFile; firFile=fopen(fHFSSFilename.c_str(),"r"); - - if(!access(fHFSSFilename.c_str(), F_OK )) - { + if(!access(fHFSSFilename.c_str(), F_OK )) + { LERROR(lmclog,"The FIR file " << fHFSSFilename.c_str() <<" doesn't exist"); return false; - } + } //logic copied from /LMCPatchSignalGenerator.cc int count=0; - while (!feof(firFile)){ + while (!feof(firFile)) + { fscanf(firFile,"%lf %lf",&firIndex,&filterMagnitude); if (count%fNSkips==0) { fFilter.push_back(filterMagnitude); - ++fFIRNBinsArray[bTE][l][m][n]; + ++fFIRNBins; } ++count; } fclose(firFile); LDEBUG( lmclog, "Finished reading FIR file"); - fIsFIRCreatedArray[bTE][l][m][n]=true; + fIsFIRCreated = true; return true; } diff --git a/Source/Core/LMCHFSSResponseFileHandler.hh b/Source/Core/LMCHFSSResponseFileHandler.hh index c887a0dd..42691b68 100644 --- a/Source/Core/LMCHFSSResponseFileHandler.hh +++ b/Source/Core/LMCHFSSResponseFileHandler.hh @@ -30,16 +30,16 @@ namespace locust // Member functions virtual bool Configure( const scarab::param_node& aNode); - virtual bool ReadHFSSFile(int bTE, int l, int m, int n); - virtual double ConvolveWithFIRFilterArray(int bTE, int l, int m, int n, std::deque);// Convolve input signal (voltage or field) with FIR - //virtual double ConvolveWithFIRFilter(std::deque);// Convolve input signal (voltage or field) with FIR + virtual bool ReadHFSSFile(); + virtual double ConvolveWithFIRFilterArray(int bTE, int l, int m, int n, std::deque);// Convolve multimode input signal with FIR + virtual double ConvolveWithFIRFilter(std::deque);// Convolve input signal with FIR virtual std::pair ConvolveWithComplexFIRFilterArray(int bTE, int l, int m, int n, std::deque inputBuffer); - //virtual std::pair ConvolveWithComplexFIRFilter(std::deque inputBuffer); - //int GetFilterSize() const;//Number of entries in the filter - int GetFilterSizeArray(int bTE, int l, int m, int n) const;//Number of entries in the filter - int GetNModes() const; - //double GetFilterResolution() const;//Get the resolution of the filter - double GetFilterResolutionArray(int bTE, int l, int m, int n) const;//Get the resolution of the filter + virtual std::pair ConvolveWithComplexFIRFilter(std::deque inputBuffer); + int GetFilterSizeArray(int bTE, int l, int m, int n) const;//Number of entries in the filter + int GetFilterSize() const;//Number of entries in the filter + int GetNModes() const; + double GetFilterResolutionArray(int bTE, int l, int m, int n) const;//Get the resolution of the filter + double GetFilterResolution() const;//Get the resolution of the filter void PrintFIR( std::vector, int nBins, std::string filename ); void PrintFIR( fftw_complex* aFilter, int nBins, std::string filename ); bool WriteRootHisto( std::vector aFilter, int nBins, bool bIQ ); @@ -47,24 +47,24 @@ namespace locust protected: // Member variables + int fNModes; std::string fHFSSFilename; std::vector fFilter; - //fftw_complex* fFilterComplex; //Updating to use multimode compatible array - std::vector< std::vector< std::vector < std::vector< fftw_complex*>>>> fFilterComplexArray; - //int fTFNBins; - std::vector < std::vector < std::vector < std::vector < int >>>> fTFNBinsArray; - //int fFIRNBins; - std::vector < std::vector < std::vector < std::vector < int >>>> fFIRNBinsArray; - int fNModes; - //double fResolution; - std::vector < std::vector < std::vector < std::vector < double >>>> fResolutionArray; + std::vector< std::vector< std::vector < std::vector< fftw_complex*>>>> fFilterComplexArray; + fftw_complex* fFilterComplex; + std::vector < std::vector < std::vector < std::vector < int >>>> fTFNBinsArray; + int fTFNBins; + std::vector < std::vector < std::vector < std::vector < int >>>> fFIRNBinsArray; + int fFIRNBins; + std::vector < std::vector < std::vector < std::vector < double >>>> fResolutionArray; + double fResolution; + std::vector < std::vector < std::vector < std::vector < bool >>>> fIsFIRCreatedArray; + bool fIsFIRCreated; int fCropIndex; double fCharacteristicImpedance; int fNSkips; bool fHFSSFiletype; ComplexFFT fComplexFFT; - //bool fIsFIRCreated; - std::vector < std::vector < std::vector < std::vector < bool >>>> fIsFIRCreatedArray; std::string fWindowName; double fWindowParam; bool fPrintFIR; @@ -80,21 +80,21 @@ namespace locust //Member functions bool ends_with(const std::string &, const std::string &); }; - /* + inline int HFSSResponseFileHandlerCore::GetFilterSize() const { return fFIRNBins; - }*/ + } inline int HFSSResponseFileHandlerCore::GetFilterSizeArray(int bTE, int l, int m, int n) const { return fFIRNBinsArray[bTE][l][m][n]; } - /* + inline double HFSSResponseFileHandlerCore::GetFilterResolution() const { return fResolution; - }*/ + } inline int HFSSResponseFileHandlerCore::GetNModes() const { @@ -121,9 +121,9 @@ namespace locust // Member functions virtual bool Configure( const scarab::param_node& aNode) override; - bool ReadHFSSFile(int bTE, int l, int m, int n) override; + bool ReadHFSSFile(); //bool ConvertAnalyticTFtoFIR(double initialFreq, std::vector> tfArray); - bool ConvertAnalyticTFtoFIR(int bTE, int l, int m, int n, double initialFreq, std::vector> tfArray); + bool ConvertAnalyticTFtoFIR(int bTE, int l, int m, int n, double initialFreq, std::vector> tfArray); bool ConvertAnalyticGFtoFIR(int bTE, int l, int m, int n, std::vector > > gfArray); @@ -134,8 +134,9 @@ namespace locust // Member functions bool ConvertTFtoFIR(int bTE, int l, int m, int n, std::vector> &, bool GeneratedTF); + bool ConvertTFtoFIR(std::vector> &, bool GeneratedTF); bool ConvertStoZ(std::vector> &tfArray, bool bConvert); - bool CropFIR(int bTE, int l, int m, int n, fftw_complex* anArray, bool bConvert); + bool CropFIR(fftw_complex* anArray, bool bConvert); protected: //Member variables @@ -158,7 +159,7 @@ namespace locust // Member functions virtual bool Configure( const scarab::param_node& aNode) override; - bool ReadHFSSFile(int bTE, int l, int m, int n) override; + bool ReadHFSSFile(); }; } /*namespace locust*/