xtensor-stack / xtensor-fftw Goto Github PK

Greeting,
I wonder this package is not supporting the m1 chip now? Do you have any plans to implement?

Best,Fu

Hey @egpbos how about a new xtensor-fftw release?

I am happy to make the tag and update the recipes if you are ok with it.

Real FFTs on arrays where the last number in the shape is an odd number require special care. In particular: an inverse rFFT cannot automatically determine whether the output should be an odd or even sized (real) array, so the user has to supply this information. I added a flag for this in the relevant functions, but this behavior is undocumented. Add to docs.

Data arrays in FFTW must be in row major format (http://www.fftw.org/fftw3_doc/Row_002dmajor-Format.html). Currently, xtensor-fftw defines output arrays as row major, but it does not check for row major-ness on the input arrays. Should this be done or do we leave it to the user?

When following the install instructions: https://github.com/xtensor-stack/xtensor-fftw#installation

The required headers do not get installed, e.g. basic_float.hpp, which is included from basic.hpp does not get installed.

This is, I believe, because ${XTENSOR_FFTW_HEADERS} is defined here without those headers:

https://github.com/xtensor-stack/xtensor-fftw/blob/master/CMakeLists.txt#L179

It's now included implicitly only if it happens to be in the same directory as xtensor... which it often does, but not necessarily. See xtensor CMakeLists.

Right now, due to my inexperience with templates, I have implemented the numpy-like interface functions using a number of (inline) functions that call a template instance, e.g.:

inline xt::xarray<std::complex<float> > fft (const xt::xarray<std::complex<float> > &input) {
  return _fft_<std::complex<float>, std::complex<float>, 1, true, false, false, fftwf_plan_dft_1d, fftwf_execute, fftwf_destroy_plan> (input);
}

I'm not sure, but it seems likely that the compiler will simply compile all these functions when a user includes this header. This makes the resulting binary a lot larger than it needs to be and increases compile time.

I think I now know how to rewrite the _fft_ and _ifft_ templates in such a way that I can write using template aliases for the three families (fft, rfft and hfft) that only depend on input precision. This will involve the following modifications to _fft_ and _ifft_:

• Putting the fftw plan creation functions (constexpr static) in templates with parameters for precision, dimensionality, "c2c vs r2c" and direction.
• Using the fftw execute and destroy_plan functions that I put in the precision dependent fftw_t template as well (these are also constexpr static).
• Putting the output type in -> notation, using decltype or declexpr to determine the type based on template parameters for input precision, "c2c vs r2c" and direction.

The latter point was the crucial knowledge gap that caused me to write the functions as I did, since I didn't get auto-type-deduction working for the output type.

Currently, the helper functions in helper.hpp only work on single dimensional arrays. We should extend this to multi-dimensional arrays. numpy.fft's helper functions use the axes keyword argument for this. This syntax could also be used in broadcasting operations (#4). Possibly, some other syntax is more appropriate in the C++ / xtensor context. We should investigate.

Like these:

• https://dsp.stackexchange.com/questions/633/what-data-should-i-use-to-test-an-fft-implementation-and-what-accuracy-should-i
• http://www.sccon.ca/sccon/fft/fft3.htm

A cling Jupyter notebook showing interactive usage. Something with an image, something with filtering, maybe something with audio. Maybe there are numpy.fft notebooks that we could draw inspiration from as well.

debug-8-7298.csv

On my machine (Rocky Linux 8.9; FFTW 3.3.5), when loading this data and execute FFTW, it reports a smallbin double linked list corruption.

The example is as following (C++ code):

fstream fs_debug;
fs_debug.open("/home/winfred/l1b-pro/data/output/debug-8-7298.csv", ios::in);
xtensor<float, 2> delta_sigx_han_t = xt::load_csv<float>(fs_debug);
xtensor<float, 1> delta_sigx_han_tt = xt::zeros<float>({delta_sigx_han_t.size()});
for (int i = 0; i < delta_sigx_han_t.size(); i++)
{
    delta_sigx_han_tt(i) = delta_sigx_han_t(0, i);
}
auto pdsig_t = xt::fftw::rfft(delta_sigx_han_tt);

FFTW is not thread safe, except for fftw_execute.
How about using c++11 locks around calls to fftw_plan_dft and fftw_destroy_plan for thread safety?

Thank you for developing xtensor-fftw.
I would like to request the following extension that will be improve its usability.

It will be better if there is a broadcasting in xtensor-fftw.
I'm studying on multi-channel acoustic signal processing, and often have to conduct FFTs to microphones.
In numpy, we can broadcast FFT like this way

np.fft.rfft(np.ndarray([M, T]), axis=1).shape

In my case M and T indicates the number of microphones and time frames.
This will be simplify codes of such signal processings.

Best,

Take inspiration from the other xtensor pages, e.g. xensor-io has a nice compact one, that would be good as a header. The rest can be something like:

1. Intro
2. Example
3. Installation
4. Dependencies
5. Usage (i.e. more examples, or links to notebooks/binder)
6. Building and running the tests
7. License

CMake fails to find FFTW (version 3.3.9) which was installed from source. I use CMake version 3.19.3 on Ubuntu 20.04.

$ cmake --version
cmake version 3.19.3

CMake suite maintained and supported by Kitware (kitware.com/cmake).

FFTW is installed in /usr/local:

/usr/local/lib/libfftw3q.a
/usr/local/lib/libfftw3q.a
/usr/local/include/fftw3.h
/usr/local/include/fftw3.f
/usr/local/include/fftw3l.f03
/usr/local/include/fftw3q.f03
/usr/local/include/fftw3.f03
/usr/local/lib/pkgconfig/fftw3q.pc
/usr/local/lib/cmake/fftw3q/FFTW3qConfig.cmake
/usr/local/lib/cmake/fftw3q/FFTW3qConfigVersion.cmake
/usr/local/lib/cmake/fftw3q/FFTW3LibraryDepends.cmake
/usr/local/lib/cmake/fftw3q/FFTW3LibraryDepends-release.cmake

I suspect that your CMake code is no longer compatible with current CMake versions.

Simple cmake call:

$ cmake ../xtensor-fftw-0.2.6
-- The C compiler identification is GNU 9.3.0
-- The CXX compiler identification is GNU 9.3.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /usr/bin/cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Building xtensor-fftw v0.2.6
-- Performing Test arch_native_supported
-- Performing Test arch_native_supported - Success
-- Found xtensor: /usr/local/include/xtensor
-- Found xtl: /usr/local/include/xtl
-- Found PkgConfig: /usr/bin/pkg-config (found version "0.29.1") 
CMake Error at /snap/cmake/769/share/cmake-3.19/Modules/FindPackageHandleStandardArgs.cmake:218 (message):
  Could NOT find FFTW (missing: FLOAT_LIB DOUBLE_LIB LONGDOUBLE_LIB)
Call Stack (most recent call first):
  /snap/cmake/769/share/cmake-3.19/Modules/FindPackageHandleStandardArgs.cmake:582 (_FPHSA_FAILURE_MESSAGE)
  cmake/Modules/findFFTW/FindFFTW.cmake:277 (find_package_handle_standard_args)
  CMakeLists.txt:88 (find_package)


-- Configuring incomplete, errors occurred!
See also "/home/bebuch/media/tmp/install/xtensor-fftw-0.2.6-build/CMakeFiles/CMakeOutput.log".

Call with FFTW_ROOT:

$ cmake -DFFTW_ROOT=/usr/local ../xtensor-fftw-${XTENSOR_FFTW_VERSION}
-- The C compiler identification is GNU 9.3.0
-- The CXX compiler identification is GNU 9.3.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /usr/bin/cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Building xtensor-fftw v0.2.6
-- Performing Test arch_native_supported
-- Performing Test arch_native_supported - Success
-- Found xtensor: /usr/local/include/xtensor
-- Found xtl: /usr/local/include/xtl
CMake Warning (dev) at CMakeLists.txt:88 (find_package):
  Policy CMP0074 is not set: find_package uses <PackageName>_ROOT variables.
  Run "cmake --help-policy CMP0074" for policy details.  Use the cmake_policy
  command to set the policy and suppress this warning.

  CMake variable FFTW_ROOT is set to:

    /usr/local

  For compatibility, CMake is ignoring the variable.
This warning is for project developers.  Use -Wno-dev to suppress it.

-- Found PkgConfig: /usr/bin/pkg-config (found version "0.29.1") 
CMake Error at /snap/cmake/769/share/cmake-3.19/Modules/FindPackageHandleStandardArgs.cmake:218 (message):
  Could NOT find FFTW (missing: FLOAT_LIB DOUBLE_LIB LONGDOUBLE_LIB)
Call Stack (most recent call first):
  /snap/cmake/769/share/cmake-3.19/Modules/FindPackageHandleStandardArgs.cmake:582 (_FPHSA_FAILURE_MESSAGE)
  cmake/Modules/findFFTW/FindFFTW.cmake:277 (find_package_handle_standard_args)
  CMakeLists.txt:88 (find_package)


-- Configuring incomplete, errors occurred!
See also "/home/bebuch/media/tmp/install/xtensor-fftw-0.2.6-build/CMakeFiles/CMakeOutput.log".

We want to implement a very simple FFT interface, mimicking the numpy.fft default behavior. This means that we want the following matrix of functions:

• fft for standard complex to complex transforms, rfft for real FFTs (real signal) and hfft for Hermitian FFTs (real Fourier spectrum);
• an inverse for each operation, prepending an i, e.g. ifft;
• a 1D (default), 2D, 3D and nD version of each function, the latter three indicated by appending 2, 3 or n, e.g. irfft3.

That's a total of 24 functions.

Additionally, we want to support single, double and long double precision arrays. This can be implemented through template specializations. The total amount of implemented functions would then be 72.

Also part of this interface would be broadcasting over multiple array dimensions, if the dimension of the function is smaller than that of the array itself; see #4.

Following #57 and #58, we should add tests that use the specific type headers basic_float.hpp etc.

The test on my last PR failed in a way that does not match the PR. Added to that, the error is the same than the last failed run on master.

My PR failure:

• https://travis-ci.org/xtensor-stack/xtensor-fftw/jobs/608699693?utm_medium=notification&utm_source=github_status

Last build on master:

• https://travis-ci.org/xtensor-stack/xtensor-fftw/builds/600631847

I may fix it if you dont before tomorrow (note Im in the US west coast so my schedule does not overlap with the main developer here)

We should benchmark whether the 1-3 dimensional special functions are actually faster than the n-dimensional one. If not, we can just use one implementation function for all ranks, with simple wrappers for the 1, 2 and 3 dimensional cases. This would reduce the number of implemented functions from 72 to 18, plus 54 simple wrapper functions.

Following #28, to make things even more general, replace shape type with xarray<T, L>::shape_type when you get the shape of an array, thus the code is more flexible and less sensitive to changes of shape / strides types in xtensor (suggestion by @JohanMabille ).

Follow some good suggestions from https://youtu.be/bsXLMQ6WgIk.

• "Use add_library for installation of header only libraries", which xtensor-fftw is (xtensor itself as well).
• Set CXX_CLANG_TIDY=ON by default.

Nevermind, misinterpreted some stuff.

adding #define XTENSOR_USE_XSIMD results in the errors below:

Configuration:
C++20
VS2019
Xtensor 0.24.0
Xsimd 8.0.3

As in https://github.com/QuantStack/xtensor-io/blob/a57e3b90ca1025a180c5813c56b8a344ec2d8165/include/xtensor-io/xaudio.hpp#L20-L21, we can put the library loading step in the headers themselves, which saves cling/notebook users from doing it.

After implementing the basic numpy.fft-like interface (#6), we should benchmark it against the actual numpy.fft functions.

We use template specializations to make sure we get the correct precision. The problem, for instance, with a non-template f(const xt::xarray<float> &input) is that it will also compile when you pass a xt::xarray<double>. Passing by const-reference in this sense behaves similarly to passing by value; it triggers the creation of a temporary variable -- input in this case -- though in the case of a reference the data is not actually copied.

Using = delete syntax on the default template, e.g.:

template<typename real_t>
xt::xarray< std::complex<real_t> > fft(const xt::xarray<real_t> &input) = delete;

we can make sure that calls to non-implemented specializations don't compile. If this is left out, the compilation will succeed, but the linker will fail, and this gives less informative error messages.

Unfortunately, clang has a bug that prevents us from using the nice = delete syntax. For the moment we're using a workaround:

template<typename real_t>
xt::xarray< std::complex<real_t> > fft(const xt::xarray<real_t> &input) {
   static_assert(sizeof(real_t) == 0, "Only specializations of fft can be used");
};

This is ugly. Once the bug in clang is fixed, we will want to replace this with the = delete syntax.

Look into which license xtensor-fftw is exactly supposed/allowed to use for what (code, conda package), given that FFTW3 itself is under GPL.

The simple numpy-style interface (#6) will by default hide some of FFTW's features, opting for simplicity and default settings for decent performance. We will also implement a high performance interface to access functions like

• Choosing FFTW_MEASURE, FFTW_PATIENT, or FFTW_EXHAUSTIVE planner flags (instead of the default FFTW_ESTIMATE);
• Reusing plans and wisdom;
• Using fftw_malloc to guarantee optimal memory alignment;
• FFTW_PRESERVE_INPUT vs FFTW_DESTROY_INPUT, the latter can sometimes be faster.

In basic.hpp, there are a few exceptions to safeguard against nullptrs from FFTW plan creation, which for instance can happen if you try to build a plan for a multi-dimensional real FFT with a FFTW_PRESERVE_INPUT flag, because such algorithms are not available in FFTW. These exceptions are there more as a developer safeguard than a user safeguard. Eventually, it would be nicer if such "illegal" combinations of plan parameters were disallowed by the code itself, e.g. by not allowing users to manually specify planner flags, but to hard code all legal possibilities. This will eliminate the need for the exceptions.

Therefore, note to self: don't try to get 100% coverage by making a test with illegal input, or something like that. Just fix the API so that the exceptions can be removed.

The current setup is to use Travis for Linux and macOS and Appveyor for Windows. We should replace both by just GitHub Actions to keep things simple. Also, Travis doesn't actually work anymore, so Actions would fix that as well.

Hi, I'm adding xtensor-fftw to our build process at work. We don't use FFTW, but use a compatible implementation through MKL. The CMakeLists file however doesn't check for this. Could you please add a check for MKL? Thanks!

Use e.g. licecap. Put it on Twitter and maybe on the documentation page under Interactive examples.

I am currently using xtensor as well as xtensor-blas. To use xtensor-fftw, Conda asks me to downgrade these packages:

  Package         Version  Build       Channel           Size
───────────────────────────────────────────────────────────────
  Install:
───────────────────────────────────────────────────────────────

  + xtensor-fftw    0.2.6  h4bd325d_1  conda-forge       12kB

  Downgrade:
───────────────────────────────────────────────────────────────

  - xtl             0.7.5  hf52228f_0  conda-forge     Cached
  + xtl            0.6.23  h4bd325d_1  conda-forge       98kB
  - xtensor        0.24.7  h00ab1b0_0  conda-forge     Cached
  + xtensor       0.21.10  h4bd325d_0  conda-forge      182kB
  - xtensor-blas   0.20.0  h4bd325d_0  conda-forge     Cached
  + xtensor-blas   0.17.2  h998ca80_1  conda-forge      264kB

  Summary:

  Install: 1 packages
  Downgrade: 3 packages

Is it possible to use xtensor-fftw with recent xtensor versions or are there some issues?

We should add items to the numpy to xtensor cheat sheet and a code snippet in the xtensor related projects section.

Some decent machine + maybe a website to publish the recorded results. Continuous benchmarks would be even better, since then we can catch performance regressions more easily.

Hi,

First of all, thank you for developing both xtensor and xtensor-fftw! I encounter an issue while trying to integrate xtensor-fftw with my current project using xtensor (for which I install both xtensor and xtensor-fftw through conda install -c conda-forge xtensor xsimd xtensor-fftw) --

Here is my CMakeLists.txt file --

cmake_minimum_required(VERSION 3.1)
project(example)

find_package(xtl REQUIRED)
find_package(xtensor REQUIRED)
find_package(xtensor-fftw REQUIRED)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED TRUE)
set(CMAKE_CXX_EXTENSIONS TRUE)

add_executable(main src/main.cpp)
target_link_libraries(main xtensor xtensor::optimize xtensor::use_xsimd)

And, the main.cpp looks like this:

#include <complex>
#include <iostream>

#include <xtensor-fftw/basic.hpp>
#include <xtensor-fftw/helper.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xmath.hpp>

int main(int argc, char *argv[]) {
  using T = std::complex<double>;
  auto size = 128;
  xt::xarray<T> arr = xt::zeros<T>({size, 1});

  for (auto i = 0; i < 10; i++) {
    arr(i, 0) = arr(size - i, 0) = 1.0;
  }
  auto arrfft = xt::fftw::fft(arr);
  std::cout << arrfft;
  return 0;
}

However, the linker complaints

✖  make                                                                                                                                                                                                                                                                                     (xt) 
Consolidate compiler generated dependencies of target main
[ 50%] Building CXX object CMakeFiles/main.dir/src/main.cpp.o
[100%] Linking CXX executable main
Undefined symbols for architecture x86_64:
  "_fftw_destroy_plan", referenced from:
      xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag> xt::fftw::_fft_<std::__1::complex<double>, std::__1::complex<double>, 1ul, -1, true, false, false, fftw_plan_dft_1d, fftw_execute, fftw_destroy_plan, void>(xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag> const&) in main.cpp.o
  "_fftw_execute", referenced from:
      xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag> xt::fftw::_fft_<std::__1::complex<double>, std::__1::complex<double>, 1ul, -1, true, false, false, fftw_plan_dft_1d, fftw_execute, fftw_destroy_plan, void>(xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag> const&) in main.cpp.o
  "_fftw_plan_dft_1d", referenced from:
      std::__1::enable_if<(dimensional::is_1<1ul, true>::value) && ((-1) != (0)), xt::fftw::fftw_t<std::__1::complex<double> >::plan>::type xt::fftw::fftw_plan_dft_caller<1ul, -1, true, std::__1::complex<double>, std::__1::complex<double>, fftw_plan_dft_1d, false, false>(xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag> const&, xt::xarray_container<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >, (xt::layout_type)1, xt::svector<xt::uvector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type, 4ul, std::__1::allocator<std::__1::vector<std::__1::complex<double>, xsimd::aligned_allocator<std::__1::complex<double>, 32ul> >::size_type>, true>, xt::xtensor_expression_tag>&, unsigned int, bool) in main.cpp.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make[2]: *** [main] Error 1
make[1]: *** [CMakeFiles/main.dir/all] Error 2
make: *** [all] Error 2

I'm wondering could you please help me with this?

The CMakeFiles have become a bit messy and seem overcomplicated. I think a modern CMake makeover would help a lot. Things like FFTW_LINK_FLAGS, for instance, should not be used as such, since linker flags and linking targets are two things that are defined separately, so combining them in one variable is not CMake-ish. Other things like set for custom variables are also code smells.

The default behavior for inverse real FFT functions (c2r) in FFTW is to destroy the input arrays during the calculation. This is probably very unexpected for most users.

We should add this to the documentation, probably multiple times (in the general howtos, and in the function docs). Maybe we should even print run-time warnings.

I don't know if it makes sense and if it's easy to achieve it, but it would be nice to not have to evaluate the xarrays for a call of fft function, meaning that the implementation of the functions would be:

template <class T>
inline xt::xarray<std::complex<float> > fft3 (const xt::xexpression<T> &input) {
      return ...;
    }

instead of

inline xt::xarray<std::complex<float> > fft3 (const xt::xarray<std::complex<float> > &input) {
      return ...;
    }

And maybe the return type could be an xexpression too..