modify a function to read ELOG-DUAL binary files

src/Analysis.cpp

@@ -100,6 +100,11 @@ void Analysis::run( std::vector<CommonParameters::DataFileSet>& 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<CommonParameters::DataFileSet>& 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<double>[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<CommonParameters::DataFileSet>& 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<CommonParameters::DataFileSet>::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<double>(i) / samplingFrequency;
+				ofs << std::setprecision(8) << std::scientific << itr->dataFile[iChan].data[i] << std::endl;
+			}
+			ofs.close();
+		}
+#endif
+	}
+	matrix[0] /= static_cast<double>(icount);
+	matrix[1] /= static_cast<double>(icount);
+	matrix[2] /= static_cast<double>(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<CommonParameters::DataFileSet>::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<double>(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<double>** cdata, const int numData, const int section, 
 	const bool afterPreprocessing, const bool afterCalibration ) const{

src/Analysis.h

@@ -195,6 +195,9 @@ class Analysis{
 		const int numSegmentsTotal, std::complex<double>** ftval, const std::vector< std::pair<std::string, std::string> >& times,
 		std::ofstream& ofsResp, std::ofstream& ofsRhoaPhs ) = 0;
 
+	// Output average spectrum
+	void modifyTimeSeriesBasedOnEOFAnalysis(std::vector<CommonParameters::DataFileSet>& dataFileSets);
+
 	// Output average spectrum
 	void outputAverageSpectrum( std::complex<double>** cdata, const int numData, const int section, 
 		const bool afterPreprocessing, const bool afterCalibration ) const;

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";
 
 }
 

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);

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;
 

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;

src/LapackInterface.cpp

@@ -42,19 +42,22 @@ void LapackInterface::calculateEigenValuesAndVectorsOfRealSymmetricMatrix( const
 
 	integer n = static_cast<integer>(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));
 	}
 
 }

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);
 
 }