30
$\begingroup$

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]

enter image description here

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

enter image description here

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}]

enter image description here

$\endgroup$
  • $\begingroup$ for 1D the operation should just swap the 2 halfs of the vector around. For 2D, it should swap the 4th quadrant with the second, and the 3rd quadrant with the first. Looking at the matrix as {{ first, second},{fourth, third}}. $\endgroup$ – Nasser Oct 7 '13 at 3:31
  • $\begingroup$ The correct answer is here: mathematica.stackexchange.com/questions/33149/… Just find the 0 in the vector ssf which shows the frequencies for all the terms. $\endgroup$ – bill s Oct 7 '13 at 4:15
  • $\begingroup$ I am afraid the answer will depend on the number of points. Is N odd or is it even? Have a look at the domain function I defined here mathematica.stackexchange.com/questions/18082/…. Try it out with a superposition of sinusoids of known frequencies and see if you are getting the peak at exactly the right position. If it does not, let me know! $\endgroup$ – Peltio Oct 7 '13 at 11:18
  • $\begingroup$ @bills that means none of the the two methods in the question is correct? Do you have reference for the answer you linked? Thanks. $\endgroup$ – xslittlegrass Oct 7 '13 at 15:41
  • $\begingroup$ The answer at the link plots all the values, those corresponding to both negative and positive frequencies. In the above, you have only shown the positive frequencies, so of course they look different. Have you tried the method in the link? If you choose to only plot the positive, I'll bet it looks more like the above. $\endgroup$ – bill s Oct 7 '13 at 17:17
31
$\begingroup$

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]

transform

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.

$\endgroup$
  • $\begingroup$ @xslittlegrass - Very nice and thorough rewrite. Thanks! $\endgroup$ – bill s Oct 26 '13 at 1:52
  • $\begingroup$ What if n is odd ? $\endgroup$ – rhermans Mar 17 '14 at 17:50
  • 1
    $\begingroup$ @rhermans I usually drop one data point if n is odd :) $\endgroup$ – xslittlegrass Apr 7 '14 at 5:12
  • 1
    $\begingroup$ @rhermans see the update for odd n :) $\endgroup$ – xslittlegrass Apr 7 '14 at 17:06
  • $\begingroup$ Why is this ssf = RotateRight[Range[-n/2, n/2 - 1]/(n dt), n/2]; correct? if you want to start from negative to positive? e.g. n can be 30. Then you have Range[-15,14] = {-15,-14,-13,....,14}. Rotating to the right would give you { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1} No?? I am very confused right now. $\endgroup$ – el psy Congroo Jun 20 '17 at 20:18
11
$\begingroup$

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}]

squared magnitude of shifted FFT

(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.stack.imgur.com/GO2jS.png"]

cameraman

Now, transform and shift:

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

Show the (logarithmically transformed) magnitude:

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

shifted squared magnitude of 2D FFT

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

ImagePeriodogram[cam]

ImagePeriodogram result

$\endgroup$
10
$\begingroup$

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}] 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.