Copyright: 2010-2023, IOhannes m zmölnig <umlaeute@debian.org>
License: Expat
+Files: debian/missing-sources/Simple-FFT/*
+Copyright: 2013-2020, Dmitry Ivanov
+License: Expat
+
License: Expat
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the
--- /dev/null
+/**
+ * Copyright (c) 2013-2020 Dmitry Ivanov
+ *
+ * This file is a part of Simple-FFT project and is distributed under the terms
+ * of MIT license: https://opensource.org/licenses/MIT
+ */
+
+#ifndef __SIMPLE_FFT__CHECK_FFT_HPP__
+#define __SIMPLE_FFT__CHECK_FFT_HPP__
+
+#include "fft_settings.h"
+#include "error_handling.hpp"
+#include "copy_array.hpp"
+#include <cstddef>
+#include <cmath>
+#include <numeric>
+
+using std::size_t;
+
+namespace simple_fft {
+namespace check_fft_private {
+
+enum CheckMode
+{
+ CHECK_FFT_PARSEVAL,
+ CHECK_FFT_ENERGY,
+ CHECK_FFT_EQUALITY
+};
+
+template <class TArray1D, class TComplexArray1D>
+void getMaxAbsoluteAndRelativeErrorNorms(const TArray1D & array1,
+ const TComplexArray1D & array2, const size_t size,
+ real_type & max_absolute_error_norm,
+ real_type & max_relative_error_norm)
+{
+ using std::abs;
+
+ real_type current_error;
+
+ // NOTE: no parallelization here, it is a completely sequential loop!
+ for(size_t i = 0; i < size; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ current_error = abs(array1[i] - array2[i]);
+#else
+ current_error = abs(array1(i) - array2(i));
+#endif
+ if (current_error > max_absolute_error_norm) {
+ max_absolute_error_norm = current_error;
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ if (abs(array1[i]) > abs(array2[i])) {
+ max_relative_error_norm = (abs(array1[i]) > 1e-20
+ ? max_absolute_error_norm / abs(array1[i])
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2[i]) > 1e-20
+ ? max_absolute_error_norm / abs(array2[i])
+ : 0.0);
+ }
+#else
+ if (abs(array1(i)) > abs(array2(i))) {
+ max_relative_error_norm = (abs(array1(i)) > 1e-20
+ ? max_absolute_error_norm / abs(array1(i))
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2(i)) > 1e-20
+ ? max_absolute_error_norm / abs(array2(i))
+ : 0.0);
+ }
+#endif
+ }
+ }
+}
+
+template <class TArray2D, class TComplexArray2D>
+void getMaxAbsoluteAndRelativeErrorNorms(const TArray2D & array1,
+ const TComplexArray2D & array2,
+ const size_t size1, const size_t size2,
+ real_type & max_absolute_error_norm,
+ real_type & max_relative_error_norm)
+{
+ using std::abs;
+
+ real_type current_error;
+
+ // NOTE: no parallelization here, it is a completely sequential loop!
+ for(int i = 0; i < static_cast<int>(size1); ++i) {
+ for(int j = 0; j < static_cast<int>(size2); ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ current_error = abs(array1[i][j] - array2[i][j]);
+#else
+ current_error = abs(array1(i,j) - array2(i,j));
+#endif
+ if (current_error > max_absolute_error_norm) {
+ max_absolute_error_norm = current_error;
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ if (abs(array1[i][j]) > abs(array2[i][j])) {
+ max_relative_error_norm = (abs(array1[i][j]) > 1e-20
+ ? max_absolute_error_norm / abs(array1[i][j])
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2[i][j]) > 1e-20
+ ? max_absolute_error_norm / abs(array2[i][j])
+ : 0.0);
+ }
+#else
+ if (abs(array1(i,j)) > abs(array2(i,j))) {
+ max_relative_error_norm = (abs(array1(i,j)) > 1e-20
+ ? max_absolute_error_norm / abs(array1(i,j))
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2(i,j)) > 1e-20
+ ? max_absolute_error_norm / abs(array2(i,j))
+ : 0.0);
+ }
+#endif
+ }
+ }
+ }
+}
+
+template <class TArray3D, class TComplexArray3D>
+void getMaxAbsoluteAndRelativeErrorNorms(const TArray3D & array1, const TComplexArray3D & array2,
+ const size_t size1, const size_t size2,
+ const size_t size3, real_type & max_absolute_error_norm,
+ real_type & max_relative_error_norm)
+{
+ using std::abs;
+
+ real_type current_error;
+
+ // NOTE: no parallelization here, it is a completely sequential loop!
+ for(int i = 0; i < static_cast<int>(size1); ++i) {
+ for(int j = 0; j < static_cast<int>(size2); ++j) {
+ for(int k = 0; k < static_cast<int>(size3); ++k) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ current_error = abs(array1[i][j][k] - array2[i][j][k]);
+#else
+ current_error = abs(array1(i,j,k) - array2(i,j,k));
+#endif
+ if (current_error > max_absolute_error_norm) {
+ max_absolute_error_norm = current_error;
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ if (abs(array1[i][j][k]) > abs(array2[i][j][k])) {
+ max_relative_error_norm = (abs(array1[i][j][k]) > 1e-20
+ ? max_absolute_error_norm / abs(array1[i][j][k])
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2[i][j][k]) > 1e-20
+ ? max_absolute_error_norm / abs(array2[i][j][k])
+ : 0.0);
+ }
+#else
+ if (abs(array1(i,j,k)) > abs(array2(i,j,k))) {
+ max_relative_error_norm = (abs(array1(i,j,k)) > 1e-20
+ ? max_absolute_error_norm / abs(array1(i,j,k))
+ : 0.0);
+ }
+ else {
+ max_relative_error_norm = (abs(array2(i,j,k)) > 1e-20
+ ? max_absolute_error_norm / abs(array2(i,j,k))
+ : 0.0);
+ }
+#endif
+ }
+ }
+ }
+ }
+}
+
+template <class TArray1D>
+real_type squareAbsAccumulate(const TArray1D & array, const size_t size,
+ const real_type init)
+{
+ int size_signed = static_cast<int>(size);
+ real_type sum = init;
+
+ using std::abs;
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for reduction(+:sum)
+#endif
+#endif
+ for(int i = 0; i < size_signed; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ sum += abs(array[i] * array[i]);
+#else
+ sum += abs(array(i) * array(i));
+#endif
+ }
+
+ return sum;
+}
+
+template <class TArray2D>
+real_type squareAbsAccumulate(const TArray2D & array, const size_t size1,
+ const size_t size2, const real_type init)
+{
+ int size1_signed = static_cast<int>(size1);
+ int size2_signed = static_cast<int>(size2);
+ real_type sum = init;
+
+ using std::abs;
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for reduction(+:sum)
+#endif
+#endif
+ for(int i = 0; i < size1_signed; ++i) {
+ for(int j = 0; j < size2_signed; ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ sum += abs(array[i][j] * array[i][j]);
+#else
+ sum += abs(array(i,j) * array(i,j));
+#endif
+ }
+ }
+
+ return sum;
+}
+
+template <class TArray3D>
+real_type squareAbsAccumulate(const TArray3D & array, const size_t size1,
+ const size_t size2, const size_t size3,
+ const real_type init)
+{
+ int size1_signed = static_cast<int>(size1);
+ int size2_signed = static_cast<int>(size2);
+ int size3_signed = static_cast<int>(size3);
+ real_type sum = init;
+
+ using std::abs;
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for reduction(+:sum)
+#endif
+#endif
+ for(int i = 0; i < size1_signed; ++i) {
+ for(int j = 0; j < size2_signed; ++j) {
+ for(int k = 0; k < size3_signed; ++k) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ sum += abs(array[i][j][k] * array[i][j][k]);
+#else
+ sum += abs(array(i,j,k) * array(i,j,k));
+#endif
+ }
+ }
+ }
+
+ return sum;
+}
+
+// Generic template for CCheckFFT struct followed by its explicit specializations
+// for certain numbers of dimensions. TArray can be either of real or complex type.
+// The technique is similar to the one applied for CFFT struct.
+template <class TArray, class TComplexArray, int NumDims>
+struct CCheckFFT
+{};
+
+template <class TArray1D, class TComplexArray1D>
+struct CCheckFFT<TArray1D,TComplexArray1D,1>
+{
+ static bool check_fft(const TArray1D & data_before,
+ const TComplexArray1D & data_after,
+ const size_t size, const real_type relative_tolerance,
+ real_type & discrepancy, const CheckMode check_mode,
+ const char *& error_description)
+ {
+ using namespace error_handling;
+
+ if(0 == size) {
+ GetErrorDescription(EC_NUM_OF_ELEMS_IS_ZERO, error_description);
+ return false;
+ }
+
+ if ( (CHECK_FFT_PARSEVAL != check_mode) &&
+ (CHECK_FFT_ENERGY != check_mode) &&
+ (CHECK_FFT_EQUALITY != check_mode) )
+ {
+ GetErrorDescription(EC_WRONG_CHECK_FFT_MODE, error_description);
+ return false;
+ }
+
+ if (CHECK_FFT_EQUALITY != check_mode)
+ {
+ real_type sum_before = squareAbsAccumulate<TArray1D>(data_before, size, 0.0);
+ real_type sum_after = squareAbsAccumulate<TComplexArray1D>(data_after, size, 0.0);
+
+ if (CHECK_FFT_PARSEVAL == check_mode) {
+ sum_after /= size;
+ }
+
+ using std::abs;
+
+ discrepancy = abs(sum_before - sum_after);
+
+ if (discrepancy / ((sum_before < 1e-20) ? (sum_before + 1e-20) : sum_before) > relative_tolerance) {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ else {
+ return true;
+ }
+ }
+ else {
+ real_type relative_error;
+ getMaxAbsoluteAndRelativeErrorNorms(data_before, data_after, size,
+ discrepancy, relative_error);
+ if (relative_error < relative_tolerance) {
+ return true;
+ }
+ else {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ }
+ }
+};
+
+template <class TArray2D, class TComplexArray2D>
+struct CCheckFFT<TArray2D,TComplexArray2D,2>
+{
+ static bool check_fft(const TArray2D & data_before,
+ const TComplexArray2D & data_after,
+ const size_t size1, const size_t size2,
+ const real_type relative_tolerance, real_type & discrepancy,
+ const CheckMode check_mode, const char *& error_description)
+ {
+ using namespace error_handling;
+
+ if( (0 == size1) || (0 == size2) ) {
+ GetErrorDescription(EC_NUM_OF_ELEMS_IS_ZERO, error_description);
+ return false;
+ }
+
+ if ( (CHECK_FFT_PARSEVAL != check_mode) &&
+ (CHECK_FFT_ENERGY != check_mode) &&
+ (CHECK_FFT_EQUALITY != check_mode) )
+ {
+ GetErrorDescription(EC_WRONG_CHECK_FFT_MODE, error_description);
+ return false;
+ }
+
+ if (CHECK_FFT_EQUALITY != check_mode)
+ {
+ real_type sum_before = squareAbsAccumulate<TArray2D>(data_before, size1, size2, 0.0);
+ real_type sum_after = squareAbsAccumulate<TComplexArray2D>(data_after, size1, size2, 0.0);
+
+ if (CHECK_FFT_PARSEVAL == check_mode) {
+ sum_after /= size1 * size2;
+ }
+
+ using std::abs;
+
+ discrepancy = abs(sum_before - sum_after);
+
+ if (discrepancy / ((sum_before < 1e-20) ? (sum_before + 1e-20) : sum_before) > relative_tolerance) {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ else {
+ return true;
+ }
+ }
+ else {
+ real_type relative_error;
+ getMaxAbsoluteAndRelativeErrorNorms(data_before, data_after, size1,
+ size2, discrepancy, relative_error);
+ if (relative_error < relative_tolerance) {
+ return true;
+ }
+ else {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ }
+ }
+};
+
+template <class TArray3D, class TComplexArray3D>
+struct CCheckFFT<TArray3D,TComplexArray3D,3>
+{
+ static bool check_fft(const TArray3D & data_before,
+ const TComplexArray3D & data_after,
+ const size_t size1, const size_t size2, const size_t size3,
+ const real_type relative_tolerance, real_type & discrepancy,
+ const CheckMode check_mode, const char *& error_description)
+ {
+ using namespace error_handling;
+
+ if( (0 == size1) || (0 == size2) || (0 == size3) ) {
+ GetErrorDescription(EC_NUM_OF_ELEMS_IS_ZERO, error_description);
+ return false;
+ }
+
+ if ( (CHECK_FFT_PARSEVAL != check_mode) &&
+ (CHECK_FFT_ENERGY != check_mode) &&
+ (CHECK_FFT_EQUALITY != check_mode) )
+ {
+ GetErrorDescription(EC_WRONG_CHECK_FFT_MODE, error_description);
+ return false;
+ }
+
+ if (CHECK_FFT_EQUALITY != check_mode)
+ {
+ real_type sum_before = squareAbsAccumulate<TArray3D>(data_before, size1, size2, size3, 0.0);
+ real_type sum_after = squareAbsAccumulate<TComplexArray3D>(data_after, size1, size2, size3, 0.0);
+
+ if (CHECK_FFT_PARSEVAL == check_mode) {
+ sum_after /= size1 * size2 * size3;
+ }
+
+ using std::abs;
+
+ discrepancy = abs(sum_before - sum_after);
+
+ if (discrepancy / ((sum_before < 1e-20) ? (sum_before + 1e-20) : sum_before) > relative_tolerance) {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ else {
+ return true;
+ }
+ }
+ else {
+ real_type relative_error;
+ getMaxAbsoluteAndRelativeErrorNorms(data_before, data_after, size1,
+ size2, size3, discrepancy, relative_error);
+ if (relative_error < relative_tolerance) {
+ return true;
+ }
+ else {
+ GetErrorDescription(EC_RELATIVE_ERROR_TOO_LARGE, error_description);
+ return false;
+ }
+ }
+ }
+};
+
+} // namespace check_fft_private
+
+namespace check_fft {
+
+template <class TArray1D, class TComplexArray1D>
+bool checkParsevalTheorem(const TArray1D & data_before_FFT,
+ const TComplexArray1D & data_after_FFT,
+ const size_t size, const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray1D,TComplexArray1D,1>::check_fft(data_before_FFT,
+ data_after_FFT, size, relative_tolerance,
+ discrepancy, check_fft_private::CHECK_FFT_PARSEVAL,
+ error_description);
+}
+
+template <class TArray2D, class TComplexArray2D>
+bool checkParsevalTheorem(const TArray2D & data_before_FFT,
+ const TComplexArray2D & data_after_FFT,
+ const size_t size1, const size_t size2,
+ const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray2D,TComplexArray2D,2>::check_fft(data_before_FFT,
+ data_after_FFT, size1, size2, relative_tolerance,
+ discrepancy, check_fft_private::CHECK_FFT_PARSEVAL,
+ error_description);
+}
+
+template <class TArray3D, class TComplexArray3D>
+bool checkParsevalTheorem(const TArray3D & data_before_FFT,
+ const TComplexArray3D & data_after_FFT,
+ const size_t size1, const size_t size2, const size_t size3,
+ const real_type relative_tolerance, real_type & discrepancy,
+ const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray3D,TComplexArray3D,3>::check_fft(data_before_FFT,
+ data_after_FFT, size1, size2, size3,
+ relative_tolerance, discrepancy,
+ check_fft_private::CHECK_FFT_PARSEVAL,
+ error_description);
+}
+
+template <class TArray1D, class TComplexArray1D>
+bool checkEnergyConservation(const TArray1D & data_before_FFT,
+ const TComplexArray1D & data_after_FFT_and_IFFT,
+ const size_t size, const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray1D,TComplexArray1D,1>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size, relative_tolerance,
+ discrepancy, check_fft_private::CHECK_FFT_ENERGY,
+ error_description);
+}
+
+template <class TArray2D, class TComplexArray2D>
+bool checkEnergyConservation(const TArray2D & data_before_FFT,
+ const TComplexArray2D & data_after_FFT_and_IFFT,
+ const size_t size1, const size_t size2,
+ const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray2D,TComplexArray2D,2>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size1, size2,
+ relative_tolerance, discrepancy,
+ check_fft_private::CHECK_FFT_ENERGY,
+ error_description);
+}
+
+template <class TArray3D, class TComplexArray3D>
+bool checkEnergyConservation(const TArray3D & data_before_FFT,
+ const TComplexArray3D & data_after_FFT_and_IFFT,
+ const size_t size1, const size_t size2, const size_t size3,
+ const real_type relative_tolerance, real_type & discrepancy,
+ const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray3D,TComplexArray3D,3>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size1, size2,
+ size3, relative_tolerance, discrepancy,
+ check_fft_private::CHECK_FFT_ENERGY,
+ error_description);
+}
+
+template <class TArray1D, class TComplexArray1D>
+bool checkEquality(const TArray1D & data_before_FFT,
+ const TComplexArray1D & data_after_FFT_and_IFFT,
+ const size_t size, const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray1D,TComplexArray1D,1>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size, relative_tolerance,
+ discrepancy, check_fft_private::CHECK_FFT_EQUALITY,
+ error_description);
+}
+
+template <class TArray2D, class TComplexArray2D>
+bool checkEquality(const TArray2D & data_before_FFT,
+ const TComplexArray2D & data_after_FFT_and_IFFT, const size_t size1,
+ const size_t size2, const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray2D,TComplexArray2D,2>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size1, size2,
+ relative_tolerance, discrepancy,
+ check_fft_private::CHECK_FFT_EQUALITY,
+ error_description);
+}
+
+template <class TArray3D, class TComplexArray3D>
+bool checkEquality(const TArray3D & data_before_FFT,
+ const TComplexArray3D & data_after_FFT_and_IFFT, const size_t size1,
+ const size_t size2, const size_t size3, const real_type relative_tolerance,
+ real_type & discrepancy, const char *& error_description)
+{
+ return check_fft_private::CCheckFFT<TArray3D,TComplexArray3D,3>::check_fft(data_before_FFT,
+ data_after_FFT_and_IFFT, size1, size2,
+ size3, relative_tolerance, discrepancy,
+ check_fft_private::CHECK_FFT_EQUALITY,
+ error_description);
+}
+
+} // namespace check_fft
+} // namespace simple_fft
+
+#endif // __SIMPLE_FFT__CHECK_FFT_HPP__
--- /dev/null
+/**\r
+ * Copyright (c) 2013-2020 Dmitry Ivanov\r
+ *\r
+ * This file is a part of Simple-FFT project and is distributed under the terms\r
+ * of MIT license: https://opensource.org/licenses/MIT\r
+ */\r
+\r
+#ifndef __SIMPLE_FFT__COPY_ARRAY_HPP\r
+#define __SIMPLE_FFT__COPY_ARRAY_HPP\r
+\r
+#include "fft_settings.h"\r
+#include "error_handling.hpp"\r
+#include <cstddef>\r
+\r
+using std::size_t;\r
+\r
+namespace simple_fft {\r
+namespace copy_array {\r
+\r
+template <class TComplexArray1D>\r
+void copyArray(const TComplexArray1D & data_from, TComplexArray1D & data_to,\r
+ const size_t size)\r
+{\r
+ int size_signed = static_cast<int>(size);\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size_signed; ++i) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i] = data_from[i];\r
+#else\r
+ data_to(i) = data_from(i);\r
+#endif\r
+ }\r
+}\r
+\r
+template <class TComplexArray1D, class TRealArray1D>\r
+void copyArray(const TRealArray1D & data_from, TComplexArray1D & data_to,\r
+ const size_t size)\r
+{\r
+ int size_signed = static_cast<int>(size);\r
+\r
+ // NOTE: user's complex type should have constructor like\r
+ // "complex(real, imag)", where each of real and imag has\r
+ // real type.\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size_signed; ++i) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i] = complex_type(data_from[i], 0.0);\r
+#else\r
+ data_to(i) = complex_type(data_from(i), 0.0);\r
+#endif\r
+ }\r
+}\r
+\r
+template <class TComplexArray2D>\r
+void copyArray(const TComplexArray2D & data_from, TComplexArray2D & data_to,\r
+ const size_t size1, const size_t size2)\r
+{\r
+ int size1_signed = static_cast<int>(size1);\r
+ int size2_signed = static_cast<int>(size2);\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size1_signed; ++i) {\r
+ for(int j = 0; j < size2_signed; ++j) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i][j] = data_from[i][j];\r
+#else\r
+ data_to(i,j) = data_from(i,j);\r
+#endif\r
+ }\r
+ }\r
+}\r
+\r
+template <class TComplexArray2D, class TRealArray2D>\r
+void copyArray(const TRealArray2D & data_from, TComplexArray2D & data_to,\r
+ const size_t size1, const size_t size2)\r
+{\r
+ int size1_signed = static_cast<int>(size1);\r
+ int size2_signed = static_cast<int>(size2);\r
+\r
+ // NOTE: user's complex type should have constructor like\r
+ // "complex(real, imag)", where each of real and imag has\r
+ // real type.\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size1_signed; ++i) {\r
+ for(int j = 0; j < size2_signed; ++j) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i][j] = complex_type(data_from[i][j], 0.0);\r
+#else\r
+ data_to(i,j) = complex_type(data_from(i,j), 0.0);\r
+#endif\r
+ }\r
+ }\r
+}\r
+\r
+template <class TComplexArray3D>\r
+void copyArray(const TComplexArray3D & data_from, TComplexArray3D & data_to,\r
+ const size_t size1, const size_t size2, const size_t size3)\r
+{\r
+ int size1_signed = static_cast<int>(size1);\r
+ int size2_signed = static_cast<int>(size2);\r
+ int size3_signed = static_cast<int>(size3);\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size1_signed; ++i) {\r
+ for(int j = 0; j < size2_signed; ++j) {\r
+ for(int k = 0; k < size3_signed; ++k) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i][j][k] = data_from[i][j][k];\r
+#else\r
+ data_to(i,j,k) = data_from(i,j,k);\r
+#endif\r
+ }\r
+ }\r
+ }\r
+}\r
+\r
+template <class TComplexArray3D, class TRealArray3D>\r
+void copyArray(const TRealArray3D & data_from, TComplexArray3D & data_to,\r
+ const size_t size1, const size_t size2, const size_t size3)\r
+{\r
+ int size1_signed = static_cast<int>(size1);\r
+ int size2_signed = static_cast<int>(size2);\r
+ int size3_signed = static_cast<int>(size3);\r
+\r
+ // NOTE: user's complex type should have constructor like\r
+ // "complex(real, imag)", where each of real and imag has\r
+ // real type.\r
+\r
+#ifndef __clang__\r
+#ifdef __USE_OPENMP\r
+#pragma omp parallel for\r
+#endif\r
+#endif\r
+ for(int i = 0; i < size1_signed; ++i) {\r
+ for(int j = 0; j < size2_signed; ++j) {\r
+ for(int k = 0; k < size3_signed; ++k) {\r
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+ data_to[i][j][k] = complex_type(data_from[i][j][k], 0.0);\r
+#else\r
+ data_to(i,j,k) = complex_type(data_from(i,j,k), 0.0);\r
+#endif\r
+ }\r
+ }\r
+ }\r
+}\r
+\r
+} // namespace copy_array\r
+} // namespace simple_fft\r
+\r
+#endif // __SIMPLE_FFT__COPY_ARRAY_HPP\r
--- /dev/null
+/**\r
+ * Copyright (c) 2013-2020 Dmitry Ivanov\r
+ *\r
+ * This file is a part of Simple-FFT project and is distributed under the terms\r
+ * of MIT license: https://opensource.org/licenses/MIT\r
+ */\r
+\r
+#ifndef __SIMPLE_FFT__ERROR_HANDLING_HPP\r
+#define __SIMPLE_FFT__ERROR_HANDLING_HPP\r
+\r
+namespace simple_fft {\r
+namespace error_handling {\r
+\r
+enum EC_SimpleFFT\r
+{\r
+ EC_SUCCESS = 0,\r
+ EC_UNSUPPORTED_DIMENSIONALITY,\r
+ EC_WRONG_FFT_DIRECTION,\r
+ EC_ONE_OF_DIMS_ISNT_POWER_OF_TWO,\r
+ EC_NUM_OF_ELEMS_IS_ZERO,\r
+ EC_WRONG_CHECK_FFT_MODE,\r
+ EC_RELATIVE_ERROR_TOO_LARGE\r
+};\r
+\r
+inline void GetErrorDescription(const EC_SimpleFFT error_code,\r
+ const char *& error_description)\r
+{\r
+ switch(error_code)\r
+ {\r
+ case EC_SUCCESS:\r
+ error_description = "Calculation was successful!";\r
+ break;\r
+ case EC_UNSUPPORTED_DIMENSIONALITY:\r
+ error_description = "Unsupported dimensionality: currently only 1D, 2D "\r
+ "and 3D arrays are supported";\r
+ break;\r
+ case EC_WRONG_FFT_DIRECTION:\r
+ error_description = "Wrong direction for FFT was specified";\r
+ break;\r
+ case EC_ONE_OF_DIMS_ISNT_POWER_OF_TWO:\r
+ error_description = "Unsupported dimensionality: one of dimensions is not "\r
+ "a power of 2";\r
+ break;\r
+ case EC_NUM_OF_ELEMS_IS_ZERO:\r
+ error_description = "Number of elements for FFT or IFFT is zero!";\r
+ break;\r
+ case EC_WRONG_CHECK_FFT_MODE:\r
+ error_description = "Wrong check FFT mode was specified (should be either "\r
+ "Parseval theorem or energy conservation check";\r
+ break;\r
+ case EC_RELATIVE_ERROR_TOO_LARGE:\r
+ error_description = "Relative error returned by FFT test exceeds specified "\r
+ "relative tolerance";\r
+ break;\r
+ default:\r
+ error_description = "Unknown error";\r
+ break;\r
+ }\r
+}\r
+\r
+} // namespace error_handling\r
+} // namespace simple_fft\r
+\r
+#endif // __SIMPLE_FFT__ERROR_HANDLING_HPP\r
--- /dev/null
+/**
+ * Copyright (c) 2013-2020 Dmitry Ivanov
+ *
+ * This file is a part of Simple-FFT project and is distributed under the terms
+ * of MIT license: https://opensource.org/licenses/MIT
+ */
+
+#ifndef __SIMPLE_FFT__FFT_H__
+#define __SIMPLE_FFT__FFT_H__
+
+#include <cstddef>
+
+using std::size_t;
+
+/// The public API
+namespace simple_fft {
+
+/// FFT and IFFT functions
+
+// in-place, complex, forward
+template <class TComplexArray1D>
+bool FFT(TComplexArray1D & data, const size_t size, const char *& error_description);
+
+template <class TComplexArray2D>
+bool FFT(TComplexArray2D & data, const size_t size1, const size_t size2,
+ const char *& error_description);
+
+template <class TComplexArray3D>
+bool FFT(TComplexArray3D & data, const size_t size1, const size_t size2, const size_t size3,
+ const char *& error_description);
+
+// in-place, complex, inverse
+template <class TComplexArray1D>
+bool IFFT(TComplexArray1D & data, const size_t size, const char *& error_description);
+
+template <class TComplexArray2D>
+bool IFFT(TComplexArray2D & data, const size_t size1, const size_t size2,
+ const char *& error_description);
+
+template <class TComplexArray3D>
+bool IFFT(TComplexArray3D & data, const size_t size1, const size_t size2, const size_t size3,
+ const char *& error_description);
+
+// not-in-place, complex, forward
+template <class TComplexArray1D>
+bool FFT(const TComplexArray1D & data_in, TComplexArray1D & data_out,
+ const size_t size, const char *& error_description);
+
+template <class TComplexArray2D>
+bool FFT(const TComplexArray2D & data_in, TComplexArray2D & data_out,
+ const size_t size1, const size_t size2, const char *& error_description);
+
+template <class TComplexArray3D>
+bool FFT(const TComplexArray3D & data_in, TComplexArray3D & data_out,
+ const size_t size1, const size_t size2, const size_t size3,
+ const char *& error_description);
+
+// not-in-place, complex, inverse
+template <class TComplexArray1D>
+bool IFFT(const TComplexArray1D & data_in, TComplexArray1D & data_out,
+ const size_t size, const char *& error_description);
+
+template <class TComplexArray2D>
+bool IFFT(const TComplexArray2D & data_in, TComplexArray2D & data_out,
+ const size_t size1, const size_t size2, const char *& error_description);
+
+template <class TComplexArray3D>
+bool IFFT(const TComplexArray3D & data_in, TComplexArray3D & data_out,
+ const size_t size1, const size_t size2, const size_t size3,
+ const char *& error_description);
+
+// not-in-place, real, forward
+template <class TRealArray1D, class TComplexArray1D>
+bool FFT(const TRealArray1D & data_in, TComplexArray1D & data_out,
+ const size_t size, const char *& error_description);
+
+template <class TRealArray2D, class TComplexArray2D>
+bool FFT(const TRealArray2D & data_in, TComplexArray2D & data_out,
+ const size_t size1, const size_t size2, const char *& error_description);
+
+template <class TRealArray3D, class TComplexArray3D>
+bool FFT(const TRealArray3D & data_in, TComplexArray3D & data_out,
+ const size_t size1, const size_t size2, const size_t size3,
+ const char *& error_description);
+
+// NOTE: There is no inverse transform from complex spectrum to real signal
+// because round-off errors during computation of inverse FFT lead to the appearance
+// of signal imaginary components even though they are small by absolute value.
+// These can be ignored but the author of this file thinks adding such an function
+// would be wrong methodogically: looking at complex result, you can estimate
+// the value of spurious imaginary part. Otherwise you may never know that IFFT
+// provides too large imaginary values due to too small grid size, for example.
+
+} // namespace simple_fft
+
+#endif // __SIMPLE_FFT__FFT_H__
+
+#include "fft.hpp"
--- /dev/null
+/**\r
+ * Copyright (c) 2013-2020 Dmitry Ivanov\r
+ *\r
+ * This file is a part of Simple-FFT project and is distributed under the terms\r
+ * of MIT license: https://opensource.org/licenses/MIT\r
+ */\r
+\r
+#ifndef __SIMPLE_FFT__FFT_HPP__\r
+#define __SIMPLE_FFT__FFT_HPP__\r
+\r
+#include "copy_array.hpp"\r
+#include "fft_impl.hpp"\r
+\r
+namespace simple_fft {\r
+\r
+// in-place, complex, forward\r
+template <class TComplexArray1D>\r
+bool FFT(TComplexArray1D & data, const size_t size, const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray1D,1>::FFT_inplace(data, size, impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray2D>\r
+bool FFT(TComplexArray2D & data, const size_t size1, const size_t size2,\r
+ const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray2D,2>::FFT_inplace(data, size1, size2, impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray3D>\r
+bool FFT(TComplexArray3D & data, const size_t size1, const size_t size2, const size_t size3,\r
+ const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray3D,3>::FFT_inplace(data, size1, size2, size3,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+// in-place, complex, inverse\r
+template <class TComplexArray1D>\r
+bool IFFT(TComplexArray1D & data, const size_t size, const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray1D,1>::FFT_inplace(data, size, impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray2D>\r
+bool IFFT(TComplexArray2D & data, const size_t size1, const size_t size2,\r
+ const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray2D,2>::FFT_inplace(data, size1, size2, impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray3D>\r
+bool IFFT(TComplexArray3D & data, const size_t size1, const size_t size2, const size_t size3,\r
+ const char *& error_description)\r
+{\r
+ return impl::CFFT<TComplexArray3D,3>::FFT_inplace(data, size1, size2, size3,\r
+ impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+// not-in-place, complex, forward\r
+template <class TComplexArray1D>\r
+bool FFT(const TComplexArray1D & data_in, TComplexArray1D & data_out,\r
+ const size_t size, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size);\r
+ return impl::CFFT<TComplexArray1D,1>::FFT_inplace(data_out, size, impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray2D>\r
+bool FFT(const TComplexArray2D & data_in, TComplexArray2D & data_out,\r
+ const size_t size1, const size_t size2, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2);\r
+ return impl::CFFT<TComplexArray2D,2>::FFT_inplace(data_out, size1, size2,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray3D>\r
+bool FFT(const TComplexArray3D & data_in, TComplexArray3D & data_out,\r
+ const size_t size1, const size_t size2, const size_t size3,\r
+ const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2, size3);\r
+ return impl::CFFT<TComplexArray3D,3>::FFT_inplace(data_out, size1, size2, size3,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+// not-in-place, complex, inverse\r
+template <class TComplexArray1D>\r
+bool IFFT(const TComplexArray1D & data_in, TComplexArray1D & data_out,\r
+ const size_t size, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size);\r
+ return impl::CFFT<TComplexArray1D,1>::FFT_inplace(data_out, size, impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray2D>\r
+bool IFFT(const TComplexArray2D & data_in, TComplexArray2D & data_out,\r
+ const size_t size1, const size_t size2, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2);\r
+ return impl::CFFT<TComplexArray2D,2>::FFT_inplace(data_out, size1, size2,\r
+ impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TComplexArray3D>\r
+bool IFFT(const TComplexArray3D & data_in, TComplexArray3D & data_out,\r
+ const size_t size1, const size_t size2, const size_t size3,\r
+ const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2, size3);\r
+ return impl::CFFT<TComplexArray3D,3>::FFT_inplace(data_out, size1, size2, size3,\r
+ impl::FFT_BACKWARD,\r
+ error_description);\r
+}\r
+\r
+// not-in-place, real, forward\r
+template <class TRealArray1D, class TComplexArray1D>\r
+bool FFT(const TRealArray1D & data_in, TComplexArray1D & data_out,\r
+ const size_t size, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size);\r
+ return impl::CFFT<TComplexArray1D,1>::FFT_inplace(data_out, size,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TRealArray2D, class TComplexArray2D>\r
+bool FFT(const TRealArray2D & data_in, TComplexArray2D & data_out,\r
+ const size_t size1, const size_t size2, const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2);\r
+ return impl::CFFT<TComplexArray2D,2>::FFT_inplace(data_out, size1, size2,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+template <class TRealArray3D, class TComplexArray3D>\r
+bool FFT(const TRealArray3D & data_in, TComplexArray3D & data_out,\r
+ const size_t size1, const size_t size2, const size_t size3,\r
+ const char *& error_description)\r
+{\r
+ copy_array::copyArray(data_in, data_out, size1, size2, size3);\r
+ return impl::CFFT<TComplexArray3D,3>::FFT_inplace(data_out, size1, size2, size3,\r
+ impl::FFT_FORWARD,\r
+ error_description);\r
+}\r
+\r
+} // simple_fft\r
+\r
+#endif // __SIMPLE_FFT__FFT_HPP__\r
--- /dev/null
+/**
+ * Copyright (c) 2013-2020 Dmitry Ivanov
+ *
+ * This file is a part of Simple-FFT project and is distributed under the terms
+ * of MIT license: https://opensource.org/licenses/MIT
+ */
+
+#ifndef __SIMPLE_FFT__FFT_IMPL_HPP__
+#define __SIMPLE_FFT__FFT_IMPL_HPP__
+
+#include "fft_settings.h"
+#include "error_handling.hpp"
+#include <cstddef>
+#include <math.h>
+#include <vector>
+
+using std::size_t;
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932
+#endif
+
+namespace simple_fft {
+namespace impl {
+
+enum FFT_direction
+{
+ FFT_FORWARD = 0,
+ FFT_BACKWARD
+};
+
+// checking whether the size of array dimension is power of 2
+// via "complement and compare" method
+inline bool isPowerOfTwo(const size_t num)
+{
+ return num && (!(num & (num - 1)));
+}
+
+inline bool checkNumElements(const size_t num_elements, const char *& error_description)
+{
+ using namespace error_handling;
+
+ if (!isPowerOfTwo(num_elements)) {
+ GetErrorDescription(EC_ONE_OF_DIMS_ISNT_POWER_OF_TWO, error_description);
+ return false;
+ }
+
+ return true;
+}
+
+template <class TComplexArray1D>
+inline void scaleValues(TComplexArray1D & data, const size_t num_elements)
+{
+ real_type mult = 1.0 / num_elements;
+ int num_elements_signed = static_cast<int>(num_elements);
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < num_elements_signed; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i] *= mult;
+#else
+ data(i) *= mult;
+#endif
+ }
+}
+
+// NOTE: explicit template specialization for the case of std::vector<complex_type>
+// because it is used in 2D and 3D FFT for both array classes with square and round
+// brackets of element access operator; I need to guarantee that sub-FFT 1D will
+// use square brackets for element access operator anyway. It is pretty ugly
+// to duplicate the code but I haven't found more elegant solution.
+template <>
+inline void scaleValues<std::vector<complex_type> >(std::vector<complex_type> & data,
+ const size_t num_elements)
+{
+ real_type mult = 1.0 / num_elements;
+ int num_elements_signed = static_cast<int>(num_elements);
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < num_elements_signed; ++i) {
+ data[i] *= mult;
+ }
+}
+
+template <class TComplexArray1D>
+inline void bufferExchangeHelper(TComplexArray1D & data, const size_t index_from,
+ const size_t index_to, complex_type & buf)
+{
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ buf = data[index_from];
+ data[index_from] = data[index_to];
+ data[index_to]= buf;
+#else
+ buf = data(index_from);
+ data(index_from) = data(index_to);
+ data(index_to)= buf;
+#endif
+}
+
+// NOTE: explicit template specialization for the case of std::vector<complex_type>
+// because it is used in 2D and 3D FFT for both array classes with square and round
+// brackets of element access operator; I need to guarantee that sub-FFT 1D will
+// use square brackets for element access operator anyway. It is pretty ugly
+// to duplicate the code but I haven't found more elegant solution.
+template <>
+inline void bufferExchangeHelper<std::vector<complex_type> >(std::vector<complex_type> & data,
+ const size_t index_from,
+ const size_t index_to,
+ complex_type & buf)
+{
+ buf = data[index_from];
+ data[index_from] = data[index_to];
+ data[index_to]= buf;
+}
+
+template <class TComplexArray1D>
+void rearrangeData(TComplexArray1D & data, const size_t num_elements)
+{
+ complex_type buf;
+
+ size_t target_index = 0;
+ size_t bit_mask;
+
+ for (size_t i = 0; i < num_elements; ++i)
+ {
+ if (target_index > i)
+ {
+ bufferExchangeHelper(data, target_index, i, buf);
+ }
+
+ // Initialize the bit mask
+ bit_mask = num_elements;
+
+ // While bit is 1
+ while (target_index & (bit_mask >>= 1)) // bit_mask = bit_mask >> 1
+ {
+ // Drop bit:
+ // & is bitwise AND,
+ // ~ is bitwise NOT
+ target_index &= ~bit_mask; // target_index = target_index & (~bit_mask)
+ }
+
+ // | is bitwise OR
+ target_index |= bit_mask; // target_index = target_index | bit_mask
+ }
+}
+
+template <class TComplexArray1D>
+inline void fftTransformHelper(TComplexArray1D & data, const size_t match,
+ const size_t k, complex_type & product,
+ const complex_type factor)
+{
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ product = data[match] * factor;
+ data[match] = data[k] - product;
+ data[k] += product;
+#else
+ product = data(match) * factor;
+ data(match) = data(k) - product;
+ data(k) += product;
+#endif
+}
+
+// NOTE: explicit template specialization for the case of std::vector<complex_type>
+// because it is used in 2D and 3D FFT for both array classes with square and round
+// brackets of element access operator; I need to guarantee that sub-FFT 1D will
+// use square brackets for element access operator anyway. It is pretty ugly
+// to duplicate the code but I haven't found more elegant solution.
+template <>
+inline void fftTransformHelper<std::vector<complex_type> >(std::vector<complex_type> & data,
+ const size_t match,
+ const size_t k,
+ complex_type & product,
+ const complex_type factor)
+{
+ product = data[match] * factor;
+ data[match] = data[k] - product;
+ data[k] += product;
+}
+
+template <class TComplexArray1D>
+bool makeTransform(TComplexArray1D & data, const size_t num_elements,
+ const FFT_direction fft_direction, const char *& error_description)
+{
+ using namespace error_handling;
+ using std::sin;
+
+ double local_pi;
+ switch(fft_direction)
+ {
+ case(FFT_FORWARD):
+ local_pi = -M_PI;
+ break;
+ case(FFT_BACKWARD):
+ local_pi = M_PI;
+ break;
+ default:
+ GetErrorDescription(EC_WRONG_FFT_DIRECTION, error_description);
+ return false;
+ }
+
+ // declare variables to cycle the bits of initial signal
+ size_t next, match;
+ real_type sine;
+ real_type delta;
+ complex_type mult, factor, product;
+
+ // NOTE: user's complex type should have constructor like
+ // "complex(real, imag)", where each of real and imag has
+ // real type.
+
+ // cycle for all bit positions of initial signal
+ for (size_t i = 1; i < num_elements; i <<= 1)
+ {
+ next = i << 1; // getting the next bit
+ delta = local_pi / i; // angle increasing
+ sine = sin(0.5 * delta); // supplementary sin
+ // multiplier for trigonometric recurrence
+ mult = complex_type(-2.0 * sine * sine, sin(delta));
+ factor = 1.0; // start transform factor
+
+ for (size_t j = 0; j < i; ++j) // iterations through groups
+ // with different transform factors
+ {
+ for (size_t k = j; k < num_elements; k += next) // iterations through
+ // pairs within group
+ {
+ match = k + i;
+ fftTransformHelper(data, match, k, product, factor);
+ }
+ factor = mult * factor + factor;
+ }
+ }
+
+ return true;
+}
+
+// Generic template for complex FFT followed by its explicit specializations
+template <class TComplexArray, int NumDims>
+struct CFFT
+{};
+
+// 1D FFT:
+template <class TComplexArray1D>
+struct CFFT<TComplexArray1D,1>
+{
+ // NOTE: passing by pointer is needed to avoid using element access operator
+ static bool FFT_inplace(TComplexArray1D & data, const size_t size,
+ const FFT_direction fft_direction,
+ const char *& error_description)
+ {
+ if(!checkNumElements(size, error_description)) {
+ return false;
+ }
+
+ rearrangeData(data, size);
+
+ if(!makeTransform(data, size, fft_direction, error_description)) {
+ return false;
+ }
+
+ if (FFT_BACKWARD == fft_direction) {
+ scaleValues(data, size);
+ }
+
+ return true;
+ }
+};
+
+// 2D FFT
+template <class TComplexArray2D>
+struct CFFT<TComplexArray2D,2>
+{
+ static bool FFT_inplace(TComplexArray2D & data, const size_t size1, const size_t size2,
+ const FFT_direction fft_direction, const char *& error_description)
+ {
+ int n_rows = static_cast<int>(size1);
+ int n_cols = static_cast<int>(size2);
+
+ // fft for columns
+ std::vector<complex_type> subarray(n_rows); // each column has n_rows elements
+
+ for(int j = 0; j < n_cols; ++j)
+ {
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < n_rows; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ subarray[i] = data[i][j];
+#else
+ subarray[i] = data(i,j);
+#endif
+ }
+
+ if(!CFFT<std::vector<complex_type>,1>::FFT_inplace(subarray, size1,
+ fft_direction,
+ error_description))
+ {
+ return false;
+ }
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < n_rows; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i][j] = subarray[i];
+#else
+ data(i,j) = subarray[i];
+#endif
+ }
+ }
+
+ // fft for rows
+ subarray.resize(n_cols); // each row has n_cols elements
+
+ for(int i = 0; i < n_rows; ++i)
+ {
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int j = 0; j < n_cols; ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ subarray[j] = data[i][j];
+#else
+ subarray[j] = data(i,j);
+#endif
+ }
+
+ if(!CFFT<std::vector<complex_type>,1>::FFT_inplace(subarray, size2,
+ fft_direction,
+ error_description))
+ {
+ return false;
+ }
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int j = 0; j < n_cols; ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i][j] = subarray[j];
+#else
+ data(i,j) = subarray[j];
+#endif
+ }
+ }
+
+ return true;
+ }
+};
+
+// 3D FFT
+template <class TComplexArray3D>
+struct CFFT<TComplexArray3D,3>
+{
+ static bool FFT_inplace(TComplexArray3D & data, const size_t size1, const size_t size2,
+ const size_t size3, const FFT_direction fft_direction,
+ const char *& error_description)
+ {
+ int n_rows = static_cast<int>(size1);
+ int n_cols = static_cast<int>(size2);
+ int n_depth = static_cast<int>(size3);
+
+ std::vector<complex_type> subarray(n_rows); // for fft for columns: each column has n_rows elements
+
+ for(int k = 0; k < n_depth; ++k) // for all depth layers
+ {
+ // fft for columns
+ for(int j = 0; j < n_cols; ++j)
+ {
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < n_rows; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ subarray[i] = data[i][j][k];
+#else
+ subarray[i] = data(i,j,k);
+#endif
+ }
+
+ if(!CFFT<std::vector<complex_type>,1>::FFT_inplace(subarray, size1,
+ fft_direction,
+ error_description))
+ {
+ return false;
+ }
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int i = 0; i < n_rows; ++i) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i][j][k] = subarray[i];
+#else
+ data(i,j,k) = subarray[i];
+#endif
+ }
+ }
+ }
+
+ subarray.resize(n_cols); // for fft for rows: each row has n_cols elements
+
+ for(int k = 0; k < n_depth; ++k) // for all depth layers
+ {
+ // fft for rows
+ for(int i = 0; i < n_rows; ++i)
+ {
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int j = 0; j < n_cols; ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ subarray[j] = data[i][j][k];
+#else
+ subarray[j] = data(i,j,k);
+#endif
+ }
+
+ if(!CFFT<std::vector<complex_type>,1>::FFT_inplace(subarray, size2,
+ fft_direction,
+ error_description))
+ {
+ return false;
+ }
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int j = 0; j < n_cols; ++j) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i][j][k] = subarray[j];
+#else
+ data(i,j,k) = subarray[j];
+#endif
+ }
+ }
+ }
+
+ // fft for depth
+ subarray.resize(n_depth); // each depth strip contains n_depth elements
+
+ for(int i = 0; i < n_rows; ++i) // for all rows layers
+ {
+ for(int j = 0; j < n_cols; ++j) // for all cols layers
+ {
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int k = 0; k < n_depth; ++k) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ subarray[k] = data[i][j][k];
+#else
+ subarray[k] = data(i,j,k);
+#endif
+ }
+
+ if(!CFFT<std::vector<complex_type>,1>::FFT_inplace(subarray, size3,
+ fft_direction,
+ error_description))
+ {
+ return false;
+ }
+
+#ifndef __clang__
+#ifdef __USE_OPENMP
+#pragma omp parallel for
+#endif
+#endif
+ for(int k = 0; k < n_depth; ++k) {
+#ifdef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR
+ data[i][j][k] = subarray[k];
+#else
+ data(i,j,k) = subarray[k];
+#endif
+ }
+ }
+ }
+
+ return true;
+ }
+};
+
+} // namespace impl
+} // namespace simple_fft
+
+#endif // __SIMPLE_FFT__FFT_IMPL_HPP__
--- /dev/null
+/**\r
+ * Copyright (c) 2013-2020 Dmitry Ivanov\r
+ *\r
+ * This file is a part of Simple-FFT project and is distributed under the terms\r
+ * of MIT license: https://opensource.org/licenses/MIT\r
+ */\r
+\r
+// In this file you can alter some settings of the library:\r
+// 1) Specify the desired real and complex types by typedef'ing real_type and complex_type.\r
+// By default real_type is double and complex_type is std::complex<real_type>.\r
+// 2) If the array class uses square brackets for element access operator, define\r
+// the macro __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+\r
+#ifndef __SIMPLE_FFT__FFT_SETTINGS_H__\r
+#define __SIMPLE_FFT__FFT_SETTINGS_H__\r
+\r
+#include <complex>\r
+\r
+typedef float real_type;\r
+typedef std::complex<real_type> complex_type;\r
+\r
+#ifndef __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+#define __USE_SQUARE_BRACKETS_FOR_ELEMENT_ACCESS_OPERATOR\r
+#endif\r
+\r
+#endif // __SIMPLE_FFT__FFT_SETTINGS_H__\r
builddir=debian/build/flavor-
+simple-fft-include=$(CURDIR)/externals/Simple-FFT/include/
export QT_SELECT=qt5
DEB_SRCDIR = .
%:
dh $@ --sourcedirectory=$(DEB_SRCDIR) --buildsystem=meson
+execute_before_dh_auto_build:
+ cp -r $(CURDIR)/debian/missing-sources/Simple-FFT/include/ \
+ $(simple-fft-include)
+
+execute_after_dh_auto_clean:
+ rm -rf $(simple-fft-include)
+
#.PHONY: $(patsubst %,configure_%,$(FLAVORS))
override_dh_auto_configure-arch: $(patsubst %,configure_%,$(FLAVORS))
configure_%:
projects / jacktrip.git / commitdiff