Skip to content

Added 2D Fourier transformation as num.fft2 and num.fft2_inv. #22

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from

Conversation

h4rm
Copy link
Contributor

@h4rm h4rm commented Feb 5, 2014

BTW: Scroll down for updated infos.

Hey Francesco,

I needed some 2D Fourier transformation in my work so I implemented a version. It differs a bit from the 1D calculation, namely I am not using the half complex vectors but I am working with complex numbers from the beginning. The function does 1D transformation for 1 dimension and then 1D transformation for the other dimension. The result is a complex matrix. The result should be comparably fast and probably faster than implementing it on top of fft and fft_inv.
I have added the tests and the docs.
Also I have cross tested the results with numpy just to make sure to have it all right. Numpy also has the 3D Fourier transformation in stock and I will also implement that for my work, following in a later commit.

Cheers,
Ben

@h4rm
Copy link
Contributor Author

h4rm commented Feb 5, 2014

One more thing:

The algorithm work for complex matrices only, so whenever real matrices are given as input I transform the matrix to a complex one with:

matrix.cnew(n1, n2, |i,j| mat:get(i,j))

Maybe there is a better solution.

@h4rm
Copy link
Contributor Author

h4rm commented Feb 5, 2014

BTW, did some performance tests with numpy doing 1e5 2D FFTs on 16x16 matrices.
Python is about 2.5 secs, the shell performs at 4.2 secs. The code is even less performant if you go to higher dimensional matrices (256x256 e.g.).

So I guess performance can still be optimized. Maybe the strides can help me to perform fast transformations.

@h4rm
Copy link
Contributor Author

h4rm commented Feb 5, 2014

So, to be honest, after all I have read about FFT, the GSL version is not the best implementation and FFTW almost always supersedes it.
It is very tempting to go for FFTW with a wrapper for https://github.com/soumith/fftw3-ffi . The FFTW libraries should be installed on every host machine due to GSL dependencies, so this should be no problem, should it?
What do you think?

@soumith
Copy link

soumith commented Feb 5, 2014

might be relevant, though not exactly what you're looking for:
https://github.com/soumith/torch-signal/blob/master/init.lua

@h4rm
Copy link
Contributor Author

h4rm commented Feb 5, 2014

Thanks soumith for the link, it helped me with the ffi.casts.
So having the fftw3-ffi in place, a one dimensional fft is quite simple in gsl-shell:

fftw = require'fftw3.init'

mat = matrix.cnew(10, 1, |i,j|i)
mat2 = matrix.cnew(10,1)

input = ffi.cast("fftw_complex*",mat.data)
output = ffi.cast("fftw_complex*",mat2.data)

plan = fftw.plan_dft_1d(10, input, output, fftw.FORWARD, fftw.MEASURE)

fftw.execute(plan)

print(mat)
print(mat2)

You can also use

plan = fftw.plan_dft_1d(10, input, input, fftw.FORWARD, fftw.MEASURE)

for in-place transformation.

@h4rm
Copy link
Contributor Author

h4rm commented Feb 5, 2014

Okay, the 2D transformation is also extremely easy in comparision to my code:

fftw = require'fftw3.init'

mat = matrix.cnew(4, 4, |i,j|i-1)
mat2 = matrix.cnew(4,4)

input = ffi.cast("fftw_complex*",mat.data)
output = ffi.cast("fftw_complex*",mat2.data)

plan = fftw.plan_dft_2d(4,4, input, output, fftw.FORWARD, fftw.MEASURE)

fftw.execute(plan)

print(mat)
print(mat2)

It produces the same results as numpys.fft2. I am going to prepare an alternative fft implementation for gsl-shell based on the ffi interface to fftw3. I will need it anyway for my projects but we can then discuss if that is a good substitute for the current implementation.

@h4rm
Copy link
Contributor Author

h4rm commented May 27, 2014

This can be cancelled.

franko added a commit that referenced this pull request Jan 20, 2016
fix method call codegen according to lj_parse.c
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants