-
Notifications
You must be signed in to change notification settings - Fork 13
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Itarative radix-2 FFT algorithm implementation
* using but-reversal permutation. * FFT block contains 3 template parameters: data type and algorithm type. * Add test of different FFT algorithms.
- Loading branch information
slebedev
committed
Sep 28, 2023
1 parent
b2d5c75
commit 2d1fad7
Showing
9 changed files
with
734 additions
and
429 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,112 @@ | ||
| #ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP | ||
| #define GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP | ||
|
|
||
| #include "fft_types.hpp" | ||
| #include "node.hpp" | ||
| #include "window.hpp" | ||
| #include <ranges> | ||
|
|
||
| namespace gr::algorithm { | ||
|
|
||
| template<typename T> | ||
| requires(ComplexType<T> || std::floating_point<T>) | ||
| struct FFT { | ||
| using PrecisionType = typename FFTInternal<T>::PrecisionType; | ||
| using OutDataType = typename FFTInternal<T>::OutDataType; | ||
|
|
||
| std::vector<OutDataType> twiddleFactors{}; | ||
| std::size_t fftSize{ 0 }; | ||
|
|
||
| FFT() = default; | ||
| FFT(const FFT &rhs) = delete; | ||
| FFT(FFT &&rhs) noexcept = delete; | ||
| FFT & | ||
| operator=(const FFT &rhs) | ||
| = delete; | ||
| FFT & | ||
| operator=(FFT &&rhs) noexcept | ||
| = delete; | ||
|
|
||
| ~FFT() = default; | ||
|
|
||
| void | ||
| initAll() { | ||
| precomputeTwiddleFactors(); | ||
| } | ||
|
|
||
| std::vector<OutDataType> | ||
| computeFFT(const std::vector<T> &in) noexcept { | ||
| std::vector<OutDataType> out(in.size()); | ||
| computeFFT(in, out); | ||
| return out; | ||
| } | ||
|
|
||
| void | ||
| computeFFT(const std::vector<T> &in, std::vector<OutDataType> &out) noexcept { | ||
| if (!std::has_single_bit(in.size())) { | ||
| throw std::runtime_error(fmt::format("Input data must have 2^N samples, input size: ", in.size())); | ||
| } | ||
| if (fftSize != in.size()) { | ||
| fftSize = in.size(); | ||
| initAll(); | ||
| } | ||
|
|
||
| // For the moment no optimization for real data inputs, just create complex with zero imaginary value. | ||
| if constexpr (!ComplexType<T>) { | ||
| std::transform(in.begin(), in.end(), out.begin(), [](const auto c) { return OutDataType(c, 0.); }); | ||
| } else { | ||
| std::copy(in.begin(), in.end(), out.begin()); | ||
| } | ||
|
|
||
| /** | ||
| * Real-valued fast fourier transform algorithms | ||
| * H.V. Sorensen, D.L. Jones, M.T. Heideman, C.S. Burrus (1987), | ||
| * in: IEEE Trans on Acoustics, Speech, & Signal Processing, 35 | ||
| */ | ||
| bitReversalPermutation(out); | ||
|
|
||
| std::size_t omega_kCounter = 0; | ||
| for (std::size_t s = 2; s <= fftSize; s *= 2) { | ||
| const auto half_s = s / 2; | ||
| for (std::size_t k = 0; k < fftSize; k += s) { | ||
| for (std::size_t j = 0; j < half_s; j++) { | ||
| const auto t{ twiddleFactors[omega_kCounter++] * out[k + j + half_s] }; | ||
| const auto u{ out[k + j] }; | ||
| out[k + j] = u + t; | ||
| out[k + j + half_s] = u - t; | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| private: | ||
| void | ||
| bitReversalPermutation(std::vector<OutDataType> &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); | ||
| rev ^= fftSize - (fftSize >> maskLen); | ||
| } | ||
| } | ||
|
|
||
| void | ||
| precomputeTwiddleFactors() { | ||
| twiddleFactors.clear(); | ||
| const auto minus2Pi = PrecisionType(-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))) }; | ||
| for (std::size_t k = 0; k < fftSize; k += s) { | ||
| OutDataType wk{ 1., 0. }; | ||
| for (std::size_t j = 0; j < m; j++) { | ||
| twiddleFactors.push_back(wk); | ||
| wk *= w; | ||
| } | ||
| } | ||
| } | ||
| } | ||
| }; | ||
|
|
||
| } // namespace gr::algorithm | ||
|
|
||
| #endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_HPP |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,34 @@ | ||
| #ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP | ||
| #define GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP | ||
|
|
||
| #include "fft_types.hpp" | ||
| #include <algorithm> | ||
| #include <cmath> | ||
| #include <vector> | ||
|
|
||
| namespace gr::algorithm { | ||
|
|
||
| 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())); | ||
| } | ||
| using PrecisionType = typename T::value_type; | ||
| std::transform(fftOut.begin(), std::next(fftOut.begin(), 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); | ||
| }); | ||
| } | ||
|
|
||
| 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())); | ||
| } | ||
| std::transform(fftOut.begin(), std::next(fftOut.begin(), phaseSpectrum.size()), phaseSpectrum.begin(), [](const auto &c) { return static_cast<U>(std::atan2(c.imag(), c.real())); }); | ||
| } | ||
|
|
||
| } // namespace gr::algorithm | ||
| #endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_COMMON_HPP |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,37 @@ | ||
| #ifndef GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP | ||
| #define GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP | ||
|
|
||
| #include "node.hpp" | ||
|
|
||
| namespace gr::algorithm { | ||
| template<typename T> | ||
| concept ComplexType = std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>; | ||
|
|
||
| template<typename T> | ||
| concept DoubleType = std::is_same_v<T, std::complex<double>> || std::is_same_v<T, double>; | ||
|
|
||
| // clang-format off | ||
| template <typename signedT> struct MakeSigned { using type = signedT;}; | ||
| template <std::integral signedT> struct MakeSigned<signedT> { using type = std::make_signed_t<signedT>; }; | ||
| template<typename evalU> struct EvalOutputType { using type1 = evalU; using type2 = typename MakeSigned<evalU>::type;}; | ||
| template<ComplexType evalU> struct EvalOutputType<evalU> { using type1 = typename evalU::value_type; using type2 = evalU;}; | ||
|
|
||
| template<typename T> | ||
| struct FFTInternal { | ||
| using PrecisionType = float; | ||
| using InDataType = std::conditional_t<ComplexType<T>, std::complex<float>, float>; | ||
| using OutDataType = std::complex<float>; | ||
| }; | ||
|
|
||
| template<DoubleType T> | ||
| struct FFTInternal<T> { | ||
| using PrecisionType = double; | ||
| using InDataType = std::conditional_t<ComplexType<T>, std::complex<double>, double>; | ||
| using OutDataType = std::complex<double>; | ||
| }; | ||
|
|
||
| template<typename T> | ||
| using FFTOutputType = std::conditional_t<ComplexType<T>, typename EvalOutputType<T>::type1, typename EvalOutputType<T>::type2>; | ||
|
|
||
| } // namespace gr::algorithm | ||
| #endif // GRAPH_PROTOTYPE_ALGORITHM_FFT_TYPES_HPP |
Oops, something went wrong.