Initial NBData, added FFT
This commit is contained in:
NaifBanana
@@ -0,0 +1,252 @@
|
||||
#pragma once
|
||||
#ifndef _NB_DATA_FFT
|
||||
#define _NB_DATA_FFT
|
||||
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
#include <NBCore/Errors.hpp>
|
||||
|
||||
#ifndef PI_VALUE
|
||||
#define PI_VALUE 3.1415926535
|
||||
#endif
|
||||
|
||||
/*! @file FFT.hpp */
|
||||
|
||||
namespace nb {
|
||||
|
||||
using namespace std::complex_literals;
|
||||
|
||||
constexpr double operator""_pi(long double val) {
|
||||
return val*PI_VALUE;
|
||||
}
|
||||
constexpr double operator""_pi(unsigned long long val) {
|
||||
return val*PI_VALUE;
|
||||
}
|
||||
|
||||
constexpr float operator""_pif(long double val) {
|
||||
return val*PI_VALUE;
|
||||
}
|
||||
constexpr float operator""_pif(unsigned long long val) {
|
||||
return val*PI_VALUE;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool isPowerOf2(const T& x) {
|
||||
bool ret = false;
|
||||
unsigned int mask=1;
|
||||
for (int i=0; i<(8*sizeof(std::declval<T>())); ++i) {ret ^= x&mask; mask<<=1;}
|
||||
return ret;
|
||||
}
|
||||
|
||||
/*!
|
||||
@brief Calculates the complex-valued twiddle factor given the `phase`
|
||||
@param phase A double representing the phase, or `(k*n)/N`, of a DFT,
|
||||
in units of 2*Pi.
|
||||
@return A std::complex<T> representing the complex twiddle factor of
|
||||
the DFT at that `phase`, calculated via `exp(-2*Pi*i*phase)`.
|
||||
|
||||
Calculates the complex-valued twiddle factor given the `phase` in units
|
||||
of the ratio around the circle (i.e. in units of 2*Pi), calculated via
|
||||
`exp(-2*Pi*i*phase)`. Note the minus sign, as this function is for the
|
||||
*forward* DFT.
|
||||
*/
|
||||
template<typename T>
|
||||
std::complex<T> calculateTwiddleFactor(double phase) {
|
||||
int sign = -1;
|
||||
while(phase>0.5) { phase-=0.5; sign*=-1; }
|
||||
return static_cast<std::complex<T>>(std::exp(2_pi*(sign)*(1i)*phase));
|
||||
}
|
||||
|
||||
/*!
|
||||
Stores the twiddle factors for Radix-2 DFT matrices for a given data type
|
||||
(i.e. stores std::complex<T> values), which can be retrieved using `TwiddleFactor<T>::get`.
|
||||
*/
|
||||
template<typename T>
|
||||
class TwiddleFactors2 {
|
||||
protected:
|
||||
static std::unordered_map<unsigned int, std::complex<T>> factors;
|
||||
|
||||
public:
|
||||
/*!
|
||||
@brief Retrieves the twiddle factor for a Radix-2 DFT
|
||||
@param numerator The numerator (i.e. `k*n` of a DFT matrix) part of the phase
|
||||
when defined in units of 2*Pi.
|
||||
@param log2N The `log2(denominator)` (i.e. `/N` of a DFT matrix) part of the phase
|
||||
when defined in units of 2*Pi.
|
||||
@return A std::complex<T> for the value of the `numerator / (2^log2N)` twiddle factor
|
||||
*/
|
||||
static std::complex<T> get(unsigned int numerator, unsigned int log2N);
|
||||
};
|
||||
|
||||
template<typename T>
|
||||
std::complex<T>* fft2_impl(
|
||||
const std::complex<T>* source,
|
||||
std::complex<T>* dest,
|
||||
unsigned int log2N,
|
||||
unsigned int stride=1
|
||||
) {
|
||||
constexpr std::complex<T> J(0, 1);
|
||||
switch(log2N) {
|
||||
case 0:
|
||||
dest[0] = source[0];
|
||||
break;
|
||||
|
||||
case 1:
|
||||
dest[0] = source[0] + source[stride];
|
||||
dest[1] = source[0] - source[stride];
|
||||
break;
|
||||
|
||||
case 2:
|
||||
dest[0] = source[0] + source[stride] + source[2*stride] + source[3*stride];
|
||||
dest[1] = source[0] - J*source[stride] - source[2*stride] + J*source[3*stride];
|
||||
dest[2] = source[0] - source[stride] + source[2*stride] - source[3*stride];
|
||||
dest[3] = source[0] + J*source[stride] - source[2*stride] - J*source[3*stride];
|
||||
break;
|
||||
|
||||
default:
|
||||
const unsigned int halfway = 1<<(log2N-1);
|
||||
fft2_impl(source, dest, log2N-1, 2*stride);
|
||||
fft2_impl(&(source[stride]), &(dest[halfway]), log2N-1, 2*stride);
|
||||
std::complex<T> twiddle;
|
||||
std::complex<T> p, q;
|
||||
for (unsigned int i=0; i < halfway; ++i) {
|
||||
twiddle = TwiddleFactors2<T>::get(i, log2N);
|
||||
p = dest[i];
|
||||
q = dest[halfway+i]*twiddle;
|
||||
dest[i] = p+q;
|
||||
dest[halfway+i] = p-q;
|
||||
}
|
||||
break;
|
||||
}
|
||||
return dest;
|
||||
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
std::complex<T>* dft_impl(
|
||||
const std::complex<T>* source,
|
||||
std::complex<T>* dest,
|
||||
unsigned int N
|
||||
) {
|
||||
for (int k=0; k < N; ++k) {
|
||||
dest[k] = 0;
|
||||
for (int n = 0; n < N; ++n) {
|
||||
dest[k] += source[n] * calculateTwiddleFactor<T>(double(n*k)/double(N));
|
||||
}
|
||||
}
|
||||
return dest;
|
||||
}
|
||||
|
||||
/*!
|
||||
@brief Computes the minimal Radix-2 DFT on a std::vector<std::complex<T>>
|
||||
@param values A std::vector<std::complex<T>> representing the time domain values
|
||||
@return A std::vector<std::complex<T>>, representing the minimal Radix-2 DFT of `values`
|
||||
|
||||
Computes the Radix-2 DFT on a std::vector<std::complex<T>>. If the input vector
|
||||
`values` is not a power of 2, then the number of points in the DFT (and thus the
|
||||
size of the returned frequency domain vector) is determined by
|
||||
`2^ceil(log2(values.size()))`
|
||||
*/
|
||||
template<typename T>
|
||||
std::vector<std::complex<T>> fft2(const std::vector<std::complex<T>>& values) {
|
||||
const auto N = values.size();
|
||||
unsigned int log2N=0;
|
||||
while(1<<log2N < N) {log2N++;}
|
||||
std::vector<std::complex<T>> temp(1<<log2N, 0);
|
||||
std::vector<std::complex<T>> ret = temp;
|
||||
for (int i=0; i < N; ++i) { temp[i] = values[i]; }
|
||||
std::complex<T>* fft_results = fft2_impl(temp.data(), ret.data(), log2N);
|
||||
return ret;
|
||||
}
|
||||
|
||||
/*!
|
||||
@brief Computes the minimal Radix-2 DFT on a `std::vector<T>`
|
||||
@param values A std::vector<T> representing the time domain values with packed real and
|
||||
imaginary parts
|
||||
@return A `std::vector<T>`, representing the minimal Radix-2 DFT of `values` with packed
|
||||
real and imaginary parts
|
||||
@throw std::string Thrown when the size of the vector is not divisible by 2 (i.e.
|
||||
there are unmatched real/imaginary pairs).
|
||||
|
||||
Computes the Radix-2 DFT on a std::vector<std::complex<T>>. Equivalent to
|
||||
`fft2(const std::vector<std::complex<T>>)` if `values[2*i]` and `values[2*i+1]`
|
||||
were mapped to the real and imaginary parts of the complex vector.
|
||||
*/
|
||||
template<typename T>
|
||||
std::vector<T> fft2(const std::vector<T>& values) {
|
||||
if (values.size()&1) {
|
||||
THROW(Error<>(Error<>::VALUE_ERROR,
|
||||
"Vector must be of even length to represent complex values."
|
||||
));
|
||||
}
|
||||
std::complex<T> complex_input(values.data(), values.size()/2);
|
||||
std::complex<T> complex_ret = fft2(complex_input);
|
||||
return std::vector<T>(complex_ret.data(), complex_ret.size()*2);
|
||||
}
|
||||
|
||||
/*!
|
||||
@brief Computes the DFT on a std::vector<std::complex<T>>
|
||||
@param values A std::vector<std::complex<T>> representing the time domain values
|
||||
@return A std::vector<std::complex<T>> of the DFT of `values` in the frequency domain
|
||||
|
||||
Computes the DFT of the std::vector<std::complex<T>> `values`. Does so using the
|
||||
naive O(N^2) matrix-multiplication method.
|
||||
*/
|
||||
template<typename T>
|
||||
std::vector<std::complex<T>> dft(const std::vector<std::complex<T>>& x) {
|
||||
std::vector<std::complex<T>> ret = x;
|
||||
std::complex<T>* fft_results = dft_impl(x.data(), ret.data(), x.size());
|
||||
return ret;
|
||||
}
|
||||
|
||||
/*!
|
||||
@brief Computes the DFT on a `std::vector<T>`
|
||||
@param values A `std::vector<T>` representing the time domain values with packed real and
|
||||
imaginary parts
|
||||
@return A `std::vector<T>` of the DFT of `values` in the frequency domain with
|
||||
packed real and imaginary parts
|
||||
@throw std::string Thrown when the size of the vector is not divisible by 2 (i.e.
|
||||
there are unmatched real/imaginary pairs).
|
||||
|
||||
Computes the DFT on a `std::vector<T>`. Equivalent to
|
||||
`dft(const std::vector<std::complex<T>>)` if `values[2*i]` and `values[2*i+1]`
|
||||
were mapped to the real and imaginary parts of the complex vector.
|
||||
*/
|
||||
template<typename T>
|
||||
std::vector<T> dft(const std::vector<T>& values) {
|
||||
if (values.size()&1) {
|
||||
THROW(Error<>(Error<>::VALUE_ERROR,
|
||||
"Vector must be of even length to represent complex values."
|
||||
));
|
||||
}
|
||||
std::complex<T> complex_input(values.data(), values.size()/2);
|
||||
std::complex<T> complex_ret = dft(complex_input);
|
||||
return std::vector<T>(complex_ret.data(), complex_ret.size()*2);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
std::unordered_map<unsigned int, std::complex<T>> TwiddleFactors2<T>::factors = {};
|
||||
|
||||
template<typename T>
|
||||
std::complex<T> TwiddleFactors2<T>::get(unsigned int a, unsigned int log2N) {
|
||||
while(a && !(a&1)) { a>>=1; log2N--; }
|
||||
const unsigned int N = 1<<log2N;
|
||||
const double phase = double(a)/double(N);
|
||||
const unsigned int key = a*N;
|
||||
try {
|
||||
return factors.at(key);
|
||||
} catch (const std::out_of_range& e) {
|
||||
if (a > N>>1 ) {
|
||||
return (factors[key] = std::conj(TwiddleFactors2<T>::get(N-a, log2N)));
|
||||
} else {
|
||||
return (factors[key] = calculateTwiddleFactor<T>(phase));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace n
|
||||
|
||||
#endif // _NB_DATA_FFT
|
||||
Reference in New Issue
Block a user