85 lines
2.5 KiB
C++
85 lines
2.5 KiB
C++
#pragma once
|
|
|
|
/**
|
|
* @file fft.hpp FFT and related facilities.
|
|
*/
|
|
|
|
#include "sample.hpp"
|
|
#include <complex>
|
|
#include <cstddef>
|
|
#include <vector>
|
|
|
|
#ifndef M_PI
|
|
#define M_PI 3.141592653589793
|
|
#endif
|
|
|
|
#define M_TAU (2.0 * M_PI)
|
|
|
|
namespace da {
|
|
|
|
// With g++ optimization -fcx-limited-range should be used for 5x performance boost.
|
|
|
|
// Based on the description of Volodymyr Myrnyy in
|
|
// http://www.dspdesignline.com/showArticle.jhtml?printableArticle=true&articleId=199903272
|
|
template<unsigned P, typename T> struct DanielsonLanczos {
|
|
static void apply(std::complex<T>* data) {
|
|
constexpr std::size_t N = 1 << P;
|
|
constexpr std::size_t M = N / 2;
|
|
// Compute even and odd halves
|
|
DanielsonLanczos<P - 1, T>().apply(data);
|
|
DanielsonLanczos<P - 1, T>().apply(data + M);
|
|
// Combine the results
|
|
std::complex<T> w(1.0);
|
|
for (std::size_t i = 0; i < M; ++i) {
|
|
const std::complex<T> temp = data[i + M] * w;
|
|
data[M + i] = data[i] - temp;
|
|
data[i] += temp;
|
|
w *= std::polar(1.0, -M_TAU / N);
|
|
}
|
|
}
|
|
};
|
|
|
|
template<typename T> struct DanielsonLanczos<0, T> { static void apply(std::complex<T>*) {} };
|
|
|
|
/** Perform FFT on data. **/
|
|
template<unsigned P, typename T> void fft(std::complex<T>* data) {
|
|
// Perform bit-reversal sorting of sample data.
|
|
constexpr std::size_t N = 1 << P;
|
|
std::size_t j = 0;
|
|
for (std::size_t i = 0; i < N; ++i) {
|
|
if (i < j) std::swap(data[i], data[j]);
|
|
std::size_t m = N / 2;
|
|
while (m > 1 && m <= j) { j -= m; m >>= 1; }
|
|
j += m;
|
|
}
|
|
// Do the actual calculation
|
|
DanielsonLanczos<P, T>::apply(data);
|
|
}
|
|
|
|
/** Perform FFT on data from floating point iterator, windowing the input. **/
|
|
template<unsigned P, typename InIt, typename Window> std::vector<std::complex<float> > fft(InIt begin, Window window) {
|
|
std::vector<std::complex<float> > data(1 << P);
|
|
// Perform bit-reversal sorting of sample data.
|
|
constexpr std::size_t N = 1 << P;
|
|
std::size_t j = 0;
|
|
for (std::size_t i = 0; i < N; ++i) {
|
|
data[j] = *begin++ * window[i];
|
|
std::size_t m = N / 2;
|
|
while (m > 1 && m <= j) { j -= m; m >>= 1; }
|
|
j += m;
|
|
}
|
|
// Do the actual calculation
|
|
DanielsonLanczos<P, float>::apply(&data[0]);
|
|
return data;
|
|
}
|
|
|
|
template<unsigned P, typename T> void ifft(std::complex<T>* data) {
|
|
constexpr std::size_t N = 1 << P;
|
|
for (std::size_t i = 0; i < N; ++i) data[i] = std::conj(data[i]); // Invert phase so that we can use FFT to do IFFT
|
|
fft<P>(data);
|
|
constexpr std::complex<T> scale(1.0/N, 0.0);
|
|
for (std::size_t i = 0; i < N; ++i) data[i] = scale * std::conj(data[i]); // Invert back, and apply IFFT scaling
|
|
}
|
|
|
|
}
|