+++ /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
projects / jacktrip.git / commitdiff