What's the correct way to shift zero frequency to the center of a Fourier Transform?

Question

I'm trying to apply a Fourier transform of a one dimensional list of a time history of some quantity using the Fourier function. I'm interested in the frequency spectrum, but the problem is that the Fourier function uses the fast Fourier transform algorithm which places the zero frequency at the beginning, complicating my analysis of the results.

So how can I shift the zero frequency to the center?

I tried to search for the solution and found two methods, which seem to give contradicting answers.

Method 1 (from course notes available here). It simply rotates the list before and after the Fourier transform:

DFT1[ls_?(EvenQ@Length[#] &), dt_] := Module[{N0, fft},
  N0 = Length[ls];
  fft = RotateRight[
    dt*Fourier[RotateLeft[ls, N0/2 - 1], FourierParameters -> {1, 1}],
     N0/2 - 1];
  fft
  ]

Method 2 (from some code of my professor). It rotates the list only after the transform, and adds a phase shift:

DFT2[ls_?(EvenQ@Length[#] &), dt_] := Module[{N0, dw, wls, fft},
  N0 = Length[ls];
  dw = (2 π)/(N0 dt);
  wls = dw Range[-(N0/2), N0/2 - 1];
  fft = dt*
    Reverse[RotateRight[Fourier[ls, FourierParameters -> {1, -1}], 
      N0/2 - 1]]*Exp[(I π)/dw wls];
  fft
  ]

If we use these two methods on an example, we get different answers:

dt = 0.05;
els = Table[Sin[ t] Sin[t/40]^2, {t, 0., 40 π, dt}];
ListPlot[els, Joined -> True]

Row[ListPlot[{#[DFT1[els, dt]], #[DFT2[els, dt]]}, Joined -> True, 
    PlotRange -> {{1200, 1300}, All}, ImageSize -> 300] & /@ {Abs, Re,
    Im}]

Questions:

1. Which of the methods is correct (if neither is correct then what is the right way)?
2. What are the frequencies corresponding to the Fourier transform results? For example, should the frequencies range from $\{-\frac{N0\Delta \omega}{2},\frac{(N0-1)\Delta \omega}{2}\}$ or $\{-\frac{(N0-1)\Delta \omega}{2},\frac{N0\Delta \omega}{2}\}$? $N0$ is the length of the data, $\Delta t$ is the time interval of the data, $\Delta \omega=\frac{2\pi }{N0 \Delta t}$.
3. What difference does it make if N0 is even or odd?

Update:

This is the result using the method in Bill's answer here to a similar question. The results appear to differ from both of the approaches mentioned above.

n = Length[els];
sampInt = dt;
data = els;
ssf = RotateRight[Range[-n/2, n/2 - 1]/(n sampInt), n/2];
fft = dt Fourier[data, FourierParameters -> {1, 1}];

Row[ListPlot[#@
     Sort[Transpose[{ssf, fft}], #1[[1]] < #2[[1]] &][[All, 2]], 
    PlotRange -> {{1200, 1300}, All}, Joined -> True, 
    ImageSize -> 300] & /@ {Abs, Re, Im}]

Answer 1

The Fourier transform is defined as:

$$ \begin{align*} H(f) &= \int_{-\infty}^\infty h(t) \mathrm e^{2\pi \mathrm i f t}\mathrm dt\\ h(t) &= \int_{-\infty}^\infty H(f) \mathrm e^{-2\pi \mathrm i f t}\mathrm df \end{align*} $$

where $ h(t) $ is the signal, and $ H(f) $ is its Fourier transform; if $ t $ is measured in seconds, then $ f $ is measured in hertz (Hz).

The discrete Fourier transform is defined as:

$$ \begin{align*} H_{f_j} &= \frac{1}{N}\sum_{k}h_{t_k}\mathrm e^{2\pi \mathrm i f_j t_k}\\ h_{t_j} &= \frac{1}{N}\sum_{k}H_{f_k}\mathrm e^{-2\pi \mathrm i f_k t_j}\\ \end{align*} $$

where the $ t_k $ are the time corresponding to my signal in the time domain $ h_{t_k} $, $ f_k $ are the corresponding frequency to my signal in the frequency domain, and $ N $ is the number of points of the signal data.

For example, if my signal data has $ N $ points, and it was taken at a time interval of $ \Delta t $, then:

• The period of this signal is $$ N*\Delta t $$

• The frequency interval of the Fourier transform is $$ \Delta f=\frac{1}{N*\Delta t} $$

• The time range of the signal is $$ \text{from}~~~0~~~\text{to}~~~(N-1)*\Delta t $$

• The frequency range of the Fourier transform is $$ \text{from}~~~0~~~\text{to}~~~(N-1)*\Delta f $$

We can also see from the definition of the discrete Fourier transform that for any frequency, we can shift it with $ \frac{1}{\Delta t} $ and get the same answer since $$ \mathrm e^{2\pi \mathrm i (f_k+1/\Delta t)t_k}=\mathrm e^{2\pi \mathrm i f_k t_k} $$ which means that frequency $ (N-1)\Delta f $ is the same as $ -\Delta f $. So we can change all the second half of the higher frequencies into the their negative counterpart, so that the frequencies are now symmetric around zero frequency:

$$ \begin{array}{cccccccccc} 0 & \Delta f & 2\Delta f & ... & (\frac{N}{2}-1)\Delta f & \frac{N}{2}\Delta f & (\frac{N}{2}+1)\Delta f & ... & (N-1)\Delta f\\ \downarrow & \downarrow &\downarrow && \downarrow & \downarrow & \downarrow & & \downarrow\\ 0 & \Delta f & 2\Delta f & ... & (\frac{N}{2}-1)\Delta f &-\frac{N}{2}\Delta f &(-\frac{N}{2}+1)\Delta f & ... &-\Delta f& \end{array} $$

If we want to arrange the frequencies from negative to positive, we can simply rotate $ N/2 $ places to the right.

Here's one way to go about it, building on this answer and using the OP's function.

dt = 0.05;
els = Table[Sin[t] Sin[t/40]^2, {t, 0., 40 π, dt}];
n = Length[els];
ssf = RotateRight[Range[-n/2, n/2 - 1]/(n dt), n/2];
fft = Fourier[els, FourierParameters -> {-1, 1}];
ListPlot[Select[Transpose[{ssf, Abs[fft]}], Abs[#[[1]]] < 0.5 &], 
     PlotRange -> All, Filling -> Axis]

The only difference here is selecting the $x$-axis values less than 0.5 Hz (needed because otherwise, the data is hard to see). You can eyeball this data to see that it's about right: the curve els has about 20 periods in 2500 samples. This is a normalized frequency of about 0.016, as seen in the graph. Looking at the ssf vector, the exact frequency is 0.159109, which occurs in location 21 of the ssf vector. Location 21 of the fft vector is the one with the magnitude at about 0.25, the peak value.

Update:

For $n$ an odd number, the frequency arrangement is slightly different:

$$ \begin{array}{ccccccccc} 0 & \Delta f & 2\Delta f & ... & (\frac{N-1}{2})\Delta f & (\frac{N+1}{2})\Delta f & ... & (N-1)\Delta f\\ \downarrow & \downarrow &\downarrow && \downarrow & \downarrow & & \downarrow\\ 0 & \Delta f & 2\Delta f & ... & (\frac{N-1}{2})\Delta f &(-\frac{N-1}{2})\Delta f & ... &-\Delta f& \end{array} $$

so for the data with an odd number of points, we only need to change the frequencies to be

ssf = RotateRight[Range[-(n-1)/2, (n-1)/2 ]/(n dt), (n+1)/2];

and we will be OK.

Answer 2

I am writing this answer to put on record how to do the equivalent of fftshift() in MATLAB, for both the one- and two-dimensional cases.

In Mathematica, the shifting function can be implemented in terms of RotateRight[], like so:

fftshift[dat_?ArrayQ, k : (_Integer?Positive | All) : All] :=
         Module[{dims = Dimensions[dat]}, 
                RotateRight[dat, If[k === All, Quotient[dims, 2], 
                                    Quotient[dims[[k]], 2] UnitVector[Length[dims], k]]]]

The function ifftshift() can be implemented similarly (spot the difference!):

ifftshift[dat_?ArrayQ, k : (_Integer?Positive | All) : All] := 
          Module[{dims = Dimensions[dat]}, 
                 RotateRight[dat, If[k === All, Ceiling[dims/2], 
                                     Ceiling[dims[[k]]/2] UnitVector[Length[dims], k]]]]

---

Using a slightly modified version of the OP's example, here's how to shift in the one-dimensional case:

h = 1/20;
els = Table[N[Sin[2 π t] Sin[π t/20]^2], {t, 0, 20 - h, h}];

ListPlot[Abs[fftshift[Fourier[els, FourierParameters -> {1, -1}]]]^2/Length[els],
         DataRange -> {-1/(2 h), 1/(2 h)}, Filling -> Axis, PlotRange -> {{-2, 2}, All}]

(Note that I am using the square of the magnitude.)

---

For the two-dimensional case, let's use a classical test image:

cam = Import["https://i.sstatic.net/GO2jS.png"]

Now, transform and shift:

tst = Abs[fftshift[Fourier[ImageData[cam], FourierParameters -> {1, -1}]]]^2;

Show the (logarithmically transformed) magnitude:

Image[Log[1 + tst]] // ImageAdjust

This is effectively what's being done by ImagePeriodogram[]:

ImagePeriodogram[cam]

Answer 3

Thanks for a complete answer bill and Ziofil.

Adding further to your answer, in the case of a 2D Fourier transform of an image:

• Case1: ImagePeriodogram[image] would give you the Fourier transform of the input image right answer, with DC centered in the middle of the resultant image.
• Case2: scaledPower = Image[PeriodogramArray[image]] would give you the Fourier transform of an image with DC peak not at the center of the image but, at the corners.

For Case2: To bring the DC peak frequency back to the center of the image, simply follow the following commands (assuming your image has even width and height):

width = ImageDimensions[scaledPower][[1]];
height = ImageDimensions[scaledPower][[2]];
widthHalf = width * 0.5;
heightHalf = height * 0.5;
firstQuad =  ImageCrop[scaledPower, {widthHalf, heightHalf}, {Left, Bottom}];
secondQuad = ImageCrop[scaledPower, {widthHalf, heightHalf}, {Right, Bottom}];
thirdQuad =  ImageCrop[scaledPower, {widthHalf, heightHalf}, {Right, Top}];
fourthQuad = ImageCrop[scaledPower, {widthHalf, heightHalf}, {Left, Top}];

ImageCollage[{fourthQuad, thirdQuad, firstQuad, secondQuad}]