Compare commits
No commits in common. "nbdata" and "main" have entirely different histories.
@ -1,70 +0,0 @@
|
|||||||
cmake_minimum_required(VERSION 3.10)
|
|
||||||
project(NBData VERSION 0.1)
|
|
||||||
|
|
||||||
include(GNUInstallDirs)
|
|
||||||
include(CMakePackageConfigHelpers)
|
|
||||||
|
|
||||||
toAbsolutePath(NB_DATA_SOURCE
|
|
||||||
./src/FFT.cpp
|
|
||||||
)
|
|
||||||
toAbsolutePath(NB_DATA_INCLUDE
|
|
||||||
./include/NBData/FFT.hpp
|
|
||||||
)
|
|
||||||
|
|
||||||
set(NB_DATA_SOURCE ${NB_DATA_SOURCE} PARENT_SCOPE)
|
|
||||||
set(NB_DATA_INCLUDE ${NB_DATA_INCLUDE} PARENT_SCOPE)
|
|
||||||
|
|
||||||
add_library(NBData ${NB_DATA_SOURCE})
|
|
||||||
add_library(NBEngine::Data ALIAS NBData)
|
|
||||||
target_link_libraries(NBData
|
|
||||||
NBCore
|
|
||||||
)
|
|
||||||
target_include_directories(NBData
|
|
||||||
PUBLIC "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>"
|
|
||||||
PUBLIC "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
|
|
||||||
)
|
|
||||||
export(
|
|
||||||
TARGETS NBData
|
|
||||||
FILE "${CMAKE_BINARY_DIR}/cmake/NBDataTargets.cmake"
|
|
||||||
NAMESPACE NBEngine::
|
|
||||||
)
|
|
||||||
configure_package_config_file(
|
|
||||||
"NBDataConfig.cmake.in"
|
|
||||||
"${CMAKE_BINARY_DIR}/cmake/NBDataConfig.cmake"
|
|
||||||
INSTALL_DESTINATION "${CMAKE_INSTALL_PREFIX}/cmake"
|
|
||||||
PATH_VARS CMAKE_INSTALL_LIBDIR
|
|
||||||
)
|
|
||||||
write_basic_package_version_file(
|
|
||||||
"${CMAKE_BINARY_DIR}/cmake/NBDataConfigVersion.cmake"
|
|
||||||
COMPATIBILITY AnyNewerVersion
|
|
||||||
)
|
|
||||||
|
|
||||||
if (NBENGINE_INSTALL)
|
|
||||||
message("Installing NBData to ${CMAKE_INSTALL_PREFIX}")
|
|
||||||
install(
|
|
||||||
TARGETS NBData
|
|
||||||
EXPORT NBDataTargets
|
|
||||||
LIBRARY DESTINATION lib
|
|
||||||
ARCHIVE DESTINATION lib
|
|
||||||
INCLUDES DESTINATION include
|
|
||||||
)
|
|
||||||
install(
|
|
||||||
DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/include/NBData"
|
|
||||||
DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}"
|
|
||||||
FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp"
|
|
||||||
)
|
|
||||||
install(
|
|
||||||
EXPORT NBDataTargets
|
|
||||||
DESTINATION "${CMAKE_INSTALL_PREFIX}/cmake"
|
|
||||||
NAMESPACE NBEngine::
|
|
||||||
)
|
|
||||||
install(FILES
|
|
||||||
"${CMAKE_BINARY_DIR}/cmake/NBDataConfig.cmake"
|
|
||||||
"${CMAKE_BINARY_DIR}/cmake/NBDataConfigVersion.cmake"
|
|
||||||
DESTINATION "${CMAKE_INSTALL_PREFIX}/CMake"
|
|
||||||
)
|
|
||||||
endif()
|
|
||||||
|
|
||||||
if (NB_BUILD_TESTS)
|
|
||||||
add_subdirectory(./tests )
|
|
||||||
endif()
|
|
||||||
@ -1,2 +0,0 @@
|
|||||||
@PACKAGE_INIT@
|
|
||||||
include("${CMAKE_CURRENT_LIST_DIR}/NBDataTargets.cmake")
|
|
||||||
@ -1,252 +0,0 @@
|
|||||||
#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
|
|
||||||
@ -1,6 +0,0 @@
|
|||||||
#include <NBCore/Errors.hpp>
|
|
||||||
#include <NBData/FFT.hpp>
|
|
||||||
|
|
||||||
namespace nb {
|
|
||||||
|
|
||||||
} // namespace nb
|
|
||||||
@ -1,22 +0,0 @@
|
|||||||
cmake_minimum_required(VERSION 3.26.0)
|
|
||||||
|
|
||||||
if (NB_BUILD_TESTS)
|
|
||||||
enable_testing()
|
|
||||||
include(GoogleTest)
|
|
||||||
|
|
||||||
add_executable(TestData
|
|
||||||
./testFFT.cpp
|
|
||||||
)
|
|
||||||
target_link_libraries(TestData
|
|
||||||
NBData
|
|
||||||
GTest::gtest_main
|
|
||||||
)
|
|
||||||
gtest_discover_tests(TestData)
|
|
||||||
|
|
||||||
add_executable(FFTSpeedTest
|
|
||||||
./fftspeedtest.cpp
|
|
||||||
)
|
|
||||||
target_link_libraries(FFTSpeedTest
|
|
||||||
NBData
|
|
||||||
)
|
|
||||||
endif()
|
|
||||||
@ -1,5 +0,0 @@
|
|||||||
#include <NBData/FFT.hpp>
|
|
||||||
|
|
||||||
int main() {
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
@ -1,27 +0,0 @@
|
|||||||
#include <gtest/gtest.h>
|
|
||||||
|
|
||||||
#include <NBData/ComplexNumbers.hpp>
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_EQ(const nb::Complex<A>& a, const nb::Complex<B>& b) {
|
|
||||||
ASSERT_EQ(a.r, b.r);
|
|
||||||
ASSERT_EQ(a.i, b.i);
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_ALMOST_EQ(const nb::Complex<A>& a, const nb::Complex<B>& b) {
|
|
||||||
ASSERT_FLOAT_EQ(a.r, b.r);
|
|
||||||
ASSERT_FLOAT_EQ(a.i, b.i);
|
|
||||||
}
|
|
||||||
|
|
||||||
TEST(ComplexNumberTests, TestConstruction) {
|
|
||||||
COMPLEX_ALMOST_EQ(nb::Complex(), nb::Complex(0,0));
|
|
||||||
COMPLEX_ALMOST_EQ(nb::Complex(5.0), nb::Complex(5, 0));
|
|
||||||
COMPLEX_ALMOST_EQ(1.5_j, nb::Complex(0.0, 1.5));
|
|
||||||
}
|
|
||||||
|
|
||||||
TEST(ComplexNumberTests, TestAdd) {
|
|
||||||
auto x = nb::Complex(1.0) - nb::Complex(2.0, 1.0);
|
|
||||||
ASSERT_FLOAT_EQ(x.r, -1);
|
|
||||||
ASSERT_FLOAT_EQ(x.i, -1);
|
|
||||||
}
|
|
||||||
@ -1,68 +0,0 @@
|
|||||||
#include <gtest/gtest.h>
|
|
||||||
|
|
||||||
#include <NBData/FFT.hpp>
|
|
||||||
|
|
||||||
using namespace nb;
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_EQ(const std::complex<A>& a, const std::complex<B>& b) {
|
|
||||||
ASSERT_EQ(a.real(), b.real());
|
|
||||||
ASSERT_EQ(a.imag(), b.imag());
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_FLOAT_EQ(
|
|
||||||
const std::complex<A>& a,
|
|
||||||
const std::complex<B>& b,
|
|
||||||
const long double epsilon=0
|
|
||||||
) {
|
|
||||||
if (epsilon) {
|
|
||||||
ASSERT_NEAR(a.real(), b.real(), epsilon);
|
|
||||||
ASSERT_NEAR(a.imag(), b.imag(), epsilon);
|
|
||||||
} else {
|
|
||||||
ASSERT_FLOAT_EQ(a.real(), b.real());
|
|
||||||
ASSERT_FLOAT_EQ(a.imag(), b.imag());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_EQ_VEC(
|
|
||||||
const std::vector<std::complex<A>>& a,
|
|
||||||
const std::vector<std::complex<A>>& b
|
|
||||||
) {
|
|
||||||
if (a.size() != b.size()) { FAIL(); }
|
|
||||||
for (int i = 0; i < a.size(); ++i) {
|
|
||||||
COMPLEX_EQ(a[i], b[i]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename A, typename B>
|
|
||||||
void COMPLEX_FLOAT_EQ_VEC(
|
|
||||||
const std::vector<std::complex<A>>& a,
|
|
||||||
const std::vector<std::complex<B>>& b,
|
|
||||||
const long double epsilon
|
|
||||||
) {
|
|
||||||
if (a.size() != b.size()) { FAIL(); }
|
|
||||||
for (int i = 0; i < a.size(); ++i) {
|
|
||||||
COMPLEX_FLOAT_EQ(a[i], b[i], epsilon);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
TEST(FFTTest, TestFFT) {
|
|
||||||
std::vector<std::vector<std::complex<double>>> input = {
|
|
||||||
{1.0, 0.5, 0},
|
|
||||||
{1,1,1,1,1,1,1,1},
|
|
||||||
|
|
||||||
};
|
|
||||||
std::vector<std::vector<std::complex<double>>> expected = {
|
|
||||||
{1.5,1.f-0.5if,0.5,1.f+0.5if},
|
|
||||||
{8, 0, 0, 0, 0, 0, 0, 0}
|
|
||||||
};
|
|
||||||
for (int i = 0; i < input.size(); ++i) {
|
|
||||||
COMPLEX_FLOAT_EQ_VEC(
|
|
||||||
nb::fft2(input[i]),
|
|
||||||
expected[i],
|
|
||||||
1e-10
|
|
||||||
);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
Loading…
x
Reference in New Issue
Block a user