diff --git a/src/Analysis.cpp b/src/Analysis.cpp index da622dd..2df27a4 100644 --- a/src/Analysis.cpp +++ b/src/Analysis.cpp @@ -100,6 +100,11 @@ void Analysis::run( std::vector& dataFileSets ){ // Merge sections mergeSections(dataFileSets); + if (ptrControl->doesPeformEOFBasedDenoising() && ptrControl->getTimingEOFBasedDenoising() == Control::BEFORE_DECIMATION) { + // Modify time-series data based on the EOF analysis + modifyTimeSeriesBasedOnEOFAnalysis(dataFileSets); + } + // Apply decimation decimation(dataFileSets); @@ -1090,7 +1095,12 @@ void Analysis::preprocessing( std::vector& dataFi } } } - + + if (ptrControl->doesPeformEOFBasedDenoising() && ptrControl->getTimingEOFBasedDenoising() == Control::AFTER_DEGITAL_FILTERS) { + // Modify time-series data based on the EOF analysis + modifyTimeSeriesBasedOnEOFAnalysis(dataFileSets); + } + if( (ptrControl->getParamsForPrewhitening()).applyPrewhitening ){ m_coefficientsOfARModel = new std::vector[numChannels]; const Control::ParamsForPrewhitening params = ptrControl->getParamsForPrewhitening(); @@ -2401,6 +2411,94 @@ void Analysis::calculateRotatedFields( const int numSegmentsTotal, std::complex< } +// Output average spectrum +void Analysis::modifyTimeSeriesBasedOnEOFAnalysis(std::vector& dataFileSets) { + + OutputFiles* ptrOutputFiles = OutputFiles::getInstance(); + ptrOutputFiles->writeLogMessage("Modify time-series data based on the EOF analysis"); + + const Control* const ptrControl = Control::getInstance(); + const double samplingFrequency = ptrControl->getSamplingFrequency(); + + int iSection(0); + int icount(0); + double matrix[3] = { 0.0, 0.0, 0.0 }; + for (std::vector::iterator itr = dataFileSets.begin(); itr != dataFileSets.end(); ++itr, ++iSection) { + const int numData = itr->numDataPoints; + for (int iChan = 0; iChan < 2; ++iChan) { + ptrOutputFiles->writeLogMessage("Secton " + Util::toString(iSection) + ", Channel " + Util::toString(iChan)); + // Make new data records by subtracting mean + const double mean = Util::calculateMeanValue(numData, itr->dataFile[iChan].data); + ptrOutputFiles->writeLogMessage("Subtract mean (" + Util::toString(mean) + ")"); + for (int i = 0; i < numData; ++i) { + itr->dataFile[iChan].data[i] -= mean; + } + } + for (int i = 0; i < numData; ++i) { + matrix[0] += pow(itr->dataFile[0].data[i], 2); + matrix[2] += pow(itr->dataFile[1].data[i], 2); + matrix[1] += itr->dataFile[0].data[i] * itr->dataFile[1].data[i]; + ++icount; + } +#ifdef _DEBUG_WRITE + const int numChannels = ptrControl->getNumberOfChannels(); + for (int iChan = 0; iChan < numChannels; ++iChan) { + std::ostringstream oss; + oss << "ts_sect_" << iSection << "_chan_" << iChan << "_before_EOF.txt"; + std::ofstream ofs; + ofs.open(oss.str().c_str(), std::ios::out); + if (ofs.fail()) { + ptrOutputFiles->writeLogMessage("File open error !! : " + oss.str()); + } + for (int i = 0; i < numData; ++i) { + const double elapsedTime = static_cast(i) / samplingFrequency; + ofs << std::setprecision(8) << std::scientific << itr->dataFile[iChan].data[i] << std::endl; + } + ofs.close(); + } +#endif + } + matrix[0] /= static_cast(icount); + matrix[1] /= static_cast(icount); + matrix[2] /= static_cast(icount); + // Calcualte all eigenvalues and eigenvectors of a real symmetric + double eigenVectors[4] = { 0.0, 0.0, 0.0, 0.0 }; + double eigenValues[2] = { 0.0, 0.0 }; + Util::calculateEigenValuesAndVectorsOfRealSymmetricMatrix(2, matrix, eigenValues, eigenVectors); + ptrOutputFiles->writeLogMessage("The first mode: eigen value = " + Util::toString(eigenValues[0]) + + ", eigen vector = (" + Util::toString(eigenVectors[0]) + ", " + Util::toString(eigenVectors[1]) + ")"); + ptrOutputFiles->writeLogMessage("The second mode: eigen value = " + Util::toString(eigenValues[1]) + + ", eigen vector = (" + Util::toString(eigenVectors[2]) + ", " + Util::toString(eigenVectors[3]) + ")"); + iSection = 0; + for (std::vector::iterator itr = dataFileSets.begin(); itr != dataFileSets.end(); ++itr, ++iSection) { + const int numData = itr->numDataPoints; + for (int i = 0; i < numData; ++i) { + const double val0 = itr->dataFile[0].data[i]; + const double val1 = itr->dataFile[1].data[i]; + itr->dataFile[0].data[i] = val0 * eigenVectors[0] + val1 * eigenVectors[1]; + itr->dataFile[1].data[i] = val0 * eigenVectors[2] + val1 * eigenVectors[3]; + } +#ifdef _DEBUG_WRITE + const int numChannels = ptrControl->getNumberOfChannels(); + for (int iChan = 0; iChan < numChannels; ++iChan) { + std::ostringstream oss; + oss << "ts_sect_" << iSection << "_chan_" << iChan << "_after_EOF.txt"; + std::ofstream ofs; + ofs.open(oss.str().c_str(), std::ios::out); + if (ofs.fail()) { + ptrOutputFiles->writeLogMessage("File open error !! : " + oss.str()); + } + for (int i = 0; i < numData; ++i) { + const double elapsedTime = static_cast(i) / samplingFrequency; + ofs << std::setprecision(8) << std::scientific << itr->dataFile[iChan].data[i] << std::endl; + } + ofs.close(); + } +#endif + } + +} + // Output average spectrum void Analysis::outputAverageSpectrum( std::complex** cdata, const int numData, const int section, const bool afterPreprocessing, const bool afterCalibration ) const{ diff --git a/src/Analysis.h b/src/Analysis.h index 7a56675..e8305b3 100644 --- a/src/Analysis.h +++ b/src/Analysis.h @@ -195,6 +195,9 @@ class Analysis{ const int numSegmentsTotal, std::complex** ftval, const std::vector< std::pair >& times, std::ofstream& ofsResp, std::ofstream& ofsRhoaPhs ) = 0; + // Output average spectrum + void modifyTimeSeriesBasedOnEOFAnalysis(std::vector& dataFileSets); + // Output average spectrum void outputAverageSpectrum( std::complex** cdata, const int numData, const int section, const bool afterPreprocessing, const bool afterCalibration ) const; diff --git a/src/CommonParameters.h b/src/CommonParameters.h index 8c9fa91..9793b1d 100644 --- a/src/CommonParameters.h +++ b/src/CommonParameters.h @@ -64,7 +64,7 @@ const static double EPS = 1.0e-20; static char programName[]="TRACMT"; -static char version[] = "v1.3.4"; +static char version[] = "v1.3.6"; } diff --git a/src/Control.cpp b/src/Control.cpp index a990d2e..5415849 100644 --- a/src/Control.cpp +++ b/src/Control.cpp @@ -65,6 +65,7 @@ Control::Control() : m_cutoffFrequencyForIIRLowPassFilter(-1.0), m_doesApplyIIRHighPassFilter(false), m_doesApplyIIRLowPassFilter(false), + m_doesPeformEOFBasedDenoising(false), m_elogMTReadingOption(Control::NOT_SPECIFIED), m_errorEstimationMethod(Control::FIXED_WEIGHTS_BOOTSTRAP), m_numOutputVariables(-1), @@ -92,6 +93,7 @@ Control::Control() : m_readAtsBinary(false), m_readElogDualBinary(false), m_readElogMTBinary(false), + m_timingEOFBasedDenoising(Control::BEFORE_DECIMATION), m_typeOfElogDual(Control::ELOG1K), m_typeOfElogMT(Control::ELOGMT_ADU_MODE) { @@ -292,6 +294,11 @@ double Control::getAzimuth( const int iChan ) const{ return m_azimuths[iChan]; } +// Get channel index +bool Control::doesPeformEOFBasedDenoising() const { + return m_doesPeformEOFBasedDenoising; +} + // Get channel index int Control::getChannelIndex( const int dataType, const int index ) const{ @@ -666,6 +673,13 @@ int Control::getProcedureType () const{ } +// Get type of ELOG-Dual +int Control::getTimingEOFBasedDenoising() const { + + return m_timingEOFBasedDenoising; + +} + // Get type of ELOG-Dual int Control::getTypeOfElogDual() const { return m_typeOfElogDual; @@ -1303,6 +1317,12 @@ void Control::readParameterFile(){ m_rangeOfSectionsForMerge.push_back ( std::make_pair(startSectionIndex, endSectionIndex ) ); } } + else if (line.find("EOF_ANALYSIS") != std::string::npos) { + m_doesPeformEOFBasedDenoising = true; + int ibuf(0); + ifs >> ibuf; + m_timingEOFBasedDenoising = ibuf; + } else if( line.find("END") != std::string::npos ){ break; } @@ -1820,6 +1840,20 @@ void Control::readParameterFile(){ } ptrOutputFiles->writeLogMessage(oss.str(), false); } + if (doesPeformEOFBasedDenoising()) { + switch (getTimingEOFBasedDenoising()) + { + case Control::BEFORE_DECIMATION: + ptrOutputFiles->writeLogMessage("EOF analysis is performed before decimation", false); + break; + case Control::AFTER_DEGITAL_FILTERS: + ptrOutputFiles->writeLogMessage("EOF analysis is performed after the application of digital filters", false); + break; + default: + ptrOutputFiles->writeErrorMessage("Unsupported timing of EOF analysis: " + Util::toString(getTimingEOFBasedDenoising())); + break; + } + } if( m_paramsForDecimation.applyDecimation ){ ptrOutputFiles->writeLogMessage("Parameters for decimation : ", false); ptrOutputFiles->writeLogMessage(" Decimation interval: " + Util::toString(m_paramsForDecimation.decimationInterval), false); diff --git a/src/Control.h b/src/Control.h index 992098d..3d10d78 100644 --- a/src/Control.h +++ b/src/Control.h @@ -106,6 +106,12 @@ class Control{ NUM_OF_THRESHOLD_TYPE_OF_DIFFERENCE_OF_RESPONSE_FUCTIONS }; + // Segment length and the index of the frequency where response functions are estimated + enum TimingEOFBasedDenoising{ + BEFORE_DECIMATION = 0 , + AFTER_DEGITAL_FILTERS, + }; + // Segment length and the index of the frequency where response functions are estimated struct SegmentInfo{ int length; @@ -291,6 +297,9 @@ class Control{ // Get flag specifing whether input file is ELOG-MT binary file bool doesReadElogMTBinary() const; + // Get azimuth + bool doesPeformEOFBasedDenoising() const; + // Get azimuth double getAzimuth( const int iChan ) const; @@ -463,7 +472,10 @@ class Control{ int getNumStartTimesSections() const; // Get procedure type - int getProcedureType () const; + int getProcedureType() const; + + // Get timing of denoising based on EOF + int getTimingEOFBasedDenoising() const; // Get type of ELOG-Dual int getTypeOfElogDual() const; @@ -512,6 +524,9 @@ class Control{ // Flag specifing whether IIR low-pass filter is applied bool m_doesApplyIIRLowPassFilter; + // Option of ELOG-Dual binary data reading + bool m_doesPeformEOFBasedDenoising; + // Option of ELOG-Dual binary data reading int m_elogDualReadingOption; @@ -671,6 +686,9 @@ class Control{ // Flag specifing whether ELOG-MT binary is read bool m_readElogMTBinary; + // Type of ELOG-Dual + int m_timingEOFBasedDenoising; + // Type of ELOG-Dual int m_typeOfElogDual; diff --git a/src/ElogDual.cpp b/src/ElogDual.cpp index 8994582..f30c864 100644 --- a/src/ElogDual.cpp +++ b/src/ElogDual.cpp @@ -130,7 +130,7 @@ void ElogDual::readElogBinaryFilesUnderADirectory(const std::string& directory, int counter(0); for (auto& p : dirItr) { const std::string fileName = p.path().string(); - if (fileName.find(stringCompared) != std::string::npos) { + if (fileName.find(stringCompared) != std::string::npos && Util::extractExtensionOfFileName(fileName).find("dat") != std::string::npos) { readElogBinaryFile(fileName, numSkipData, numDataPoints, counter, ex, ey); if (counter >= numSkipData + numDataPoints){ break; diff --git a/src/LapackInterface.cpp b/src/LapackInterface.cpp index 91c5853..f480070 100644 --- a/src/LapackInterface.cpp +++ b/src/LapackInterface.cpp @@ -42,19 +42,22 @@ void LapackInterface::calculateEigenValuesAndVectorsOfRealSymmetricMatrix( const integer n = static_cast(dimension); integer lda = n; - integer lwork = n * n; + integer dum1 = 1; + integer dum2 = -1; + integer nb = ilaenv_(&dum1, "DSYTRD", "L", &n, &dum2, &dum2, &dum2); + integer lwork = (nb + 2) * n; double* work = new double[lwork]; integer info = 0; - dsyev_("V", "L", &n, matrix, &lda, vectors, work, &lwork, &info); - delete [] work; + delete[] work; - if( info < 0 ){ + if (info < 0) { OutputFiles* ptrOutputFiles = OutputFiles::getInstance(); - ptrOutputFiles->writeErrorMessage("An argument had an illegal value : info="+Util::toString(info) ); - }else if( info > 0 ){ + ptrOutputFiles->writeErrorMessage("An argument had an illegal value : info=" + Util::toString(info)); + } + else if (info > 0) { OutputFiles* ptrOutputFiles = OutputFiles::getInstance(); - ptrOutputFiles->writeErrorMessage("Eigenvalue calculation is not converged : info="+Util::toString(info) ); + ptrOutputFiles->writeErrorMessage("Eigenvalue calculation is not converged : info=" + Util::toString(info)); } } diff --git a/src/Util.cpp b/src/Util.cpp index ab48a15..0c4f852 100644 --- a/src/Util.cpp +++ b/src/Util.cpp @@ -467,15 +467,28 @@ double Util::calculateDeterminantOfMatrix( const int dimension, const double* co } // Calcualte all eigenvalues and eigenvectors of a real symmetric -void Util::calculateEigenValuesAndVectorsOfRealSymmetricMatrix( const int dimension, const double* const matrix, - double* eigenValues, double* eigenVectors ){ +void Util::calculateEigenValuesAndVectorsOfRealSymmetricMatrix(const int dimension, const double* const matrix, + double* eigenValues, double* eigenVectors) { OutputFiles* ptrOutputFiles = OutputFiles::getInstance(); - if( dimension < 1 ){ - ptrOutputFiles->writeErrorMessage("Dimension of linear equation is less than 1" ); + if (dimension < 1) { + ptrOutputFiles->writeErrorMessage("Dimension of linear equation is less than 1"); + } + + int icount(0); + // Column major + for (int col = 0; col < dimension; ++col) { + for (int row = 0; row < col; ++row) { + const int index = col * dimension + row; + eigenVectors[index] = 0.0; + } + for (int row = col; row < dimension; ++row) { + const int index = col * dimension + row; + eigenVectors[index] = matrix[icount]; + ++icount; + } } - memcpy(eigenVectors, matrix, sizeof(double)*dimension*dimension); LapackInterface::calculateEigenValuesAndVectorsOfRealSymmetricMatrix(dimension, eigenVectors, eigenValues); }