Skip to content

Commit

Permalink
FFT block and algorithms refactoring:
Browse files Browse the repository at this point in the history 
* Change signature of the FFT algorithm classs. Add TInput and TOutput template parameters
* Precision of the calculation is defined by the TOutput precision
* Add options to calculate phase in deg and unwrapped
* Clean up FFT unittests
  • Loading branch information
@drslebedev
drslebedev committed Oct 17, 2023
1 parent f491f52 commit bf360a2
Show file tree
Hide file tree
Showing 9 changed files with 447 additions and 320 deletions.
49 changes: 28 additions & 21 deletions algorithm/include/gnuradio-4.0/algorithm/fourier/fft.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,17 @@

#include <ranges>

#include "fft_types.hpp"
#include "window.hpp"

namespace gr::algorithm {

template<typename T>
requires(ComplexType<T> || std::floating_point<T>)
template<typename TInput, typename TOutput = std::conditional<gr::meta::complex_like<TInput>, TInput, std::complex<typename TInput::value_type>>>
requires((gr::meta::complex_like<TInput> || std::floating_point<TInput>) && (gr::meta::complex_like<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 +32,16 @@ struct FFT {
precomputeTwiddleFactors();
}

std::vector<OutDataType>
computeFFT(const std::vector<T> &in) {
std::vector<OutDataType> out(in.size());
computeFFT(in, out);
return out;
}
auto
compute(const std::ranges::input_range auto &in, std::ranges::output_range<TOutput> auto &&out) {
if constexpr (requires(std::size_t n) { out.resize(n); }) {
if (out.size() != in.size()) {
out.resize(in.size());
}
} else {
static_assert(std::tuple_size_v<decltype(in)> == std::tuple_size_v<decltype(out)>, "Size mismatch for fixed-size container.");
}

void
computeFFT(const std::vector<T> &in, std::vector<OutDataType> &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,9 +51,10 @@ 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 (!gr::meta::complex_like<TInput>) {
std::ranges::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return TOutput(c, 0.); });
} else {
// precision is defined by output type, if cast is needed, let `std::copy` do the job
std::ranges::copy(in.begin(), in.end(), out.begin());
}

Expand All @@ -77,11 +77,18 @@ struct FFT {
}
}
}

return out;
}

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

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 +99,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
107 changes: 90 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 @@ -3,36 +3,109 @@

#include <algorithm>
#include <cmath>
#include <numbers>
#include <vector>

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

#include "fft_types.hpp"
namespace gr::algorithm::fft {

namespace gr::algorithm {
struct ConfigMagnitude {
bool computeHalfSpectrum = false;
bool outputInDb = false;
};

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()));
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>
struct ConfigPhase {
bool computeHalfSpectrum = false;
bool outputInDeg = false;
bool unwrapPhase = false;
};

template<std::ranges::input_range TContainerInOut, typename T = TContainerInOut::value_type>
requires(std::floating_point<T>)
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()));
unwrapPhase(TContainerInOut &phase) {
const auto pi = std::numbers::pi_v<T>;
auto prev = phase.front();
std::transform(phase.begin() + 1, phase.end(), phase.begin() + 1, [&prev, pi](T &current) {
T diff = current - prev;
while (diff > pi) {
current -= 2 * pi;
diff = current - prev;
}
while (diff < -pi) {
current += 2 * pi;
diff = current - prev;
}
prev = current;
return current;
});
}

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
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(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()); });

if (config.unwrapPhase) {
unwrapPhase(phaseOut);
}

if (config.outputInDeg) {
std::transform(phaseOut.begin(), phaseOut.end(), phaseOut.begin(),
[](const auto &phase) { return phase * static_cast<typename T::value_type>(180.) * std::numbers::inv_pi_v<typename T::value_type>; });
}
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())); });

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
28 changes: 0 additions & 28 deletions algorithm/include/gnuradio-4.0/algorithm/fourier/fft_types.hpp
View file

This file was deleted.

Loading

0 comments on commit bf360a2

Please sign in to comment.