# The need of (i)fftshift's to compute a PDE using Fourier Transform

**URL:** <https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587>\
**Category:** Numerics\
**Tags:** question, fftw, diffeq, pde, dsp\
**Created:** [March 5, 2021, 9:05pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587 "2021-03-05T21:05:00Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 5, 2021, 9:05pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/1 "2021-03-05T21:05:00Z")

</div>

Hello,  
I’m currently writing an algorithm to solve a specific kind of PDE, via Fourier Transform.  
In the several examples that I saw, the general idea is to transform everything to the frequency domain using an `fft`, do the computational stuff (Forward Euler, etc) and in the end transform it back to the natural domain using `ifft`, thus obtaining the solution to our PDE.

My difficulties come when I, naively, apply the `fft` and `ifft` to my equation. For a toy example, where the analytical solution is known I get the following:

 ![sol1derrada](https://global.discourse-cdn.com/julialang/original/3X/6/8/68c8d4eef74e3e3b991aed4980fb3d3ad7b93b29.png)  
It seems to me that some kind of boundary conditions is being wrongly satisfied. The solution doesn’t have the domain periodicity (in this case [-10;10]).

However, if I replace `fft` and `ifft` by `fftshift(fft(ifftshift(U)))` and `fftshift(ifft(ifftshift(U)))`, respectively. I obtain the correct answer:

 ![sol1d](https://global.discourse-cdn.com/julialang/original/3X/8/f/8f14a0b0fd27be7a5987f6ec42c1c9fcbfda8986.png)

Well, this seems strange to me, because the `(i)fftshifts` commands are only to shift the Fourier coefficients to their regular order when transforming to the frequency domain. But this is  
reestablished when we perform the inverse transform:

```julia
using FFTW
julia> A=trunc.(Int64,rand(6)*10)
6-element Array{Int64,1}:
 0
 2
 3
 7
 1
 2

julia> â=fft(A) # F. coefs with shifted order
6-element Array{Complex{Float64},1}:
 15.0 + 0.0im
 -7.0 - 1.7320508075688772im
  3.0 + 1.7320508075688772im
 -7.0 + 0.0im
  3.0 - 1.7320508075688772im
 -7.0 + 1.7320508075688772im

julia> ifft(â) # equals to original vector A
6-element Array{Complex{Float64},1}:
 0.0 + 0.0im
 2.0 + 0.0im
 3.0 + 0.0im
 7.0 + 0.0im
 1.0 + 0.0im
 2.0 + 0.0im

```

Someone know why this behavior happens when solving PDE using the Fourier Transform?  
I really need that my algorithm be fast as possible. For each time step computing 2 shifts for every time I compute an `fft` it’s too much computations spent on this issue, I guess 😕

It’s the first time that I’m using this method to solve PDEs, so I don’t know that much about signal processing and fourier transforms. Any help?

Best regards,  
Tiago

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 5, 2021, 9:16pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/2 "2021-03-05T21:16:28Z")

</div>

> [@tiago](#):
>
> However, if I replace `fft` and `ifft` by `fftshift(fft(ifftshift(U)))` and `fftshift(ifft(ifftshift(U)))` , respectively. I obtain the correct answer:

Probably that means that you are taking the derivative incorrectly. See [these notes](https://math.mit.edu/~stevenj/fft-deriv.pdf) (which use 0-based indexing, so you will need to subtract 1 from Julia’s usual 1-based indices to convert).

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 6, 2021, 12:24pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/3 "2021-03-06T12:24:39Z")

</div>

Thanks for the reference!! After reading carefully the notes and yes, I’m not treating the time derivative term correctly.

However, I’m a bit confused when applying it to my case. As I said, I’m dealing with a PDE, so my solution `V`, depends on time and space, so I’m doing two discretizations here, in space I’m using a Galerkin method (to use the Fourier Transforms) and in time I’m approximating the time derivative via an explicit Euler scheme, so I can compute `V_i+1`:

```julia
# Explicit Euler scheme (in Fourier domain):
# α(v_{i+1}-v_i)/Δt = v_i+(î_i-v_i+a_i) <=> v_i+1 = v_i+(Δt/α)*(î_i-v_i+a_i)
@. v = v + (dt/α)*(î - v + a)

```

where `v`, `î` and `a` are the Fourier coefficients of the respective functions `V`, `I` and `A`.

But, when looking to Algorithm 1 of the notes, we no longer have the idea of a finite difference, right? Because the Fourier coefficients of the derivative, `Y_k'`, are obtained in step 2:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/6/46d770f5d9fb9b0ded62b045d754c9698c9ac815.png)

Since I’m using the finite difference `(v_{i+1}-v_i)/Δt` (fourier coefficients) this actually corresponds to `Y_k'` of the step 3 of Algorithm 1, or I’m viewing this the wrong way? 😕 Even if I deal with the derivative as a “whole” as proposed in Alg. 1 how it’s supposed to compute `v_i+1`?

Sorry for the newbie questions.

Really appreciate the help!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 6, 2021, 2:02pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/4 "2021-03-06T14:02:16Z")

</div>

> [@tiago](#):
>
> Since I’m using the finite difference `(v_{i+1}-v_i)/Δt` (fourier coefficients) this actually corresponds to `Y_k'` of the step 3 of Algorithm 1

The notes are only about differentiating in the dimensions you are Fourier transforming. If you are Fourier transforming in space but not time, then the notes have nothing to do with your time derivatives.

(You might consider using DifferentialEquations.jl for your timestepping, i.e. using the method of lines.)

---

<div class="post-metadata">

**Author:** ![jgidi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgidi/32/26206_2.png) [@jgidi](https://discourse.julialang.org/u/jgidi)\
**Post date:** [March 6, 2021, 6:30pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/5 "2021-03-06T18:30:53Z")

</div>

This is just a guess of what the problem might be, but I assume that your algorithm at some point makes use of the conjugate variable in the Fourier space (usually called frequency or wavenumber). If so, are you taking into account the FFTW convention that the first component of the transformed arrays correspond to the origin (zero) of the frequency/wavenumber?

What I mean is that, assuming a space discretization step `dx` and an even array size `N`, the wavenumber has to be defined as

```julia
k = vcat( 0:div(N, 2)-1, -div(N, 2):-1 ) * 2pi/(N*dx)

```

or, equivalently

```julia
k = fftshift( -div(N, 2):(div(N, 2)-1) ) * 2pi/(N*dx)

```

but not simply as a monotonic increasing variable.

However, [FFTW.jl](https://github.com/JuliaMath/FFTW.jl) provides the function `fftfreq` to make this easier,

```julia
k = fftfreq(N, 2pi/dx)

```

which works for even as well as odd `N` (You can see the docs to know the difference).

Let us know if it works!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 6, 2021, 9:10pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/6 "2021-03-06T21:10:10Z")

</div>

> [@jgidi](#):
>
> the FFTW convention that the first component of the transformed arrays correspond to the origin (zero) of the frequency/wavenumber?

(Note that this is not just an FFTW convention. It is the standard definition of the [discrete Fourier transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform).)

---

<div class="post-metadata">

**Author:** ![jgidi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgidi/32/26206_2.png) [@jgidi](https://discourse.julialang.org/u/jgidi)\
**Post date:** [March 6, 2021, 9:17pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/7 "2021-03-06T21:17:53Z")

</div>

> [@stevengj](#):
>
> (Note that this is not just an FFTW convention. It is the standard definition of the [discrete Fourier transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform).)

That is a good point. Thanks for the clarification!

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 8, 2021, 12:17pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/8 "2021-03-08T12:17:25Z")

</div>

@stevengj, Oh I get it, of course, the differentiating is with respect to the space (where I’m doing the Fourier transform) and not with time, sorry my bad 😕. About the timestepping I’m following a specific paper in order to implement this numerical scheme, so I’m a bit stitched to what the authors wrote.

> [@jgidi](#):
>
> at some point makes use of the conjugate variable in the Fourier space (usually called frequency or wavenumber)

Well… I don’t quite use the k, i.e. k is the variable of the Fourier domain. So, for example, when I apply the `fft` to the initial condition `V0`, I simply obtain `v0`, the variable k is kind of implied when using `fft`, right? Here’s a screenshot taken from [Solving PDEs with the FFT, Part 2 [Python]](https://www.youtube.com/watch?v=mMdIxa5qC9Y)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/4/f4d5c8c94137bba2a5a6deb5e988b167353cb17d.png)  
From what I understand, the variable k is only necessary when we differentiate in the dimension where the Fourier transform takes place, that is not my case.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/9/b99684b36df30735d3c905f6ca30584867082a89.png)

To be more specific, I’m dealing with a Neural Field Equation, here’s a simplified version (in 2D):  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/4/c488c53794b684328c83dc4995512e8c324ff324.png)  
Basically, the algorithm that I’m implementing is quite straightforward:

1. Apply `fft` to `V0` to obtain `v0`;
2. Apply `fft` to `I` to obtain `î` in each time step, like this:

```julia
I = [prob.I(i,ti) for i in x] # variable x is 'natural' space
î = fft(I)

```

1. According to the [authors](https://hal.inria.fr/hal-00872132/document) the integral term can be computed in the frequency domain by simply taking into consideration the `fft` of `K` and the non-linear functional `SV`, in each time step;
2. Compute `v_{i+1}` (using a forward Euler scheme) then apply `ifft` to transform the solution back to the natural domain.

I don’t use the frequency variable k, that can be given by the command `fftfreq`. Is this my problem? If so, where enters k?  
(one small remark, I’m using the `rfft`, since my data is real)

Thanks a lot for this help, now I’m trying to know more about the `fftfreq` and trying to understand if it fits for what I’m trying to do.

Cheers!

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 8, 2021, 12:29pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/9 "2021-03-08T12:29:34Z")

</div>

The wavenumber would appear if you had spatial derivatives and used FFT along that variable, which is not the case.

It seems that 1D FFT along time is all that is required?  
Unless for the spatial integral the algorithm you are following uses the FFT in space?

---

<div class="post-metadata">

**Author:** ![jgidi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgidi/32/26206_2.png) [@jgidi](https://discourse.julialang.org/u/jgidi)\
**Post date:** [March 8, 2021, 2:26pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/10 "2021-03-08T14:26:46Z")

</div>

> [@tiago](#):
>
> I don’t use the frequency variable k, that can be given by the command `fftfreq`.

You are right. I can see now that you are using the Fourier Transform to calculate the spatial convolution given by the integral term (a convolution in real space turns into multiplication in the Fourier space). You don’t use `k` and that’s correct.

> [@tiago](#):
>
> (one small remark, I’m using the `rfft` , since my data is real)

Now I find specially weird that using `fftshift` helps you at all, since the real transforms only compute one half of the spectrum, and `ffshift` will swap the first and second half of the arrays wrongfully assuming that they correspond to different halves of the spectrum.

However, you said your problem is multidimensional. I gave a quick look at the [paper you linked](https://hal.inria.fr/hal-00872132/document), and from what I see, you should only transform the space dimensions (not the time!) and you should not use any `fftshift`. Can we know exactly which dimensions are you transforming, and which dimensions are you `ffftshift`-ing?

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 8, 2021, 3:12pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/11 "2021-03-08T15:12:23Z")

</div>

> [@rafael.guerra](#):
>
> It seems that 1D FFT along time is all that is required?  
> Unless for the spatial integral the algorithm you are following uses the FFT in space?

I don’t use the `fft` in time, only in space, which can be one or two-dimensional space (and this behavior happens in 2 dimensions as well, if I put the `(i)fftshift`’s every time that I perform the `fft` and its inverse).

> [@jgidi](#):
>
> Now I find specially weird that using `fftshift` helps you at all, since the real transforms only compute one half of the spectrum, and `ffshift` will swap the first and second half of the arrays wrongfully assuming that they correspond to different halves of the spectrum.

I get what you’re saying, but the same behavior happens if I use the `rfft` or `fft`, with or without the use of `(i)fftshift`.

> [@jgidi](#):
>
> However, you said your problem is multidimensional … and from what I see, you should only transform the space dimensions (not the time!)

Yes, I’m only doing the transforms along the spacial variables, not in the time variable. If the domain is 1D, I perform `fft` to vectors and if it’s 2D I perform `fft` to matrices, so the `(i)fftshift` are applied to vectors or matrices, according to the dimension of the space.  
Also, the authors of that paper, made available a script in python of that algorithm where they use the `(i)fftshift` when computing the transforms.

> [@jgidi](#):
>
> I can see now that you are using the Fourier Transform to calculate the spatial convolution given by the integral term (a convolution in real space turns into multiplication in the Fourier space)

I think the problem might be here, the use of the Fourier transform to compute the convolutional term. One thing, the authors treat a slightly different equation that I posted, they work with the delayed form. The whole new idea of that paper is treating the integral term in two convolutions, one in space and other in time,  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/a/0a0ac857fa3ea9d4d965a36db911df49c814332f.png)  
But in practical terms, the convolution in space is the only one computed using `fft`, the other is treated in a different way.

---

<div class="post-metadata">

**Author:** ![jgidi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgidi/32/26206_2.png) [@jgidi](https://discourse.julialang.org/u/jgidi)\
**Post date:** [March 8, 2021, 4:47pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/12 "2021-03-08T16:47:24Z")

</div>

> [@tiago](#):
>
> in practical terms, the convolution in space is the only one computed using `fft`, the other is treated in a different way.

I understand that they solve the integral term as shown in Eq. (10) from the paper you posted

 ![Screenshot_20210308_132219](https://global.discourse-cdn.com/julialang/original/3X/0/3/03260db35ecee9dda827c4d164d94b7706ebe52d.png)

where the kernel L is known beforehands, and its Fourier transform, DFT[L(T\_u)] is precomputed for all u.

Looking for the authors, I found two equivalent python scripts for an [older paper (2010)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.82.055701) in which they solve the integral term. In [one of the scripts](https://github.com/rougier/delayed-neural-field/blob/master/DNF.py), they perform the `fft`’s as usual, but they `fftshift` the kernel when they first compute it. However, in the other script (‘dnf\_template.py’ distributed in the [DNF GUI simulator](https://gforge.inria.fr/projects/nfsimulator/)) they use the Fourier Transforms as you do; concatenating `fft`’s and `fftshift`’s.

Is it possible that the shifting becomes necessary due to the way the kernel is defined on each script? You could try to use the `rfft` without shifting in general, and only shift the kernel when it is first computed.

Of course, I can not know with certainty what the problem is, as I just gave a fast look over the material. I am just guessing here, and it’s very likely that I have got something wrong. However, hope it helps!

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 8, 2021, 8:54pm UTC](https://discourse.julialang.org/t/the-need-of-i-fftshifts-to-compute-a-pde-using-fourier-transform/56587/13 "2021-03-08T20:54:50Z")

</div>

> [@jgidi](#):
>
> In [one of the scripts](https://github.com/rougier/delayed-neural-field/blob/master/DNF.py), they perform the `fft` 's as usual, but they `fftshift` the kernel when they first compute it

I knew the older paper from 2010, but I was unaware of this script. And in fact, if I just perform the `fftshift` to the kernel `k`, when it’s pre-computed for all `u` and the correct result is obtained! So, I no longer need to perform the chain of `(i)fftshift` everytime I compute the transform.

**However** , in 2D the output is still wrong and this is actually really strange, since the algorithm is exactly the same as the 1D algorithm. The only difference is that one deals with matrices and the other with vectors.  
Since I’m using the multidimensional `rfft`, this means dealing with matrices of dimension `div(N,2)+1 X N`, probably it’s not a good idea to apply the `fftshift` to non-square matrices. This problem doesn’t exist in 1D because the vectors are just `div(N,2)` long. I will try to use the ‘normal’ `fft` to see what I can get.

I was reading about computing convolutions using `fft`, and, as you said, that corresponds to a multiplication in the Fourier domain. However, computing the product of two `fft` corresponds to a circular convolution, not linear.  
Probably that’s why the need of shifting the kernel first

Thanks for the hint!!