mborgerding / kissfft

a Fast Fourier Transform (FFT) library that tries to Keep it Simple, Stupid
Other
1.39k stars 278 forks source link

Data stride in N-D FFTs #8

Open GPMueller opened 5 years ago

GPMueller commented 5 years ago

Your library works very nicely for a spatial distribution of 1D (i.e. scalar) data! However, the transformation of 3D (i.e. vector) data can be more tricky, if the data is arranged as (x0,y0,z0, x1,y1,z1, ...).

Maybe I am misunderstanding how the APIs are supposed to be used, but it looks to me like one has to (for N vectors) create 3 arrays (x, y and z) of size N and subsequently do 3 FFTs. This would, in my use case, mean copying data back and forth.

This would be solved by adding a stride to the config or to the APIs. The corresponding initial offset is trivial, as one can simply pass the corresponding data pointer. We tried adding a stride to the APIs, but got stuck in errors in the recursive function calls, which we do not sufficiently understand. It would be great if you could help out with this!

Might also be related to this feature request on source forge, but I don't want to create an account there just for this.

mborgerding commented 5 years ago

The kiss_fft_stride function allow you to grab strided input data. One restriction is that it places each output contiguously. In your case, you'd make 3 calls with fin_stride=3, but the outputs would be packed e.g. fft(x)

MSallermann commented 5 years ago

But kiss_fft_stride only does 1D transforms, right? For some use-cases it would be very nice to also be able to perform multidimensional transforms with an input/output stride. That would allow us to transform over threedimensional arrangements of vectors with much less awkward restrictions on memory layout and improve the usefulness of kiss_fft!

As @GPMueller mentioned we tried adding an input stride to kiss_fft_ndr and it was ... not a nice experience. So if you could give us some pointers that would be very helpful!

mborgerding commented 5 years ago

The basic thing to remember is that multi-dimensional FFTs are simply 1d FFTs done along each dimension. e.g. for a 2d FFT of a matrix X, the fft2(X) is FXF^H, where F represents the DFT operator. Does that help?

related: https://stackoverflow.com/questions/14407522/3d-fft-decomposition-in-2d-fft/14418892#14418892