Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to Fourier transform 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 and makes the zero frequency at the beginning of the results, which makes me difficult to analysis.

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

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

Method 1(from a course notes at here). It simply rotate 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 rotate only after the transform, and attaches some phase:

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 one is correct (if none is correct then what is the right way)?
  2. What is 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 is the difference between N0 is even or odd?

Update:

This is the result using the method in bill's answer here to a similar question. The results looks different than any of the two 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

share|improve this question
    
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}}. –  Nasser Oct 7 '13 at 3:31
    
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. –  bill s Oct 7 '13 at 4:15
    
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! –  Peltio Oct 7 '13 at 11:18
    
@bills that means none of the the two methods in the question is correct? Do you have reference for the answer you linked? Thanks. –  xslittlegrass Oct 7 '13 at 15:41
    
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. –  bill s Oct 7 '13 at 17:17

1 Answer 1

up vote 14 down vote accepted

The Fourier transform is defined as:

$$ H(f)=\int_\infty^\infty h(t) e^{2\pi i f t}dt\\ h(t)=\int_\infty^\infty H(f) e^{-2\pi i f t}df $$

where $h(t)$ is the signal, and $H(f)$ is it's Fourier transform, if $t$ is meassured in second, then $f$ is measured in Hz.

The discrete Fourier transform is defined as:

$$ H_{f_j}=\frac{1}{N}\sum_{k}h_{t_k}e^{2\pi i f_j t_k}\\ h_{t_j}=\frac{1}{N}\sum_{k}H_{f_k}e^{-2\pi i f_k t_j}\\ $$

where $t_k$ are the time corresponding to my signal in time domain $h_{t_k}$, $f_k$ are the corresponding frequency to my signal in 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 $$ e^{2\pi i (f_k+1/\Delta t)t_k}=e^{2\pi 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 frequency from negative to positive, we can simplify rotate N/2 to the right.

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

dt = 0.05;
els = Table[Sin[t] Sin[t/40]^2, {t, 0., 40 \[Pi], 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]

enter image description here

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

share|improve this answer
    
@xslittlegrass - Very nice and thorough rewrite. Thanks! –  bill s Oct 26 '13 at 1:52
    
What if n is odd ? –  rhermans Mar 17 at 17:50
1  
@rhermans I usually drop one data point if n is odd :) –  xslittlegrass Apr 7 at 5:12
1  
@rhermans see the update for odd n :) –  xslittlegrass Apr 7 at 17:06

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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