modify functoins to read filters from MTH5 files

Author: yoshiya-usui

docs/Manual_Of_TRACMT.pdf

@@ -111,6 +111,12 @@ void Analysis::run( std::vector<CommonParameters::DataFileSet>& dataFileSets ){
 	// Apply decimation
 	decimation(dataFileSets);
 
+#ifdef _MTH5
+	if (ptrControl->doesReadMTH5() && ptrControl->doesReadMTH5Filters()) {
+		(MTH5::getInstance())->readFiltersAll(ptrControl->getNumberOfChannels(), dataFileSets);
+	}
+#endif
+
 	// Read calibration files for main analysis
 	readCalibrationFiles(freqAll);
 
@@ -923,22 +929,6 @@ void Analysis::readCalibrationFiles( const std::vector<double>& freq ){
 			}
 		}
 	}
-#ifdef _MTH5
-	else if (ptrControl->doesReadMTH5()) {
-		const int numFilterInfo = ptrControl->getNumFilterInfoMTH5();
-		if (numFilterInfo > 0) {
-			if (numFilterInfo != ptrControl->getNumberOfChannels()) {
-				(OutputFiles::getInstance())->writeErrorMessage("Number of the filter info should be equal to channel number");
-			}
-			const MTH5* ptrMTH5 = MTH5::getInstance();
-			for (int iChan = 0; iChan < numFilterInfo; ++iChan) {
-				const std::string fileName = ptrControl->getFileNameForFilterInMTH5(iChan);
-				const std::string path = ptrControl->getPathFilterInMTH5(iChan);
-				ptrMTH5->makeCalibrationFile(iChan, freq);
-			}
-		}
-	}
-#endif
 
 	if (m_calibrationFunctions != NULL) {
 		delete[] m_calibrationFunctions;
@@ -1178,19 +1168,9 @@ void Analysis::convertToFrequencyData( const int segmentLength, const std::vecto
 				}else{
 					memcpy(dataSegments[iChan][counterSegment], &(itr->dataFile[iChan].data[index]), sizeof(double)*segmentLength);
 				}
-//#ifdef _DEBUG_WRITE
-//				std::ostringstream oss;
-//				oss << "sect_" << section << "_seg_" << counterSegment << "_chan_" << iChan << ".csv"; 
-//				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 < segmentLength; ++i ){
-//					ofs << std::setprecision(12) << std::scientific << dataSegments[iChan][counterSegment][i] << std::endl;
-//				}
-//				ofs.close();
-//#endif
+#ifdef _MTH5
+				m_segmentIndexToSectionIndex.insert(std::make_pair(counterSegment, section));
+#endif
 			}
 			// Start time
 			const int index1 = iSeg * shiftLength;
@@ -1219,19 +1199,6 @@ void Analysis::convertToFrequencyData( const int segmentLength, const std::vecto
 	for( int iChan = 0; iChan < numChannels; ++iChan ){
 		for( int iSeg = 0; iSeg < numSegmentsTotal; ++iSeg ){
 			Util::hanningWindow(segmentLength, dataSegments[iChan][iSeg]);
-//#ifdef _DEBUG_WRITE
-//			std::ostringstream oss;
-//			oss << "seg_" << iSeg << "_chan_" << iChan << "_hanning.csv"; 
-//			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 < segmentLength; ++i ){
-//				ofs << std::setprecision(12) << std::scientific << dataSegments[iChan][iSeg][i] << std::endl;
-//			}
-//			ofs.close();
-//#endif
 		}
 	}
 
@@ -1364,6 +1331,21 @@ void Analysis::calibrationCorrection(const int iChan, const int numSegmentsTotal
 		ftval[iSeg] *= calCorrFunc;
 	}
 
+#ifdef _MTH5
+	if (ptrControl->doesReadMTH5() && ptrControl->doesReadMTH5Filters()) {
+		OutputFiles* ptrOutputFiles = OutputFiles::getInstance();
+		const MTH5* ptrMTH5 = MTH5::getInstance();
+		for(int iSeg = 0; iSeg < numSegmentsTotal; ++iSeg) {
+			std::map<int, int>::const_iterator itrFind = m_segmentIndexToSectionIndex.find(iSeg);
+			if (itrFind == m_segmentIndexToSectionIndex.end()) {
+				ptrOutputFiles->writeErrorMessage("Segment " + Util::toString(iSeg) + " is not stored in the map");
+			}
+			const int iSection = itrFind->second;
+			ftval[iSeg] /= ptrMTH5->calcResponse(iSection, iChan, freq);
+		}
+	}
+#endif
+
 }
 
 // Calculate apparent resistivity

src/Analysis.cpp

@@ -35,6 +35,7 @@
 #include "RobustWeightThomson.h"
 #include "RobustWeightTukeysBiweights.h"
 #include <vector>
+#include <map>
 
 // Class of analysis
 class Analysis{
@@ -272,6 +273,11 @@ class Analysis{
 	// M-estimators
 	RobustWeight* m_robustWeight[2];
 
+#ifdef _MTH5
+	// Map segment index to sectionindex
+	std::map<int, int> m_segmentIndexToSectionIndex;
+#endif
+
 };
 
 #endif

src/Analysis.h

@@ -95,6 +95,7 @@ Control::Control() :
 	m_procedureType(Control::ORDINARY_REMOTE_REFERENCE),
 	m_readAtsBinary(false),
 	m_readMTH5(false),
+	m_readMTH5Filters(false),
 	m_readElogDualBinary(false),
 	m_readElogMTBinary(false),
 	m_timingEOFBasedDenoising(Control::BEFORE_DECIMATION),
@@ -210,16 +211,12 @@ Control::~Control(){
 }
 
 // Run analysis
-void Control::run(const bool outputToConsole){
+void Control::run(const bool outputToConsole) {
 
 	OutputFiles* ptrOutputFiles = OutputFiles::getInstance();
 	ptrOutputFiles->setOutputToConsole(outputToConsole);
 	readParameterFile();
 
-#ifdef _MTH5
-	(MTH5::getInstance())->createChannelResponses(m_filterInfoForMTH5);
-#endif
-
 	m_analysis->run(m_dataFileSets);
 
 	ptrOutputFiles->writeLogMessage("End " + Util::toString(CommonParameters::programName) );
@@ -281,6 +278,11 @@ bool Control::doesReadMTH5() const {
 	return m_readMTH5;
 }
 
+// Get flag specifing whether MTH5 filters are read or not
+bool Control::doesReadMTH5Filters() const {
+	return m_readMTH5Filters;
+}
+
 // Get flag specifing whether input file is ELOG-Dual binary file
 bool Control::doesReadElogDualBinary() const {
 	return m_readElogDualBinary;
@@ -427,26 +429,11 @@ int Control::getNumCalibrationFilesForMFS() const {
 	return static_cast<int>(m_calibrationFilesForMFS.size());
 }
 
-// Get number of filter info (pair of file name and path) for MTH5
-int Control::getNumFilterInfoMTH5() const {
-	return static_cast<int>(m_filterInfoForMTH5.size());
-}
-
 // Get name of calibration file for MFS
 std::string Control::getCalibrationFileNameForMFS(const int iFile) const {
 	return m_calibrationFilesForMFS[iFile];
 }
 
-// Get name of MTH5 file storing filter
-std::string Control::getFileNameForFilterInMTH5(const int i) const {
-	return m_filterInfoForMTH5[i].first;
-}
-
-// Get path of the filter in MTH5 file
-std::string Control::getPathFilterInMTH5(const int i) const {
-	return m_filterInfoForMTH5[i].second;
-}
-
 // Get numebur of calibration files
 int Control::getNumCalibrationFiles() const{
 	return static_cast<int>( m_calibrationFiles.size() );
@@ -863,18 +850,7 @@ void Control::readParameterFile(){
 		}
 #ifdef _MTH5
 		else if (line.find("MTH5_FILTERS") != std::string::npos) {
-			const int numChannels = getNumberOfChannels();
-			m_filterInfoForMTH5.clear();
-			m_filterInfoForMTH5.reserve(numChannels);
-			m_calibrationFiles.clear();
-			m_calibrationFiles.reserve(numChannels);
-			for (int iChan = 0; iChan < numChannels; ++iChan) {
-				std::string fileName;
-				std::string channelPath;
-				ifs >> fileName >> channelPath;
-				m_filterInfoForMTH5.push_back(std::make_pair(fileName, channelPath));
-				m_calibrationFiles.push_back(MTH5::getCalibrationFileName(iChan));
-			}
+			m_readMTH5Filters = true;
 		}
 		else if (line.find("MTH5") != std::string::npos) {
 			m_readMTH5 = true;
@@ -1748,15 +1724,8 @@ void Control::readParameterFile(){
 	}
 	if( doesReadMTH5() ) {
 		ptrOutputFiles->writeLogMessage("Read MTH5 files", false);
-		ptrOutputFiles->writeLogMessage("Filters stored in MTH5 files are used for the calibration", false);
-		ptrOutputFiles->writeLogMessage("     File       Path", false);
-		for (std::vector< std::pair< std::string, std::string> >::const_iterator itr = m_filterInfoForMTH5.begin(); itr != m_filterInfoForMTH5.end(); ++itr) {
-			std::ostringstream oss;
-			oss << std::setw(5) << "";
-			oss << itr->first;
-			oss << std::setw(5) << " ";
-			oss << itr->second;
-			ptrOutputFiles->writeLogMessage(oss.str(), false);
+		if (doesReadMTH5Filters()) {
+			ptrOutputFiles->writeLogMessage("Read MTH5 filters", false);
 		}
 	}
 	if (doesMakeCalibrationFileForMFS()) {

src/Control.cpp

@@ -295,6 +295,9 @@ class Control{
 	// Get flag specifing whether input file is MTH5 file
 	bool doesReadMTH5() const;
 
+	// Get flag specifing whether MTH5 filters are read or not
+	bool doesReadMTH5Filters() const;
+
 	// Get flag specifing whether input file is ELOG-Dual binary file
 	bool doesReadElogDualBinary() const;
 
@@ -373,12 +376,6 @@ class Control{
 	// Get name of calibration file for MFS
 	std::string getCalibrationFileNameForMFS(const int iFile) const;
 
-	// Get name of MTH5 file storing filter
-	std::string getFileNameForFilterInMTH5(const int i) const;
-
-	// Get path of the filter in MTH5 file
-	std::string getPathFilterInMTH5(const int i) const;
-
 	// Get numebur of calibration files
 	int getNumCalibrationFiles() const;
 
@@ -603,9 +600,6 @@ class Control{
 	// Calibration file for MFS
 	std::vector<std::string> m_calibrationFilesForMFS;
 
-	// Filter information in MTH5 files
-	std::vector< std::pair< std::string, std::string> > m_filterInfoForMTH5;
-
 	// Percentage of ommited data in subset deletion jackknife
 	double m_percentageOfOmmitedDataSubsetDeletionJackknife;
 
@@ -696,6 +690,9 @@ class Control{
 	// Flag specifing whether input file is MTH5 file
 	bool m_readMTH5;
 
+	// Flag specifing whether input file is MTH5 filters
+	bool m_readMTH5Filters;
+
 	// Flag specifing whether ELOG-Dual binary is read
 	bool m_readElogDualBinary;
 

src/Control.h

@@ -34,19 +34,23 @@
 #include <algorithm>
 #include <sstream>
 #include <iomanip>
+#include <assert.h>
 
 // Default constructer
 MTH5::MTH5():
-	m_numOfChannelRespones(0),
-	m_channelResponses(NULL)
+	m_numOfChannelRespones(0)
 {
 }
 
 // Destructer
 MTH5::~MTH5(){
-	if (m_channelResponses != NULL) {
-		delete[] m_channelResponses;
+	
+	for (std::vector<MTH5ChannelResponse*>::iterator itr = m_channelResponses.begin(); itr != m_channelResponses.end(); ++itr) {
+		if ( *itr != NULL) {
+			delete[] * itr;
+		}
 	}
+
 }
 
 // Return the instance of the class
@@ -99,32 +103,30 @@ std::string MTH5::getCalibrationFileName(const int channelIndex){
 
 }
 
-// Get all filters and combine filters into a complete channel response for each channel
-void MTH5::createChannelResponses(const std::vector< std::pair< std::string, std::string> >& fileNameAndPath) {
+// Calculate frequency response functions using the frequency response functions of all filter
+std::complex<double> MTH5::calcResponse(const int sectionIndex, const int channelIndex, const double freq) const {
 
-	if (fileNameAndPath.empty()) {
-		return;
-	}
-	m_numOfChannelRespones = static_cast<int>(fileNameAndPath.size());
-	m_channelResponses = new MTH5ChannelResponse[m_numOfChannelRespones];
-	int iChan(0);
-	for (std::vector< std::pair< std::string, std::string> >::const_iterator itr = fileNameAndPath.begin(); itr != fileNameAndPath.end(); ++itr, ++iChan) {
-		m_channelResponses[iChan].createChannelResponse(itr->first, itr->second);
-	}
+	assert(sectionIndex >= 0 && sectionIndex < m_channelResponses.size());
+	assert(channelIndex >= 0);
+
+	return m_channelResponses[sectionIndex][channelIndex].calcResponse(freq);
 
 }
 
-// Make calibration files using the requency response functions of all filter 
-void MTH5::makeCalibrationFile(const int channelIndex, const std::vector<double>& freqs) const {
+// Read filters for indivial sections and channels
+void MTH5::readFiltersAll(const int numChannels, const std::vector<CommonParameters::DataFileSet>& dataFileSets) {
 
-	if (m_channelResponses[channelIndex].isAppliedForwardly()) {
-		return;
-	}
-	
-	if (channelIndex < 0 && channelIndex >= m_numOfChannelRespones) {
-		OutputFiles::getInstance()->writeErrorMessage("Channel index is out of range: " + Util::toString(channelIndex));
+	m_numOfChannelRespones = numChannels;
+
+	for(std::vector<CommonParameters::DataFileSet>::const_iterator itr = dataFileSets.begin(); itr != dataFileSets.end(); ++itr){
+		const std::vector<CommonParameters::DataFile>& dataFileList = itr->dataFile;
+		MTH5ChannelResponse* channelResponses = new MTH5ChannelResponse[numChannels];
+		for (int iChan = 0; iChan < numChannels; ++iChan) {
+			const std::string fileName = dataFileList[iChan].fileName;
+			const std::string channelPath = dataFileList[iChan].mth5GroupName;
+			channelResponses[iChan].createChannelResponse(fileName, channelPath);
+		}
+		m_channelResponses.push_back(channelResponses);
 	}
-	
-	m_channelResponses[channelIndex].makeCalibrationFile(getCalibrationFileName(channelIndex), channelIndex, freqs);
 
 }

src/MTH5.cpp

@@ -34,6 +34,7 @@
 #include <H5Cpp.h>
 
 #include "MTH5ChannelResponse.h"
+#include "CommonParameters.h"
 
 // Class of MTH5 file
 class MTH5{
@@ -49,22 +50,19 @@ class MTH5{
 	// Get name of the calibration file name made from the channel responses 
 	static std::string getCalibrationFileName(const int channelIndex);
 
-	// Get all filters and combine filters into a complete channel response for each channel
-	void createChannelResponses(const std::vector< std::pair< std::string, std::string> >& fileNameAndPath);
+	// Read filters for indivial sections and channels
+	void readFiltersAll(const int numChannels, const std::vector<CommonParameters::DataFileSet>& dataFileSets);
 
-	//// Correct frequency response functions using the requency response functions of all filter
-	//void correctResponse(const int channelIndex, const double freq, const double samplingFreq, std::complex<double>& response) const;
-
-	// Make calibration files using the requency response functions of all filter 
-	void makeCalibrationFile(const int channelIndex, const std::vector<double>& freqs) const;
+	// Calculate frequency response functions using the frequency response functions of all filter
+	std::complex<double> calcResponse(const int sectionIndex, const int channelIndex, const double freq) const;
 
 private:
 
 	// Number of channel responses
 	int m_numOfChannelRespones;
 
 	// List of channel response (combination of all filters)
-	MTH5ChannelResponse* m_channelResponses;
+	std::vector<MTH5ChannelResponse*> m_channelResponses;
 
 	// Constructer
 	MTH5();