From f2ad3bd04ab96b57097e8c5c10896533799db89d Mon Sep 17 00:00:00 2001
From: Robert <robert.trololo@gmail.com>
Date: Thu, 29 Apr 2021 23:21:52 +0200
Subject: [PATCH] added mertz method

---
 src/FFT.cpp  | 9 +++++++--
 src/FFT.hpp  | 3 ++-
 src/main.cpp | 9 +++++++--
 3 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/FFT.cpp b/src/FFT.cpp
index 2171da9..140e334 100644
--- a/src/FFT.cpp
+++ b/src/FFT.cpp
@@ -58,7 +58,8 @@ FFT(const std::vector<double>::const_iterator& begin,
 	const std::vector<double>::const_iterator& end,
 	size_t sampleRate,
 	double minFreq, double maxFreq,
-	unsigned int zeropadding)
+	unsigned int zeropadding,
+	bool mertz)
 {
 	std::vector<double> signal(begin, end);
 	size_t N = signal.size();
@@ -89,7 +90,11 @@ FFT(const std::vector<double>::const_iterator& begin,
 
 	for (int k = freq / freqRes; freq < nyquistLimit && freq < maxFreq; k++)
 	{
-		output.push_back(std::make_pair(freq, 2.0f * std::abs(spectrum[k]) / (double)N));
+		if(!mertz)
+			output.push_back(std::make_pair(freq, 2.0f * std::abs(spectrum[k]) / (double)N));
+		else
+			output.push_back(std::make_pair(freq, 2.0f * (spectrum[k] * std::exp(-1i * std::arg(spectrum[k]))).real() / (double)N));
+
 		freq += freqRes;
 	}
 
diff --git a/src/FFT.hpp b/src/FFT.hpp
index 0a04842..25c0e8f 100644
--- a/src/FFT.hpp
+++ b/src/FFT.hpp
@@ -14,6 +14,7 @@ extern std::vector<std::pair<double, double>> FFT(const std::vector<double>::con
 	const std::vector<double>::const_iterator& end,
 	size_t sampleRate,
 	double minFreq, double maxFreq,
-	unsigned int zeropadding);
+	unsigned int zeropadding,
+	bool mertz);
 
 extern void SetWindowFunction(WindowFunctions func, unsigned int width);
\ No newline at end of file
diff --git a/src/main.cpp b/src/main.cpp
index 8daff5a..0e76896 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -26,6 +26,7 @@ struct Settings {
 	double minFreq, maxFreq;
 	unsigned int analyzeChannel;
 	unsigned int zeropadding;
+	bool mertz;
 	WindowFunctions window;
 };
 
@@ -72,7 +73,8 @@ int main(int argc, char** argv)
 						audioFile.samples[c-1].cend(), 
 						sampleRate,
 						setts.minFreq, setts.maxFreq,
-						setts.zeropadding
+						setts.zeropadding,
+						setts.mertz
 					);
 
 				output[chName] = nlohmann::json::array();
@@ -96,7 +98,8 @@ int main(int argc, char** argv)
 							),
 							sampleRate,
 							setts.minFreq, setts.maxFreq,
-							setts.zeropadding
+							setts.zeropadding,
+							setts.mertz
 						);
 
 					output[chName].push_back({
@@ -140,6 +143,7 @@ Settings Parse(int argc, char** argv)
 			("f,frequency", "Defines the frequency range of the output spectrum (Default: all the frequencies)", cxxopts::value<std::vector<double>>())
 			("p,pad", "Add extra zero-padding. By default, the program will pad the signals with 0s until the number of samples is a power of 2 (this would be equivalent to -p 1). With this option you can tell the program to instead pad until the power of 2 after the next one (-p 2) etc. This increases frequency resolution", cxxopts::value<unsigned int>())
 			("w,window", "Specify the window function used (rectangle (default), von-hann, gauss, triangle, blackman (3-term))", cxxopts::value<std::string>()->default_value("rectangle"))
+			("mertz", "Use the Mertz method to phase-correct the complex Fourier spectrum")
 			("m,mono", "Analyze only the given channel", cxxopts::value<unsigned int>()->default_value("0"))
 			("files", "Files to fourier transform", cxxopts::value<std::vector<std::filesystem::path>>())
 			("h,help", "Print usage")
@@ -176,6 +180,7 @@ Settings Parse(int argc, char** argv)
 		setts.splitInterval = (result.count("interval") ? result["interval"].as<float>() : 0.0f);
 		setts.analyzeChannel = (result.count("mono") ? result["mono"].as<unsigned int>() : 0);
 		setts.zeropadding = (result.count("pad") ? result["pad"].as<unsigned int>() : 1);
+		setts.mertz = (result.count("mertz") ? true : false);
 
 		if (!result.count("window"))
 		{