commit ea3a39f55050392fa84eb0eed50bbe708693edf9 Author: Andreas Tsouchlos Date: Wed Feb 18 21:32:22 2026 +0100 Initial commit diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..806794c --- /dev/null +++ b/.clang-format @@ -0,0 +1,28 @@ +BasedOnStyle: LLVM +Language: Cpp + +IndentWidth: 4 +UseTab: Never +#NamespaceIndentation: All + +PointerAlignment: Left +AccessModifierOffset: -4 +AlwaysBreakTemplateDeclarations: true +LambdaBodyIndentation: Signature + +MaxEmptyLinesToKeep: 3 +#ColumnLimit: 128 + +CompactNamespaces: true +FixNamespaceComments: true + +AllowShortFunctionsOnASingleLine: false +AllowShortIfStatementsOnASingleLine: true + +AlignConsecutiveAssignments: true +AlignConsecutiveBitFields: true +AlignConsecutiveDeclarations: true +AlignConsecutiveMacros: true + +#BraceWrapping: +# BeforeElse: true diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..a867038 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,11 @@ +cmake_minimum_required(VERSION 3.26) +project(R_S) + +set(CMAKE_CXX_STANDARD 23) + +add_executable(recursive_fft + fft/recursive_fft.cpp) +add_executable(recursive_fft_cplx + fft/recursive_fft_cplx.cpp) +add_executable(iterative_fft_cplx + fft/iterative_fft_cplx.cpp) diff --git a/fft/.idea/.gitignore b/fft/.idea/.gitignore new file mode 100644 index 0000000..26d3352 --- /dev/null +++ b/fft/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/fft/.idea/fft.iml b/fft/.idea/fft.iml new file mode 100644 index 0000000..8b8c395 --- /dev/null +++ b/fft/.idea/fft.iml @@ -0,0 +1,12 @@ + + + + + + + + + + \ No newline at end of file diff --git a/fft/.idea/inspectionProfiles/Project_Default.xml b/fft/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 0000000..700f76d --- /dev/null +++ b/fft/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,23 @@ + + + + \ No newline at end of file diff --git a/fft/.idea/inspectionProfiles/profiles_settings.xml b/fft/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/fft/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/fft/.idea/misc.xml b/fft/.idea/misc.xml new file mode 100644 index 0000000..e3fd5f7 --- /dev/null +++ b/fft/.idea/misc.xml @@ -0,0 +1,4 @@ + + + + \ No newline at end of file diff --git a/fft/.idea/modules.xml b/fft/.idea/modules.xml new file mode 100644 index 0000000..2e6856f --- /dev/null +++ b/fft/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/fft/fft.py b/fft/fft.py new file mode 100644 index 0000000..ed9da34 --- /dev/null +++ b/fft/fft.py @@ -0,0 +1,31 @@ +import numpy as np + + +def fft(x: np.array) -> np.array: + N = len(x) + + if N == 1: + return x + else: + even_fft = fft(x[::2]) + odd_fft = fft(x[1::2]) + + factors = np.exp(-2j * np.pi * np.arange(N) / N) + + result = np.concatenate([even_fft + factors[:N // 2] * odd_fft, + even_fft + factors[N // 2:] * odd_fft]) + + return result + + +def main(): + x = np.arange(8) + fft(x) + # # x = np.random.normal(size=4) + # x = np.array([1, 0, 1, 0, 0, 0, 0, 0]) + # print(f"own:\t{fft(x)}") + # print(f"numpy:\t{np.fft.fft(x)}") + + +if __name__ == "__main__": + main() diff --git a/fft/iterative_fft.py b/fft/iterative_fft.py new file mode 100644 index 0000000..7829681 --- /dev/null +++ b/fft/iterative_fft.py @@ -0,0 +1,48 @@ +import numpy as np + + +def bit_reverse(val: int, num_bits: int) -> None: + b = '{:0{num_bits}b}'.format(val, num_bits=num_bits) + return int(b[::-1], 2) + + +def fft_reorder(x: np.array) -> np.array: + result = np.zeros(x.size) + for i in range(x.size): + result[i] = x[bit_reverse(i, int(np.log2(x.size)))] + return result + + +def fft(x: np.array) -> np.array: + x = fft_reorder(x).astype('complex64') + + N0 = x.size + N = 2 + # while N <= N0: + for N in [2, 4, 8]: + first_sub_fft_index = 0 + while first_sub_fft_index < N0: + for k in range(N // 2): + factor = np.exp(-2j * np.pi * k / N) + + temp = factor * x[first_sub_fft_index + N // 2 + k] + + x[first_sub_fft_index + N//2 + k]\ + = x[first_sub_fft_index + k] - temp + x[first_sub_fft_index + k]\ + = x[first_sub_fft_index + k] + temp + + first_sub_fft_index += N + # N *= 2 + + return x + + +def main(): + x = np.arange(8) + print(fft(x)) + print(np.fft.fft(np.arange(8))) + + +if __name__ == "__main__": + main() diff --git a/fft/iterative_fft_cplx.cpp b/fft/iterative_fft_cplx.cpp new file mode 100644 index 0000000..47a2a28 --- /dev/null +++ b/fft/iterative_fft_cplx.cpp @@ -0,0 +1,212 @@ +#include +#include +#include +#include +#include +#include +#include +#include + + +/* + * + * Helper classes and functions + * + */ + + +template +class Complex { +public: + constexpr Complex() noexcept : mRe(), mIm() { + } + + constexpr Complex(const value_t& re, const value_t& im) noexcept + : mRe(re), mIm(im) { + } + + constexpr Complex operator-() const { + return {-mRe, -mIm}; + } + + constexpr Complex operator+(const Complex& other) const { + return {mRe + other.mRe, mIm + other.mIm}; + } + + constexpr Complex operator-(const Complex& other) const { + return *this + (-other); + } + + constexpr Complex operator*(const Complex& other) const { + return {mRe * other.mRe - mIm * other.mIm, + mIm * other.mRe + mRe * other.mIm}; + } + + operator std::string() { + return std::format("{:.2f} + j*{:.2f}", mRe, mIm); + } + + constexpr value_t& real() { + return mRe; + } + + constexpr value_t& imag() { + return mIm; + } + +private: + value_t mRe; + value_t mIm; +}; + +using ComplexFloat = Complex; + + +template +constexpr void print(std::vector> vec) { + for (auto& elem : vec) { + std::cout << "\t" << elem.real() << "; " << elem.imag() << std::endl; + } +} + + +/* + * + * FFT implementation + * + */ + + +template +constexpr value_t bit_reverse(value_t n, std::size_t b = sizeof(value_t) * 8) { + assert(b <= std::numeric_limits::digits); + + value_t rv = 0; + + for (size_t i = 0; i < b; ++i, n >>= 1) { + rv = (rv << 1) | (n & 0x01); + } + + return rv; +} + +constexpr std::vector +fft_reorder(const std::vector& x) { + std::vector result; + + for (int n = 0; n < x.size(); ++n) { + const unsigned rev_n = bit_reverse(n, log2(x.size())); + result.push_back(x[rev_n]); + } + + return result; +} + +//// TODO: Precompute Factors +//// TODO: Avoid unnecessary initialization cost of second buffer +//// TODO: Stream- or Vector-Hardware +// constexpr std::vector fft(const std::vector& x) { +// const unsigned N0 = x.size(); +// +// std::array, 2> buffers; +// buffers[0] = fft_reorder(x); +// buffers[1].resize(N0); +// +// for (int i_iter = 1; i_iter <= log2(N0); ++i_iter) { +// std::vector& input = buffers[(i_iter - 1) % 2]; +// std::vector& output = buffers[i_iter % 2]; +// +// const unsigned N = std::pow(2, i_iter); +// const unsigned num_sub_ffts = N0 / N; +// +// for (int i_sub_fft = 0; i_sub_fft < num_sub_ffts; ++i_sub_fft) { +// for (int k = 0; k < N; ++k) { +// float factor_re = std::cos((2 * std::numbers::pi * k) / N); +// float factor_im = -std::sin((2 * std::numbers::pi * k) / N); +// ComplexFloat factor = {factor_re, factor_im}; +// +// const unsigned index1 = i_sub_fft * N + k % (N / 2); +// const unsigned index2 = i_sub_fft * N + N / 2 + k % (N / 2); +// +// output[i_sub_fft * N + k] = +// input[index1] + factor * input[index2]; +// } +// } +// } +// +// return buffers[static_cast(log2(N0)) % 2]; +// } + + +std::vector precompute_factors(unsigned N0) { + std::vector result; + + for (int N = 2; N <= N0; N *= 2) { + for (int k = 0; k < N / 2; k++) { + float factor_re = std::cos((2 * std::numbers::pi * k) / N); + float factor_im = -std::sin((2 * std::numbers::pi * k) / N); + + result.push_back({factor_re, factor_im}); + } + } + + return result; +} + +constexpr unsigned get_factor_index(unsigned N, unsigned k) { + return (N - 2) / 2 + k; +} + + +const std::vector factors = precompute_factors(8); + + +// TODO: Stream- or Vector-Hardware +std::vector fft(const std::vector& input) { + const unsigned N0 = input.size(); + + std::vector x = fft_reorder(input); + + for (int N = 2; N <= N0; N *= 2) { + for (int i_sub_fft = 0; i_sub_fft < N0; i_sub_fft += N) { + for (int k = 0; k < N / 2; k++) { + const ComplexFloat& factor = factors[get_factor_index(N, k)]; + + ComplexFloat z = factor * x[i_sub_fft + k + N / 2]; + x[i_sub_fft + k + N / 2] = x[i_sub_fft + k] - z; + x[i_sub_fft + k] = x[i_sub_fft + k] + z; + } + } + } + + return x; +} + + +/* + * + * Usage + * + */ + + +int main() { + std::vector x = {{1, 0}, {2, 0}, {3, 0}, {4, 0}, + {5, 0}, {6, 0}, {7, 0}, {8, 0}}; + // std::vector x = {{1, 0}, {2, 0}, {3, 0}, {4, 0}}; + + const auto& X = fft(x); + + // int N0 = 8; + // for (int N = 2; N <= N0; N *= 2) { + // for (int k = 0; k < N / 2; k++) { + // std::cout << get_factor_index(N, k) << std::endl; + // } + // } + + std::cout << "================" << std::endl; + std::cout << "result: " << std::endl; + print(X); + + return 0; +} diff --git a/fft/real_fft_test.py b/fft/real_fft_test.py new file mode 100644 index 0000000..66ddb57 --- /dev/null +++ b/fft/real_fft_test.py @@ -0,0 +1,83 @@ +"""This script is used to visualize the behavior of the fft if the input is +purely real.""" + +import numpy as np +import seaborn as sns +import matplotlib.pyplot as plt + + +def bit_reverse(val: int, num_bits: int) -> None: + b = '{:0{num_bits}b}'.format(val, num_bits=num_bits) + return int(b[::-1], 2) + + +def fft_reorder(x: np.array) -> np.array: + result = np.zeros(x.size) + for i in range(x.size): + result[i] = x[bit_reverse(i, int(np.log2(x.size)))] + return result + + +def fft(x: np.array) -> np.array: + x = fft_reorder(x).astype('complex64') + + N0 = x.size + + # fig, axes = plt.subplots(2, int(np.log2(N0))) + + N = 2 + while N <= N0: + first_subfft_index = 0 + while first_subfft_index < N0: + temp = x[first_subfft_index + N // 2] + x[first_subfft_index + N // 2] = x[first_subfft_index] - temp + x[first_subfft_index] = x[first_subfft_index] + temp + + for k in range(1, N // 2): + factor = np.exp(-2j * np.pi * k / N) + + x[first_subfft_index + k] \ + = x[first_subfft_index + k] \ + + factor * x[first_subfft_index + k + N // 2] + + for k in range(1, N // 2): + x[first_subfft_index + N - k] \ + = x[first_subfft_index + k].conj() + # \ + # = x[first_subfft_index + k] + temp + # for k in range(N // 2): + # factor = np.exp(- 2j * np.pi * k / N) + # + # temp = factor * x[first_subfft_index + k + N // 2] + # x[first_subfft_index + k + N // 2] \ + # = x[first_subfft_index + k] - temp + # x[first_subfft_index + k] \ + # = x[first_subfft_index + k] + temp + + first_subfft_index += N + # axes[0][int(np.log2(N)) - 1].plot(x.real[1:]) + # axes[1][int(np.log2(N)) - 1].plot(x.imag[1:]) + N *= 2 + + plt.show() + + return x + + +def main(): + sns.set_theme() + + x = np.arange(16) + 1 + X = fft(x) + X_reference = np.fft.fft(x) + + for X_i, X_ref_i in zip(X, X_reference): + print( + f"{X_i.real:06.2f} + j({X_i.imag:06.2f})\t\t{X_ref_i.real:06.2f} " + f"+ j({X_ref_i.imag:06.2f})") + + print(np.allclose(X, X_reference, rtol=1e-05, atol=1e-08)) + + +if __name__ == "__main__": + main() diff --git a/fft/recursive_fft.cpp b/fft/recursive_fft.cpp new file mode 100644 index 0000000..6448bf4 --- /dev/null +++ b/fft/recursive_fft.cpp @@ -0,0 +1,85 @@ +#include +#include +#include +#include +#include +#include +#include + + +constexpr void print(auto& x_re, auto& x_im) { + for (std::tuple elem : + std::views::zip(x_re, x_im)) { + const auto re = std::get<0>(elem); + const auto im = std::get<1>(elem); + + std::cout << "\t" << re << "; " << im << std::endl; + } +} + + +constexpr std::pair, std::vector> +fft(const std::vector& x_re, const std::vector& x_im) { + // TODO: Make sure this is ignored in release mode + assert(x_re.size() == x_im.size()); + + const unsigned N = x_re.size(); + + if (N == 1) { + return {x_re, x_im}; + } else { + const auto& even_x_re = x_re | std::views::stride(2); + const auto& even_x_im = x_im | std::views::stride(2); + const auto& odd_x_re = + x_re | std::views::drop(1) | std::views::stride(2); + const auto& odd_x_im = + x_im | std::views::drop(1) | std::views::stride(2); + + const auto& [even_fft_re, even_fft_im] = + fft(std::vector(even_x_re.begin(), even_x_re.end()), + std::vector(even_x_im.begin(), even_x_im.end())); + const auto& [odd_fft_re, odd_fft_im] = + fft(std::vector(odd_x_re.begin(), odd_x_re.end()), + std::vector(odd_x_im.begin(), odd_x_im.end())); + + /// Combine even and odd ffts + + std::vector factors_re(N); + std::vector factors_im(N); + + for (int k = 0; k < N; ++k) { + factors_re[k] = std::cos((2 * std::numbers::pi * k) / N); + factors_im[k] = -std::sin((2 * std::numbers::pi * k) / N); + } + + std::vector result_re(N); + std::vector result_im(N); + + for (int k = 0; k < N / 2; ++k) { + result_re[k] = even_fft_re[k] + factors_re[k] * odd_fft_re[k] - + factors_im[k] * odd_fft_im[k]; + result_im[k] = even_fft_im[k] + factors_im[k] * odd_fft_re[k] + + factors_re[k] * odd_fft_im[k]; + } + + for (int k = 0; k < N / 2; ++k) { + result_re[N / 2 + k] = even_fft_re[k] + + factors_re[N / 2 + k] * odd_fft_re[k] - + factors_im[k] * odd_fft_im[k]; + result_im[N / 2 + k] = even_fft_im[k] + + factors_im[N / 2 + k] * odd_fft_re[k] + + factors_re[k] * odd_fft_im[k]; + } + + return {result_re, result_im}; + } +} + + +int main() { + std::vector x_re = {1, 2, 3, 4}; + std::vector x_im = {0, 0, 0, 0}; + const auto [X_re, X_im] = fft(x_re, x_im); + + return 0; +} diff --git a/fft/recursive_fft_cplx.cpp b/fft/recursive_fft_cplx.cpp new file mode 100644 index 0000000..3bcc68a --- /dev/null +++ b/fft/recursive_fft_cplx.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include +#include + + +/* + * + * Helper classes and functions + * + */ + + +template +class Complex { +public: + constexpr Complex() noexcept : mRe(), mIm() { + } + + constexpr Complex(const value_t& re, const value_t& im) noexcept + : mRe(re), mIm(im) { + } + + constexpr Complex operator-() const { + return {-mRe, -mIm}; + } + + constexpr Complex operator+(const Complex& other) const { + return {mRe + other.mRe, mIm + other.mIm}; + } + + constexpr Complex operator-(const Complex& other) const { + return *this + (-other); + } + + constexpr Complex operator*(const Complex& other) const { + return {mRe * other.mRe - mIm * other.mIm, + mIm * other.mRe + mRe * other.mIm}; + } + + constexpr value_t& real() { + return mRe; + } + + constexpr value_t& imag() { + return mIm; + } + +private: + value_t mRe; + value_t mIm; +}; + +using ComplexFloat = Complex; + + +template +constexpr void print(std::vector> vec) { + for (auto& elem : vec) { + std::cout << "\t" << elem.real() << "; " << elem.imag() << std::endl; + } +} + + +/* + * + * FFT implementation + * + */ + + +constexpr std::vector fft(const std::vector& x) { + const unsigned N = x.size(); + + if (N == 1) { + return x; + } else { + /// Compute even and odd ffts + + const auto& even_x = x | std::views::stride(2); + const auto& odd_x = x | std::views::drop(1) | std::views::stride(2); + + const auto& even_fft = + fft(std::vector(even_x.begin(), even_x.end())); + const auto& odd_fft = + fft(std::vector(odd_x.begin(), odd_x.end())); + + /// Combine even and odd ffts + + std::vector factors(N); + + for (int k = 0; k < N; ++k) { + factors[k].real() = std::cos((2 * std::numbers::pi * k) / N); + factors[k].imag() = -std::sin((2 * std::numbers::pi * k) / N); + } + + std::vector result(N); + + for (int k = 0; k < N; ++k) { + result[k] = + even_fft[k % (N / 2)] + factors[k] * odd_fft[k % (N / 2)]; + } + + return result; + } +} + + +/* + * + * Usage + * + */ + + +int main() { + std::vector x = {{1, 0}, {2, 0}, {3, 0}, {4, 0}}; + + const auto& X = fft(x); + + std::cout << "================" << std::endl; + std::cout << "result: " << std::endl; + print(X); + + return 0; +} diff --git a/fft/vector_fft.py b/fft/vector_fft.py new file mode 100644 index 0000000..2258ec9 --- /dev/null +++ b/fft/vector_fft.py @@ -0,0 +1,222 @@ +import numpy as np +import sys +import typing + +buffer = np.zeros(1024).astype('float32') + + +# TODO temp memory +class MemoryWrapper: + """Class allowing for easier access to the memory buffer.""" + + def __init__(self, N0: int, imag_offset: int): + self._current_offset = 0 + self._num_sub_ffts = 0 + self._N0 = N0 + self._imag_offset = imag_offset + + def set_kN(self, k, num_sub_ffts): + N = self._N0 / num_sub_ffts + i = bit_reverse(k, int(np.log2(N))) + self._current_offset = i * num_sub_ffts + self._num_sub_ffts = num_sub_ffts + + @property + def OddReal(self) -> np.array: + start_index = self._current_offset + self._num_sub_ffts + end_index = self._current_offset + self._num_sub_ffts * 2 + return buffer[start_index:end_index] + + @OddReal.setter + def OddReal(self, o): + start_index = self._current_offset + self._num_sub_ffts + end_index = self._current_offset + self._num_sub_ffts * 2 + buffer[start_index:end_index] = o + + @property + def OddImag(self) -> np.array: + start_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts + end_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts * 2 + return buffer[start_index:end_index] + + @OddImag.setter + def OddImag(self, o): + start_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts + end_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts * 2 + buffer[start_index: end_index] = o + + @property + def EvenReal(self) -> np.array: + start_index = self._current_offset + end_index = self._current_offset + self._num_sub_ffts + return buffer[start_index:end_index] + + @EvenReal.setter + def EvenReal(self, e): + start_index = self._current_offset + end_index = self._current_offset + self._num_sub_ffts + buffer[start_index:end_index] = e + + @property + def EvenImag(self) -> np.array: + start_index = self._imag_offset + self._current_offset + end_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts + return buffer[start_index:end_index] + + @EvenImag.setter + def EvenImag(self, e): + start_index = self._imag_offset + self._current_offset + end_index = self._imag_offset \ + + self._current_offset + self._num_sub_ffts + buffer[start_index:end_index] = e + + +# def mul_n(l_start: int, r_start: int, out_start, num: int): +# for i in range(num): +# buffer[out_start + i] = buffer[l_start + i] * buffer[r_start * i] + +def mul_32_factor(l_start: int, factor: float, out_start, num: int): + if num > 32: + raise ValueError('A maximum of 32 elements can be multiplied') + + for i in range(num): + buffer[out_start + i] = buffer[l_start + i] * factor + + +def add_32(l_start: int, r_start: int, out_start, num: int): + if num > 32: + raise ValueError('A maximum of 32 elements can be added') + + for i in range(num): + buffer[out_start + i] = buffer[l_start + i] + buffer[r_start * i] + + +def sub_32(l_start: int, r_start: int, out_start: int, num: int): + if num > 32: + raise ValueError('A maximum of 32 elements can be subtracted') + + for i in range(num): + buffer[out_start + i] = buffer[l_start + i] - buffer[r_start * i] + + +def bit_reverse(val: int, num_bits: int): + """Reverse the bits of an input integer. + + :param val: Value to reverse + :param num_bits: Number of bits to consider + """ + b = '{:0{width}b}'.format(val, width=num_bits) + return int(b[::-1], 2) + + +def fft_reorder(x: np.array, N: int) -> np.array: + """Reorder the elements of an array so that the indices are bit-inverse + values of the original ones. + """ + result = np.zeros(x.size) + for i in range(x.size): + result[i] = x[bit_reverse(i, int(np.log2(N)))] + return result + + +def fft(x_re: np.array, x_im: np.array) -> typing.Tuple[np.array, np.array]: + """Perform the fft algorithm on a given vector. + """ + imag_mem_offset = x_re.size + temp_mem_offset = x_re.size * 2 + + # Copy data to hardware buffer + for i in range(x_re.size): + buffer[i] = x_re[i] + # for i in range(x_im.size): + # buffer[i + imag_offset] = x_im[i] + + N0 = x_re.size + mem_wrapper = MemoryWrapper(N0, x_re.size) + + num_sub_ffts = N0 + while num_sub_ffts > 0: + N = N0 / num_sub_ffts + for k in range(int(N // 2)): + # mem_wrapper.set_kN(k, num_sub_ffts) + # + # factor_re = np.exp(-2j * np.pi * k / N).real + # factor_im = np.exp(-2j * np.pi * k / N).imag + # + # temp_re = factor_re * mem_wrapper.OddReal \ + # - factor_im * mem_wrapper.OddImag + # temp_im = factor_re * mem_wrapper.OddImag \ + # + factor_im * mem_wrapper.OddReal + # + # mem_wrapper.OddReal = mem_wrapper.EvenReal - temp_re + # mem_wrapper.OddImag = mem_wrapper.EvenImag - temp_im + # mem_wrapper.EvenReal = mem_wrapper.EvenReal + temp_re + # mem_wrapper.EvenImag = mem_wrapper.EvenImag + temp_im + + i = bit_reverse(k, int(np.log2(N))) + fft_offset = i * num_sub_ffts + odd_offset = num_sub_ffts + + factor_re = np.exp(-2j * np.pi * k / N).real + factor_im = np.exp(-2j * np.pi * k / N).imag + + # Compute temp variables + # TODO: Don't hardcode 32 + temp_re_start = temp_mem_offset + temp_im_start = temp_re_start + imag_mem_offset + even_re_start = fft_offset + even_im_start = fft_offset + imag_mem_offset + odd_re_start = fft_offset + odd_offset + odd_im_start = fft_offset + imag_mem_offset + + mul_32_factor(l_start=odd_re_start, factor=factor_re, + out_start=temp_re_start, num=4) + mul_32_factor(l_start=odd_im_start, factor=-factor_im, + out_start=temp_im_start, num=4) + add_32(l_start=temp_re_start, r_start=temp_mem_offset + 32, + out_start=temp_mem_offset, num=4) + + mul_32_factor(l_start=(fft_offset + odd_offset + imag_mem_offset), + factor=factor_re, + out_start=temp_mem_offset + 32, num=4) + mul_32_factor(l_start=(fft_offset + odd_offset + imag_mem_offset), + factor=factor_im, + out_start=temp_mem_offset + 64, num=4) + add_32(l_start=temp_mem_offset + 32, r_start=temp_mem_offset + 64, + out_start=temp_mem_offset + 32, num=4) + + num_sub_ffts = num_sub_ffts // 2 + + buffer[0:imag_mem_offset] = fft_reorder(buffer[0:imag_mem_offset], + imag_mem_offset) + buffer[imag_mem_offset:2 * imag_mem_offset] = fft_reorder( + buffer[imag_mem_offset:2 * imag_mem_offset], imag_mem_offset) + + return buffer[0:imag_mem_offset], buffer[ + imag_mem_offset:imag_mem_offset * 2] + + +def main(): + x_re = np.arange(8) + x_im = np.zeros(x_re.size) + X_re, X_im = fft(x_re, x_im) + + # for k, (X_k, X_exp_k) in enumerate(zip(X_re, np.fft.fft(x_re))): + # print( + # f"X[{k:02}]: {X_k.real:07.2f} + j({X_k.imag:07.2f})\t\tX_exp[" + # f"{k:02}]: {X_exp_k.real:07.2f} + j({X_exp_k.imag:07.2f})") + + for k, ((X_re_k, X_im_k), X_exp_k) in enumerate( + zip(zip(X_re, X_im), np.fft.fft(x_re))): + print( + f"X[{k:02}]: {X_re_k:07.2f} + j({X_im_k:07.2f})\t\tX_exp[" + f"{k:02}]: {X_exp_k.real:07.2f} + j({X_exp_k.imag:07.2f})") + + +if __name__ == "__main__": + main()