From 1a219839034e9b11a4771fb84c90d4a2667365ce Mon Sep 17 00:00:00 2001
From: Roland Stigge <stigge@antcom.de>
Date: Sat, 16 Feb 2019 23:37:21 +0100
Subject: Added IFFT + AutoCorrelation

---
 fft.cpp | 73 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 68 insertions(+), 5 deletions(-)

(limited to 'fft.cpp')

diff --git a/fft.cpp b/fft.cpp
index 42b0acf..7d76e23 100644
--- a/fft.cpp
+++ b/fft.cpp
@@ -43,7 +43,7 @@ RIT::FFT::FFT(int size, bool halfOnly): mSize(size), order(size), expLUT(size/2)
 }
 
 
-std::vector<std::complex<double>> RIT::FFT::operator()(const std::vector<std::complex<double>> &v) {
+std::vector<std::complex<double>> RIT::FFT::operator()(const std::vector<std::complex<double>> &v) const {
 	if (v.size() != mSize)
 		throw std::length_error("Bad input size");
 
@@ -64,7 +64,7 @@ RIT::FFT& RIT::FFT::SetHalfOnly(bool enable) {
 	return *this;
 }
 
-int RIT::FFT::bitreverse(int i) {
+int RIT::FFT::bitreverse(int i) const {
 	int size{mSize};
 	int result{0};
 
@@ -78,7 +78,7 @@ int RIT::FFT::bitreverse(int i) {
 	return result;
 }
 
-void RIT::FFT::reorder(const std::vector<std::complex<double>>& src, std::vector<std::complex<double>>& dst) {
+void RIT::FFT::reorder(const std::vector<std::complex<double>>& src, std::vector<std::complex<double>>& dst) const {
 	int size = src.size();
 
 	dst.resize(size);
@@ -94,7 +94,7 @@ void RIT::FFT::reorder(const std::vector<std::complex<double>>& src, std::vector
 // Because of Nyquist theorem, N samples means 
 // only first N/2 FFT results in X[] are the answer.
 // (upper half of X[] is a reflection with no new information).
-void RIT::FFT::fft_recursive(std::vector<std::complex<double>>::iterator X, int N) {
+void RIT::FFT::fft_recursive(std::vector<std::complex<double>>::iterator X, int N) const {
 	if(N > 2) {
 		fft_recursive(X,     N/2);   // recurse even items
 		fft_recursive(X+N/2, N/2);   // recurse odd  items
@@ -111,7 +111,7 @@ void RIT::FFT::fft_recursive(std::vector<std::complex<double>>::iterator X, int
 }
 
 // Same as fft_recursive, but results in only half the result due to symmetry in real input
-void RIT::FFT::fft_half(std::vector<std::complex<double>>::iterator X, int N) {
+void RIT::FFT::fft_half(std::vector<std::complex<double>>::iterator X, int N) const {
 	if(N > 2) {
 		fft_recursive(X,     N/2);   // recurse even items
 		fft_recursive(X+N/2, N/2);   // recurse odd  items
@@ -121,3 +121,66 @@ void RIT::FFT::fft_half(std::vector<std::complex<double>>::iterator X, int N) {
 		X[k] += expLUT[k * mSize / N] * X[k+N/2];
 	}
 }
+
+struct RIT::IFFT::Impl {
+public:
+	using transform_function = std::complex<double> (*)(const std::complex<double>&);
+
+	Impl(int size): mFft(std::make_shared<RIT::FFT>(size)) {}
+	Impl(std::shared_ptr<RIT::FFT> fft): mFft(fft) {}
+	~Impl(){}
+	std::shared_ptr<RIT::FFT> mFft;
+};
+
+RIT::IFFT::IFFT(int size): mImpl(std::make_unique<RIT::IFFT::Impl>(size))
+{
+}
+
+RIT::IFFT::IFFT(std::shared_ptr<RIT::FFT> fft): mImpl(std::make_unique<RIT::IFFT::Impl>(fft))
+{
+}
+
+RIT::IFFT::~IFFT()
+{
+}
+
+std::vector<std::complex<double>> RIT::IFFT::operator()(const std::vector<std::complex<double>> &v) const
+{
+	std::vector<std::complex<double>> input(v.size());
+
+	std::transform(std::begin(v), std::end(v), std::begin(input), static_cast<Impl::transform_function>(std::conj));
+
+	auto result = (*mImpl->mFft)(input);
+	
+	const double factor = 1.0 / v.size();
+	std::transform(std::begin(result), std::end(result), std::begin(result),
+		       [=](const auto& x){return std::conj(x) * factor;});
+
+	return result;
+}
+
+struct RIT::AutoCorrelation::Impl {
+public:
+	Impl(int size): mFft(std::make_shared<RIT::FFT>(size)), mIfft(mFft) {}
+	std::shared_ptr<RIT::FFT> mFft;
+	RIT::IFFT mIfft;
+};
+
+RIT::AutoCorrelation::AutoCorrelation(int size): mImpl(std::make_unique<RIT::AutoCorrelation::Impl>(size))
+{
+}
+
+RIT::AutoCorrelation::~AutoCorrelation()
+{
+}
+
+std::vector<std::complex<double>> RIT::AutoCorrelation::operator()(const std::vector<std::complex<double>> &v) const
+{
+	auto result = (*mImpl->mFft)(v);
+
+	std::transform(std::begin(result), std::end(result), std::begin(result),
+		       [=](const auto& x){return x * std::conj(x);});
+	
+	return mImpl->mIfft(result);
+}
+
-- 
cgit v1.2.3