# Compute fourrier transform in Julia

**URL:** <https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664>\
**Category:** Numerics\
**Tags:** question\
**Created:** [May 31, 2023, 10:35am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664 "2023-05-31T10:35:16Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 10:35am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/1 "2023-05-31T10:35:16Z")

</div>

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

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

```

how can I compute its fourrier transform

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

```

?

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [May 31, 2023, 10:48am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/2 "2023-05-31T10:48:38Z")

</div>

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](https://juliahub.com/ui/Search?q=Fourier%20transform&type=packages)), 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](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757) has some good tips for asking questions.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:00am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/3 "2023-05-31T11:00:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 11:03am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/4 "2023-05-31T11:03:39Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/d/3/d3b776e141e6e8ab6b629ae35369085acf37a658.png)

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.

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 11:09am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/5 "2023-05-31T11:09:07Z")

</div>

> [@kellertuer](#):
>
> Could you be a bit more precise which Fourier transform you mean and in what sense you mean “numerically”?

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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:10am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/6 "2023-05-31T11:10:00Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:13am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/7 "2023-05-31T11:13:15Z")

</div>

Ah found it, for the Fourier transform see here

> **[Fourier transform | Other conventions](https://en.wikipedia.org/wiki/Fourier_transform#Other_conventions)**
>
> The Fourier transform can also be written in terms of angular frequency:

where your factor mit come from and the same for the DFT/FFT see e.g. [DFT matrix - Wikipedia](https://en.wikipedia.org/wiki/DFT_matrix)

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 11:20am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/8 "2023-05-31T11:20:56Z")

</div>

> [@kellertuer](#):
>
> 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).

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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:25am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/9 "2023-05-31T11:25:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:31am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/10 "2023-05-31T11:31:07Z")

</div>

> [@Fourier](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 11:52am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/11 "2023-05-31T11:52:46Z")

</div>

I think you are actually just missing an `ifftshift`. See

```plaintext
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](https://global.discourse-cdn.com/julialang/original/3X/9/0/909efc3648c5dfb894f6c200255d6d1aba574848.png)

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 🙂

_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`).

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 8:25pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/12 "2023-05-31T20:25:44Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/7/4/747ade3b020f4a0de57f3c0c245fff61397ca566.png)

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 😅 and somehow I cant find much information regarding fft vs actual fourrier transform on the internet.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 8:39pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/13 "2023-05-31T20:39:46Z")

</div>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/c/0/c09a867d30e3d05e826178f86e520f50f4d19b0b.png)

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](https://link.springer.com/book/10.1007/978-3-030-04306-3) (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).

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 8:58pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/14 "2023-05-31T20:58:33Z")

</div>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/d/0/d0d13206095c9db2c8db00cfa3b4bc4f299f1a3c.png)

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.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [May 31, 2023, 9:05pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/15 "2023-05-31T21:05:17Z")

</div>

> [@Fourier](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 9:10pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/16 "2023-05-31T21:10:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [May 31, 2023, 9:21pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/17 "2023-05-31T21:21:46Z")

</div>

Hey, yes, thats maybe a good starting point, thanks for the idea. 😄  
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.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [May 31, 2023, 9:32pm UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/18 "2023-05-31T21:32:46Z")

</div>

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 😉

---

<div class="post-metadata">

**Author:** ![TheLateKronos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thelatekronos/32/12824_2.png) [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)\
**Post date:** [June 3, 2023, 11:27am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/19 "2023-06-03T11:27:30Z")

</div>

I made a package to automatically perform this normalization and shifting, which you may find interesting. It is called [EasyFFTs.jl](https://github.com/KronosTheLate/EasyFFTs.jl).

---

<div class="post-metadata">

**Author:** ![Fourier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fourier/32/38176_2.png) [@Fourier](https://discourse.julialang.org/u/Fourier)\
**Post date:** [June 5, 2023, 8:34am UTC](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664/20 "2023-06-05T08:34:24Z")

</div>

Oh, very nice! I will take a look.

[Next page](https://discourse.julialang.org/t/compute-fourrier-transform-in-julia/99664.md?page=2)