Skip to content

Commit

Permalink
WIP: refactoring of FFT.
Browse files Browse the repository at this point in the history
  • Loading branch information
@drslebedev
drslebedev committed Oct 11, 2023
1 parent a4da95b commit d4ebfd5
Show file tree
Hide file tree
Showing 7 changed files with 300 additions and 237 deletions.
33 changes: 16 additions & 17 deletions algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,13 @@

namespace gr::algorithm {

template<typename T>
requires(ComplexType<T> || std::floating_point<T>)
template<typename TInput, typename TOutput>
requires((ComplexType<TInput> || std::floating_point<TInput>) && (ComplexType<TOutput>) )
struct FFT {
using PrecisionType = FFTAlgoPrecision<T>::type;
using OutDataType = std::conditional_t<ComplexType<T>, T, std::complex<T>>;
using Precision = TOutput::value_type;

std::vector<OutDataType> twiddleFactors{};
std::size_t fftSize{ 0 };
std::vector<TOutput> twiddleFactors{};
std::size_t fftSize{ 0 };

FFT() = default;
FFT(const FFT &rhs) = delete;
Expand All @@ -34,15 +33,15 @@ struct FFT {
precomputeTwiddleFactors();
}

std::vector<OutDataType>
computeFFT(const std::vector<T> &in) {
std::vector<OutDataType> out(in.size());
computeFFT(in, out);
auto
compute(const std::ranges::input_range auto &in) {
std::vector<TOutput> out(in.size());
compute(in, out);
return out;
}

void
computeFFT(const std::vector<T> &in, std::vector<OutDataType> &out) {
compute(const std::ranges::input_range auto &in, std::ranges::output_range<TOutput> auto &out) {
if (!std::has_single_bit(in.size())) {
throw std::invalid_argument(fmt::format("Input data must have 2^N samples, input size: ", in.size()));
}
Expand All @@ -52,8 +51,8 @@ struct FFT {
}

// For the moment no optimization for real data inputs, just create complex with zero imaginary value.
if constexpr (!ComplexType<T>) {
std::ranges::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return OutDataType(c, 0.); });
if constexpr (!ComplexType<TInput>) {
std::ranges::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return TOutput(c, 0.); });
} else {
std::ranges::copy(in.begin(), in.end(), out.begin());
}
Expand Down Expand Up @@ -81,7 +80,7 @@ struct FFT {

private:
void
bitReversalPermutation(std::vector<OutDataType> &vec) const noexcept {
bitReversalPermutation(std::vector<TOutput> &vec) const noexcept {
for (std::size_t j = 0, rev = 0; j < fftSize; j++) {
if (j < rev) std::swap(vec[j], vec[rev]);
auto maskLen = static_cast<std::size_t>(std::countr_zero(j + 1) + 1);
Expand All @@ -92,12 +91,12 @@ struct FFT {
void
precomputeTwiddleFactors() {
twiddleFactors.clear();
const auto minus2Pi = PrecisionType(-2. * std::numbers::pi);
const auto minus2Pi = Precision(-2. * std::numbers::pi);
for (std::size_t s = 2; s <= fftSize; s *= 2) {
const std::size_t m{ s / 2 };
const OutDataType w{ std::exp(OutDataType(0., minus2Pi / static_cast<PrecisionType>(s))) };
const TOutput w{ std::exp(TOutput(0., minus2Pi / static_cast<Precision>(s))) };
for (std::size_t k = 0; k < fftSize; k += s) {
OutDataType wk{ 1., 0. };
TOutput wk{ 1., 0. };
for (std::size_t j = 0; j < m; j++) {
twiddleFactors.push_back(wk);
wk *= w;
Expand Down
75 changes: 58 additions & 17 deletions algorithm/include/gnuradio-4.0/algorithm/fourier/fft_common.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,33 +6,74 @@
#include <vector>

#include <fmt/format.h>
#include <ranges>

#include "fft_types.hpp"

namespace gr::algorithm {
namespace gr::algorithm { namespace fft {

template<ComplexType T, typename U>
void
computeMagnitudeSpectrum(const std::vector<T> &fftOut, std::vector<U> &magnitudeSpectrum, std::size_t fftSize, bool outputInDb) {
if (fftOut.size() < magnitudeSpectrum.size()) {
throw std::invalid_argument(fmt::format("FFT vector size ({}) must be more or equal than magnitude vector size ({}).", fftOut.size(), magnitudeSpectrum.size()));
struct ConfigMagnitude {
bool computeHalfSpectrum = false;
bool outputInDb = false;
};

template<std::ranges::input_range TContainerIn, std::ranges::output_range<typename TContainerIn::value_type::value_type> TContainerOut = std::vector<typename TContainerIn::value_type::value_type>,
typename T = TContainerIn::value_type>
requires(std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>)
auto
computeMagnitudeSpectrum(const TContainerIn &fftIn, TContainerOut &&magOut = {}, ConfigMagnitude config = {}) {
std::size_t magSize = config.computeHalfSpectrum ? fftIn.size() : (fftIn.size() / 2);
if constexpr (requires(std::size_t n) { magOut.resize(n); }) {
if (magOut.size() != magSize) {
magOut.resize(magSize);
}
} else {
static_assert(std::tuple_size_v<TContainerIn> == std::tuple_size_v<TContainerOut>, "Size mismatch for fixed-size container.");
}

using PrecisionType = typename T::value_type;
std::transform(fftOut.begin(), std::next(fftOut.begin(), static_cast<std::ptrdiff_t>(magnitudeSpectrum.size())), magnitudeSpectrum.begin(), [fftSize, outputInDb](const auto &c) {
const auto mag{ std::hypot(c.real(), c.imag()) * PrecisionType(2.0) / static_cast<PrecisionType>(fftSize) };
return static_cast<U>(outputInDb ? PrecisionType(20.) * std::log10(std::abs(mag)) : mag);
std::transform(fftIn.begin(), std::next(fftIn.begin(), static_cast<std::ptrdiff_t>(magSize)), magOut.begin(), [fftSize = fftIn.size(), outputInDb = config.outputInDb](const auto &c) {
const auto mag{ std::hypot(c.real(), c.imag()) * PrecisionType(2.) / static_cast<PrecisionType>(fftSize) };
return outputInDb ? PrecisionType(20.) * std::log10(std::abs(mag)) : mag;
});

return magOut;
}

template<std::ranges::input_range TContainerIn, typename T = TContainerIn::value_type>
requires(std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>)
auto
computeMagnitudeSpectrum(const TContainerIn &fftIn, ConfigMagnitude config) {
return computeMagnitudeSpectrum(fftIn, {}, config);
}

template<ComplexType T, typename U>
void
computePhaseSpectrum(const std::vector<T> &fftOut, std::vector<U> &phaseSpectrum) {
if (fftOut.size() < phaseSpectrum.size()) {
throw std::invalid_argument(fmt::format("FFT vector size ({}) must be more or equal than phase vector size ({}).", fftOut.size(), phaseSpectrum.size()));
struct ConfigPhase {
bool computeHalfSpectrum = false;
};

template<std::ranges::input_range TContainerIn, std::ranges::output_range<typename TContainerIn::value_type::value_type> TContainerOut = TContainerIn, typename T = TContainerIn::value_type>
requires(std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>)
auto
computePhaseSpectrum(const TContainerIn &fftIn, TContainerOut &&phaseOut = {}, ConfigPhase config = {}) {
std::size_t phaseSize = config.computeHalfSpectrum ? fftIn.size() : (fftIn.size() / 2);
if constexpr (requires(std::size_t n) { phaseOut.resize(n); }) {
if (phaseOut.size() != phaseSize) {
phaseOut.resize(phaseSize);
}
} else {
static_assert(std::tuple_size_v<TContainerIn> == std::tuple_size_v<TContainerOut>, "Size mismatch for fixed-size container.");
}
std::transform(fftOut.begin(), std::next(fftOut.begin(), static_cast<std::ptrdiff_t>(phaseSpectrum.size())), phaseSpectrum.begin(),
[](const auto &c) { return static_cast<U>(std::atan2(c.imag(), c.real())); });
std::transform(fftIn.begin(), std::next(fftIn.begin(), static_cast<std::ptrdiff_t>(phaseOut.size())), phaseOut.begin(), [](const auto &c) { return std::atan2(c.imag(), c.real()); });

return phaseOut;
}

template<std::ranges::input_range TContainerIn, typename T = TContainerIn::value_type>
requires(std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>)
auto
computePhaseSpectrum(const TContainerIn &fftIn, ConfigPhase config) {
return computePhaseSpectrum(fftIn, {}, config);
}

} // namespace gr::algorithm
}} // namespace gr::algorithm::fft
#endif // GNURADIO_ALGORITHM_FFT_COMMON_HPP
65 changes: 33 additions & 32 deletions algorithm/include/gnuradio-4.0/algorithm/fourier/fftw.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,18 @@ namespace gr::algorithm {
template<typename T>
concept FFTwDoubleType = std::is_same_v<T, std::complex<double>> || std::is_same_v<T, double>;

template<typename T>
requires(ComplexType<T> || std::floating_point<T>)
template<typename TInput, typename TOutput>
requires((ComplexType<TInput> || std::floating_point<TInput>) && (ComplexType<TOutput>) )
struct FFTw {
private:
inline static std::mutex fftw_plan_mutex;

public:
// clang-format off
template<typename FftwT>
struct fftwImpl {
template<typename TDataInput, typename TDataOutput>
struct FFTwImpl {
using PlanType = fftwf_plan;
using InAlgoDataType = std::conditional_t<ComplexType<FftwT>, fftwf_complex, float>;
using InAlgoDataType = std::conditional_t<ComplexType<TDataInput>, fftwf_complex, float>;
using OutAlgoDataType = fftwf_complex;
using InUniquePtr = std::unique_ptr<InAlgoDataType [], decltype([](InAlgoDataType *ptr) { fftwf_free(ptr); })>;
using OutUniquePtr = std::unique_ptr<OutAlgoDataType [], decltype([](OutAlgoDataType *ptr) { fftwf_free(ptr); })>;
Expand Down Expand Up @@ -51,10 +51,11 @@ struct FFTw {
static void forgetWisdom() {fftwf_forget_wisdom();}
};

template<FFTwDoubleType FftwDoubleT>
struct fftwImpl<FftwDoubleT> {
template<typename TDataInput, FFTwDoubleType TDataOutput>
// requires (std::is_same_v<TDataOutput, std::complex<double>> || std::is_same_v<TDataOutput, double>)
struct FFTwImpl<TDataInput, TDataOutput> {
using PlanType = fftw_plan;
using InAlgoDataType = std::conditional_t<ComplexType<FftwDoubleT>, fftw_complex, double>;
using InAlgoDataType = std::conditional_t<ComplexType<TDataInput>, fftw_complex, double>;
using OutAlgoDataType = fftw_complex;
using InUniquePtr = std::unique_ptr<InAlgoDataType [], decltype([](InAlgoDataType *ptr) { fftwf_free(ptr); })>;
using OutUniquePtr = std::unique_ptr<OutAlgoDataType [], decltype([](OutAlgoDataType *ptr) { fftwf_free(ptr); })>;
Expand Down Expand Up @@ -84,12 +85,12 @@ struct FFTw {
};

// clang-format on
using OutDataType = std::conditional_t<ComplexType<T>, T, std::complex<T>>;
using InAlgoDataType = typename fftwImpl<T>::InAlgoDataType;
using OutAlgoDataType = typename fftwImpl<T>::OutAlgoDataType;
using InUniquePtr = typename fftwImpl<T>::InUniquePtr;
using OutUniquePtr = typename fftwImpl<T>::OutUniquePtr;
using PlanUniquePtr = typename fftwImpl<T>::PlanUniquePtr;
// Precision of the algorithm is defined by the output type `TOutput`
using InAlgoDataType = typename FFTwImpl<TInput, TOutput>::InAlgoDataType;
using OutAlgoDataType = typename FFTwImpl<TInput, TOutput>::OutAlgoDataType;
using InUniquePtr = typename FFTwImpl<TInput, TOutput>::InUniquePtr;
using OutUniquePtr = typename FFTwImpl<TInput, TOutput>::OutUniquePtr;
using PlanUniquePtr = typename FFTwImpl<TInput, TOutput>::PlanUniquePtr;

std::size_t fftSize{ 0 };
std::string wisdomPath{ ".gr_fftw_wisdom" };
Expand All @@ -111,19 +112,19 @@ struct FFTw {

~FFTw() { clearFftw(); }

std::vector<OutDataType>
computeFFT(const std::vector<T> &in) {
auto
compute(const std::ranges::input_range auto &in) {
if (fftSize != in.size()) {
fftSize = in.size();
initAll();
}
std::vector<OutDataType> out(getOutputSize());
computeFFT(in, out);
std::vector<TOutput> out(getOutputSize());
compute(in, out);
return out;
}

void
computeFFT(const std::vector<T> &in, std::vector<OutDataType> &out) {
compute(const std::ranges::input_range auto &in, std::ranges::output_range<TOutput> auto &out) {
if (!std::has_single_bit(in.size())) {
throw std::invalid_argument(fmt::format("Input data must have 2^N samples, input size: ", in.size()));
}
Expand All @@ -137,46 +138,46 @@ struct FFTw {
throw std::out_of_range(fmt::format("Output vector size ({}) is not enough, at least {} needed. ", out.size(), getOutputSize()));
}

static_assert(sizeof(InAlgoDataType) == sizeof(T), "Input std::complex<T> and T[2] must have the same size.");
static_assert(sizeof(InAlgoDataType) == sizeof(TOutput), "Input std::complex<T> and T[2] must have the same size.");
std::memcpy(fftwIn.get(), &(*in.begin()), sizeof(InAlgoDataType) * fftSize);

fftwImpl<T>::execute(fftwPlan.get());
FFTwImpl<TInput, TOutput>::execute(fftwPlan.get());

static_assert(sizeof(OutDataType) == sizeof(OutAlgoDataType), "Output std::complex<T> and T[2] must have the same size.");
std::memcpy(out.data(), fftwOut.get(), sizeof(OutDataType) * getOutputSize());
static_assert(sizeof(TOutput) == sizeof(OutAlgoDataType), "Output std::complex<T> and T[2] must have the same size.");
std::memcpy(out.data(), fftwOut.get(), sizeof(TOutput) * getOutputSize());
}

[[nodiscard]] inline int
importWisdom() const {
// lock file while importing wisdom?
return fftwImpl<T>::importWisdomFromFilename(wisdomPath);
return FFTwImpl<TInput, TOutput>::importWisdomFromFilename(wisdomPath);
}

[[nodiscard]] inline int
exportWisdom() const {
// lock file while exporting wisdom?
return fftwImpl<T>::exportWisdomToFilename(wisdomPath);
return FFTwImpl<TInput, TOutput>::exportWisdomToFilename(wisdomPath);
}

[[nodiscard]] inline int
importWisdomFromString(const std::string &wisdomString) const {
return fftwImpl<T>::importWisdomFromString(wisdomString);
return FFTwImpl<TInput, TOutput>::importWisdomFromString(wisdomString);
}

[[nodiscard]] std::string
exportWisdomToString() const {
return fftwImpl<T>::exportWisdomToString();
return FFTwImpl<TInput, TOutput>::exportWisdomToString();
}

inline void
forgetWisdom() const {
return fftwImpl<T>::forgetWisdom();
return FFTwImpl<TInput, TOutput>::forgetWisdom();
}

private:
[[nodiscard]] constexpr std::size_t
getOutputSize() const {
if constexpr (ComplexType<T>) {
if constexpr (ComplexType<TInput>) {
return fftSize;
} else {
return 1 + fftSize / 2;
Expand All @@ -186,14 +187,14 @@ struct FFTw {
void
initAll() {
clearFftw();
fftwIn = InUniquePtr(static_cast<InAlgoDataType *>(fftwImpl<T>::malloc(sizeof(InAlgoDataType) * fftSize)));
fftwOut = OutUniquePtr(static_cast<OutAlgoDataType *>(fftwImpl<T>::malloc(sizeof(OutAlgoDataType) * getOutputSize())));
fftwIn = InUniquePtr(static_cast<InAlgoDataType *>(FFTwImpl<TInput, TOutput>::malloc(sizeof(InAlgoDataType) * fftSize)));
fftwOut = OutUniquePtr(static_cast<OutAlgoDataType *>(FFTwImpl<TInput, TOutput>::malloc(sizeof(OutAlgoDataType) * getOutputSize())));

{
std::lock_guard lg{ fftw_plan_mutex };
// what to do if error is returned
std::ignore = importWisdom();
fftwPlan = PlanUniquePtr(fftwImpl<T>::plan(static_cast<int>(fftSize), fftwIn.get(), fftwOut.get(), sign, flags));
fftwPlan = PlanUniquePtr(FFTwImpl<TInput, TOutput>::plan(static_cast<int>(fftSize), fftwIn.get(), fftwOut.get(), sign, flags));
std::ignore = exportWisdom();
}
}
Expand Down
Loading

0 comments on commit d4ebfd5

Please sign in to comment.