Solving FFT Problems with MATLAB and FFTW in C

  • Context: MATLAB 
  • Thread starter Thread starter birdhen
  • Start date Start date
Click For Summary
SUMMARY

This discussion addresses issues encountered while calculating the potential of a system from overdensity using Fast Fourier Transforms (FFT) in MATLAB and FFTW in C. The equation phi(k) = -del(k)/(k^2) is central to the problem, particularly when k=0, where division by zero occurs. The recommended solution involves setting up a grid of points that excludes k=0, potentially requiring a phase shift to maintain the integrity of the FFT results. The user reports unexpected outcomes, such as a sine wave instead of the anticipated (1/6)*x^3 when performing the inverse FFT.

PREREQUISITES
  • Understanding of Fast Fourier Transform (FFT) principles
  • Familiarity with MATLAB and FFTW in C
  • Knowledge of Fourier space equations and potential calculations
  • Basic concepts of numerical methods and grid setups
NEXT STEPS
  • Explore techniques for handling k=0 in Fourier transforms
  • Research phase shifting methods in FFT applications
  • Learn about MATLAB's FFT functions and their parameters
  • Investigate common pitfalls in inverse FFT implementations
USEFUL FOR

Researchers and developers working with signal processing, physicists modeling potential fields, and anyone implementing FFT algorithms in MATLAB or C.

birdhen
Messages
31
Reaction score
0
Hi there,

I am having trouble with some Fast Fourier Transforms.
I have been using both MATLAB and FFTW in c.
The problem is that I want to calculate the potential of a system from the overdensity. The equation in Fourier space is simple

phi(k)= -del(k)/(k^2) ,

where del is the overdensity, phi is the potential and k is the wavenumber.

To test this I use a simple del(x) = x, and feed this into the fft. I then run through the k values and divide through by k squared. When k =0, I am not sure what to do as obviously you can't divide by zero.
I have tried a variety of things such as setting phi(k) = 0 when k=0, or just leaving phi(k) = del(k)..both of which I am sure are wrong.

When I reverse FFT this phi(k), my result is always a sine wave with a wavelength the size of the x range, where as I am expecting (1/6)*x^3.

I am having trouble finding the answer online. Maybe I am searching for the wrong thing. Can anyone help?
Best wishes,
 
Physics news on Phys.org
birdhen said:
When k =0, I am not sure what to do as obviously you can't divide by zero.
I have tried a variety of things such as setting phi(k) = 0 when k=0, or just leaving phi(k) = del(k)..both of which I am sure are wrong.
The best choice here is to set up a grid of points such that ##k=0## is not included. This may require a shift that will add a complex phase to the result of the FFT, which can then be removed.
 

Similar threads

  • · Replies 9 ·
Replies
9
Views
3K
  • · Replies 8 ·
Replies
8
Views
3K
  • · Replies 4 ·
Replies
4
Views
3K
  • · Replies 16 ·
Replies
16
Views
3K
Replies
5
Views
3K
  • · Replies 2 ·
Replies
2
Views
4K
  • · Replies 16 ·
Replies
16
Views
15K
  • · Replies 6 ·
Replies
6
Views
2K
  • · Replies 10 ·
Replies
10
Views
3K
  • · Replies 1 ·
Replies
1
Views
3K