#include "DiscreteFourierTransform.hh" #include <bits/stdc++.h> template <class T, class E> DiscreteFourierTransform<T, E>::DiscreteFourierTransform(vector<T> signal_, double sampleSpacing_) { signal = signal_; sampleSpacing = sampleSpacing_; } template <class T, class E> void DiscreteFourierTransform<T, E>::compute(bool computeFFT) { if (computeFFT) { const int N = signal.size(); const int powerOf2 = DiscreteFourierTransform<T, E>::highestPowerof2(N); // Zero padding is a simple concept; it simply refers to adding zeros to end of a time-domain signal to increase its length to reach power of 2. while (signal.size() < powerOf2) signal.push_back((T)0); phasors = DiscreteFourierTransform<T, E>::fft(signal); } else phasors = DiscreteFourierTransform<T, E>::dft(signal); } template <class T, class E> std::vector<E> DiscreteFourierTransform<T, E>::createTestPoints(int N, double sampleSpacing) { std::vector<E> out; const E t = (E)sampleSpacing; for (int i = 0; i < N; i++) { out.push_back((E)t * i); } return out; } template <class T, class E> std::vector<E> DiscreteFourierTransform<T, E>::createTestFunction(std::vector<E> x) { const int N = x.size(); std::vector<E> out; for (int i = 0; i < N; i++) { const E val = (E)sin(50.0 * 2.0 * M_PI * x[i]) + 0.5 * sin(80.0 * 2.0 * M_PI * x[i]); out.push_back(val); } return out; } /** * @brief We use Cooley-Tukey FFT Algorithm if the size of the signal is a power of 2 otherwise you should DFT brut. In The future we will use bluestein algorithm instead of DFT. * Cooley-Tukey FFT Algorithms: http://people.scs.carleton.ca/~maheshwa/courses/5703COMP/16Fall/FFT_Report.pdf * * @tparam T template type of input data * @tparam E template type of output data * @param sig the signal as an input * @return std::vector<std::complex<E>> the phasors (array of complex number) */ template <class T, class E> std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::fft(std::vector<T> sig) { const int N = sig.size(); if (N == 1) { std::vector<std::complex<E>> out; const std::complex<E> temp((E)sig[0], (E)0); out.push_back(temp); return out; } const std::complex<E> WN = (complex<E>)std::polar(1.0, 2 * M_PI / N); std::complex<E> W((E)1, (E)0); // divide and conquer: // Recurse: all even samples std::vector<std::complex<E>> x_evens = fft(getEven(sig)); // Recurse: all odd samples std::vector<std::complex<E>> x_odds = fft(getOdd(sig)); // Now, combine and perform N/2 operations! std::complex<E> zeroComplex((E)0, (E)0); std::vector<std::complex<E>> x(N, zeroComplex); for (int k = 0; k < N / 2; k++) { x[k] = x_evens[k] + W * x_odds[k]; x[k + (N / 2)] = x_evens[k] - W * x_odds[k]; W = W * WN; } return x; } template <class T, class E> std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::dft(std::vector<T> x) { const int N = x.size(); std::complex<E> zeroComplex((E)0, (E)0); std::vector<std::complex<E>> out(N, zeroComplex); for (int k = 0; k < N; k++) { for (int n = 0; n < N; n++) { const std::complex<E> expVal = (complex<E>)std::polar(1.0, -2 * k * n * M_PI / N); out[k] += ((E)x[n]) * expVal; } } return out; } template <class T, class E> std::vector<E> DiscreteFourierTransform<T, E>::computeDSP(std::vector<std::complex<E>> x) { const int N = x.size(); std::vector<E> out; for (int k = 0; k < N / 2; k++) { const E magnitude = (E)pow(abs(x[k]), 1); const E dsp = (E)(2.0 / N) * magnitude; out.push_back(dsp); } return out; } template <class T, class E> std::vector<E> DiscreteFourierTransform<T, E>::getFreq(std::vector<std::complex<E>> x, double frequency_) { const int N = x.size(); std::vector<E> out; for (int k = 0; k < N / 2; k++) { const E freq = (E)k * frequency_ / N; out.push_back(freq); } return out; } template <class T, class E> std::vector<E> DiscreteFourierTransform<T, E>::getPeriods(std::vector<std::complex<E>> x, double frequency_) { const int N = x.size(); const double stepPeriod = N / sampleSpacing_; std::vector<E> out; for (int k = 0; k < N / 2; k++) { const E period = (E)k * frequency_; out.push_back(period); } return out; } template <class T, class E> std::vector<T> DiscreteFourierTransform<T, E>::getOdd(std::vector<T> x) { std::vector<T> odd; for (int i = 0; i < x.size(); i++) { if (i % 2 != 0) odd.push_back(x[i]); } return odd; } template <class T, class E> std::vector<T> DiscreteFourierTransform<T, E>::getEven(std::vector<T> x) { std::vector<T> even; for (int i = 0; i < x.size(); i++) { if (i % 2 == 0) even.push_back(x[i]); } return even; } template <class T, class E> bool DiscreteFourierTransform<T, E>::isPowerOfTwo(int N) { return (N & (N - 1)) == 0; } template <class T, class E> int DiscreteFourierTransform<T, E>::highestPowerof2(int n) { int p = n; while (!DiscreteFourierTransform<T, E>::isPowerOfTwo(p)) { p += 1; } return p; }