Compute fourrier transform in Julia

Specific Domains Numerics
1

Given a function evaluated to a grid, how can I compute its fourrier transform nummericially?
For example, given the box function

function step(x)
if abs(x) < a
    return 1
else 
    return 0
end
end

how can I compute its fourrier transform

function stepff(p)
return 2*sin(a*p)/p
end

?

2

I don’t know the answer to your question, but both of these functions reference an a variable that isn’t passed as an argument. You typically want to avoid referencing globals in Julia (unless a is a const).

Also, forgive me, but this reads a bit like a homework question. There are surely people here that can help (have you seen FFTW.jl? Or JuliaHub), but I think being a bit more specific with what you’re asking might help. Eg, what is your input (using a MWE)? What have you tried? What is the error you’re seeing?
This has some good tips for asking questions.

1 Like
3

Could you be a bit more precise which Fourier transform you mean and in what sense you mean “numerically”?

  • sampling step, you could do a DFT (FFT)
  • if you continue the periodic continuation of step you could ask whether there is a computational way to return its Fourier coefficients (though that will probably still include some integration
  • If keep step, you could askl the Fourier transform of step (kind of like Zygote computes derivatives), i.e. in a symbolic sense

For 2 and 3 I do not have a good idea what to do in Julia here. But for the first, sure, FFTW.jl is probably your best choice.

1 Like
4

Hi, sure I can provide more detail. a can be any factor, but regarding the implementation I used a=1.

As to what I have tried, I wanted to offrain from using too much detail because it will streamline the discussion into a certain direction, eg using fftw, which implements the discrete fourrier transform and I am not sure about the relationship between discrete fourrier transform and the actual fourrier transform.

I tried the following:

I tried using julia for the computation and doing some plots

x = -5:dx:5
p = fftfreq(length(x),2pi/dx)
scatter(p, real.(fft(step.(x))))
scatter!(p, imag.(fft(step.(x))))

however the real and imaginary parts dont really seem to match the expected result at all.

However the absolute value only seems to be off by a constant

plotfft

Also when looking at the formula, the discrete fourrier transform is given by

X_k=\sum_{n=0}^{N-1} x_n \cdot e^{-\frac{i 2 \pi}{N} k n}

while the fourrier transform is given by

\hat{f}(\xi)=\int_{-\infty}^{\infty} f(x) e^{-i 2 \pi \xi x} d x

I am however unsure how they are related.
I do see some relationship between the two, in particular if I would to approximate the integral by a Riemann sum

\int_{-\infty}^{\infty} f(x) \, dx \approx \frac{1}{N}\sum_{N}f(x[N])

this reletionship would suggest divideing by N to obtain the correct fourrier transform but that does’t lead to a correct result eather, also it doesn’t help with the feact that the real and complex components don’t match.

I also looked at [this post][2] from maths stackexchange, but it did not really prove to be useful, because the answers dont answer my question.

5

What I want to do is the following: Given a function f evaluated on a grid x, how can I obtain the fourrier transform \hat f? Where the fourrier transform is given matematicially by the equation in my post above. So reformulating the question you could ask something like this: What is the nummerical equivalent to performing a fourrier transform? Because just applying the fft and using fftfreq doesnt seem to get the desired result.

6

Keep in mind that the 1/N is placed on different places depending on who you ask (I for example prefer a 1/sqrt(N) upfront both the DFT and its inverse).
The same discussion with the factor upfront holds for the Fourier transform on the real line, which I usually know introduced with a \frac{1}{2\pi} upfront compared to your notation.

Not also that the DFT stems from periodisation of your original function, and then discretising the Fourier coefficients, and not necessarily from the Fourier integral and a Riemann sum as you wrote (I am not sure where to stop said sum for example).

7

Ah found it, for the Fourier transform see here

where your factor mit come from and the same for the DFT/FFT see e.g. DFT matrix - Wikipedia

8

Yes, that makes sense. Is there any relationship between the fft and the fourrier transform I would do in maths? Or is the DFT actually more related to a fourrier series?

I dont perticularily care about convention as of now, because I can just choose some parameters accordingly. The fundamental problem is that the real and imaginary part of the fft dont match the values obtained anaticially at all (except for the magnitude), for example the analytical solution says there should not be an imaginary component, but the value obtained using fft does.

9

The main steps are

  1. you restrict your function on the real line f to some interval (maybe [-\pi,\pi)) and periodise. Let’s call the result g
  2. For periodic functions you have the Fourier series with Fourier coefficients c_k(g)
  3. Discretising c_k(g) e.g. with trapezoidal rule yields the DFT (FFT is a speedup in computations for this

So if you are off by a factor, that is just the conventions you miss.
The complex part might come from the fftshift-problem.

10

I think the magnitude is a problem of scaling from the convention, the imaginary part probably from the fftshift. Besides that I am not sure what “at all” means here.

11

I think you are actually just missing an ifftshift. See

step(x) = abs(x) < 1 ? 1.0 : 0.0
x = range(-5.0,5.0,101)
y = step.(x)
Y = fft(ifftshift(y)) # image part zero up to machine precision
plot(-50:50,imag.(Y)); plot!(-50:50, real.(Y)) #I was a bit lazy corresponding to your p

yields

Screenshot 2023-05-31 at 13.52.15
Screenshot 2023-05-31 at 13.52.151180Ă—778 36.4 KB

where the -50,50 corresponds a bit more to the c_k(g) (their index k to be exact), since I usually do more periodic stuff and more used to that.

To be more precise: the sampling for the FFT/DFT should be starting at x=0, the ifftshift does exactly that. then all fits :slight_smile:

edit: I am also not 100% sure your sinc is correct (stepff) since I would have thought the 2 is too much and it should just be sin(x)/x (and of course 1 at x=0).

4 Likes
12

Thanks a lot that looks much better, but not quite like a \frac{\sin(x)}{x} function, because that function has a global maximum around 0, while your function does not.

Also here can you see the derivation that the fourrier transform of the box is really the \mathrm{sinc} function, the factor of two comes from the definition of the \sin in terms of complex exponential. In the bottom right you can see what I expect the fourrier transform to look like. Maybe the fft needs to be shifted as well.

Also I would like to state that calculating p (or whatever you want to call the coordinates in fourrier space) is quite important to me, because I want to integrate over them and apply fourrier theorems like for the derivative.

image
image1449Ă—943 87.6 KB

I am very sorry that I bother you with such things that should be trivial, but I never had a course in it and now I need it for my thesis which is due by end of the week :sweat_smile: and somehow I cant find much information regarding fft vs actual fourrier transform on the internet.

13

that depends a lot where you want to see the zero. You just have to (again) ffthisft mine and there you are, or in other words:

For the FFT (in time and domain) zero is always on the very left.

For the factor 2 – that again depends a lot which Fourier transform definition you take.

So, my example bringing you the “bump” into the center – note the outer fftshift

using Plots, FFTW
step(x) = abs(x) < 1 ? 1.0 : 0.0
x = range(-5.0,5.0,101)
y = step.(x)
Y = fftshift(fft(ifftshift(y))) # image part zero up to machine precision
plot(-50:50,imag.(Y)); plot!(-50:50, real.(Y))

yields

Bildschirmfoto 2023-05-31 um 22.38.55
Bildschirmfoto 2023-05-31 um 22.38.551178Ă—790 37.1 KB

The main things to remember with fftshift (starting continuous again) its moving between [0,2\pi) and [-\pi,\pi), and there is two variants (shift and its inverse), because for an odd number of values one is exactly a single index shift more than the other.

edit: Well for a thesis due end of this week, that is a bit late. You will have to find out about the p yourself, because I am not sure what your interval is (was it -5 to 5?)not what your sampling rate is, nor your actual width a.

My recommended reference is Numerical Fourier Analysis | SpringerLink (I might be biased, the third author is my former work group leader for my habilitation – and to some extend I worked in a similar area for some time).

1 Like
14

Okay, so if I normalize by sqrt(N), I can get something that looks more or less as expected.

dx = 0.1
step(x) = abs(x) < 1 ? 1.0 : 0.0
x = -5:dx:5
p = fftfreq(length(x), 2pi/dx) |> ifftshift # 101 is odd
y = step.(x)
Y = fft(ifftshift(y))./sqrt(length(p)) |> ifftshift # image part zero up to machine precision
plot(p,imag.(Y)); plot!(p, real.(Y))

plot!(p, 2*sin.(p)./p, linestyle=:dash, linewidth=2, color=:red)
savefig("fft.png")

fft

PS: Dont worry I dont want to calculate some box functions for my thesis, this is just a toy problem. I have some result for gauss orbitals that was derived analyticially and I want to check if it is actually correct by performing the fourrier transform. Thats also why the normalisation factor is important. But its going to be allright : D.

15

Yes, they are different things indeed, and with the DFT you obtain the coefficients in a Fourier series expansion of a function.

Numerical Fourier Transform is a more difficult problem, as esencially you have to do numerical integration, and if you need high frequencies your integral can become highly oscillatory.

If your thesis is due to next week or so, I’d advise you truncate the integration domain to a large number, and compute the integral numerically with QuadGK, and just wait for the computation to finish.

If you need more advanced fast methods that exploit the FFT, it is possible to obtain them, but you’ll need to do some extra work. Basically, you can apply a DFT to obtain a Fourier series of a function (making sure that you are under the conditions for this series to converge fast, which is non-trivial) and then you can integrate each term analytically. In this way you can go from a Fourier series to a Fourier transform.

1 Like
16

Uff, that looks for me like an off-by-one error as close as they both are. As an idea: if you have N points you only have N-1 intervals between them, in the same idea 2\pi / dx and 2pi*10*N (if the late our does not mess up my ±1) are off by one in that sense, and I think either your p or your scaling of the sink – not in height but in x are slightly off.

But it has been a long day and these ±1 are then a bit hard to track down that late (since the screenshot was in German: ein bischen tüdelig dat Ganze).

edit: First thing to note – your Y has a double ifftshift, I think the second one (back) should be fftshift.

1 Like
17

Hey, yes, thats maybe a good starting point, thanks for the idea. :smile:
Also thanks for the book recomendation. I guess I am going to read it a bit, because knowing how fft works in this day and age is relevant I guess. If I do find a final answer, Ill post it here.

1 Like
18

Good luck. A final idea – since everything is periodic, you have to also bet careful not to sample a point twice.
For example, if you are 2\pi-periodic, you make a mistake if you e.g. sample at x=\pi and x=-\pi. In your case, if you are 10-periodic (for the DFT/discrete case, not the your original rect/sinc in the continuous case) then ±5 is the same point, so you can not have that twice in your samples. Also one of these technicalities with an off-by-one error :wink:

19

I made a package to automatically perform this normalization and shifting, which you may find interesting. It is called EasyFFTs.jl.

1 Like
20

Oh, very nice! I will take a look.

1 Like