0

Connected with this topic , I would like to have some clarifications about the Fourier Transform and the Python tools that make it. Having my data set whose plot is showed in the above link, my array of data is stored in rate. Then, if I check the result from the Fourier transform:

yfft = np.fft.rfft(rate)
y_smooth = np.fft.irfft(yfft)
print y_smooth

I get:

[  3.20000028e+01+0.j          -7.92633533e+00-7.75032655j
-5.34245688e+00-3.82390808j  -2.39880935e+00-0.60037314j
-8.71249166e-01-0.04362804j  -8.52507194e-02+1.09973189j
-1.58332047e-01+1.15256295j  -3.38909974e-01+0.96327143j
 5.51323015e-02+0.43936749j   3.81205886e-01-0.02489853j
 2.99110153e-01-0.26319101j   3.28283652e-01+0.11668814j
 1.08849205e-02-0.11317937j   3.79952011e-01+0.22278952j
 7.37720998e-02-0.26810083j  -2.14247586e-01-0.25950865j
-1.32242152e-01+0.j        ]

So, I am quite confused about this result, and I would need some tips how to interpret these numbers.

1) The array contains $17$ complex number : does it mean that my signal can be approximated by a superimposition of $17$ harmonics?

2) Are these numbers the power of each harmonic, or what?

3) are they ordered by the frequency order, or by power? Let's say I want to put to zero the higher orders, which ones should I put to zero?

Community
  • 1
Py-ser
  • 235
  • 1
  • 2
  • 7

2 Answers2

3

The array contains 17 complex number : does it mean that my signal can be approximated by a superimposition of 17 harmonics?

16.

The first component is the DC component.

Note that this is not specific to your signal. Any real sequence of length 32 can be decomposed as 1 real DC component, 16 complex amplitudes, and 1 real amplitude (32 real number ins, 32 real numbers out).

Are these numbers the power of each harmonic, or what?

They are the complex amplitudes.

are they ordered by the frequency order, or by power? Let's say I want to put to zero the higher orders

They are ordered by frequency.

which ones should I put to zero?

The last ones, but why would you want to do that?

pichenettes
  • 19,413
  • 1
  • 50
  • 69
  • Thanks. Because I want to smooth a signal leaving off the higher orders. – Py-ser Apr 16 '14 at 07:26
  • Then you are probably doing it wrong: http://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins . Design a filter with all the properties you want, and apply it to the time-domain data. – pichenettes Apr 16 '14 at 07:47
  • I think in my case (see the linked discussion in the main topic) it is fine. At least, people more experienced than me do it in such cases. If you can, please, take a look at it and tell me if you still think it is not the case to proceed this way. I will not use it anyway (since I will use the Savitzky-Golay), but this post serves me to understand the principles, so it would be great to understand if the next time I don't even have to think about FFT smoothing. – Py-ser Apr 16 '14 at 08:08
3

Your data consist of 32 real values. Its discrete Fourier transform (DFT), which you computed using an FFT (Fast Fourier Transform) algorithm also has 32 (complex) values. However, since your data are real, the Fourier transform is symmetric and the last 15 values are redundant (i.e. can be computed from the other values). This is why you see 17 remaining non-redundant values. These values are no approximation but they are an exact representation of your data, just in a different form.

The DFT is given by

$$X_k=\sum_{n=0}^{N-1}x_ne^{-jnk2\pi/N},\quad k=0,1,\ldots,N-1$$

where $N$ is the number of data points (32 in your case). If $x_n$ is real-valued you can show that

$$X_k=X^*_{N-k},\quad k=1,\ldots, N-1$$

which means that all information about your data is contained in the values $X_k$ for $k=0,1,\ldots,N/2$. Note that $X_0$ and $X_{N/2}$ are real-valued. So we have $N/2+1$ non-redundant complex number, two of which are real-valued. So you have $N/2-1$ real and imaginary parts, and 2 real numbers, which in total are exactly $N$ real numbers. No surprise that they are sufficient (and necessary) to represent your $N$ data points.

In order to see how you can interpret the values $X_k$ you need to look at the inverse transform

$$x_n=\frac{1}{N}\sum_{n=0}^{N-1}X_ke^{jnk2\pi/N}$$

So the values $X_k$ are the complex amplitudes of the complex exponentials representing your signal. They are ordered according to frequency. $X_0$ represents the DC value, and $X_{N/2}$ (assuming that $N$ is even) represents the value at half the sampling rate.

Matt L.
  • 89,963
  • 9
  • 79
  • 179
  • Thanks! Put in this way, it seems that I need to do as many DFT as much are the component of my signal, instead of just one DFT that gives back all the components. Is it correct to say the amplitudes are the square root of the Power of each considered harmonic? – Py-ser Apr 16 '14 at 07:24
  • @Py-ser The amplitudes are complex, so you could say that the squared magnitude of the amplitudes are equivalent to the energy in the $k^{th}$ frequency component. I do not understand the first part of your comment. – Matt L. Apr 16 '14 at 08:24
  • Forget about the first part. Why do you use the term "energy" instead of power? – Py-ser Apr 16 '14 at 08:27
  • @Py-ser $\sum_n|x_n|^2=\frac{1}{N}\sum_k|X_k|^2$ is the total energy of the signal $x_n$ (look up Parseval's theorem). So $|X_k|^2$ can be considered the energy in bin $k$. – Matt L. Apr 16 '14 at 10:21
  • ... I should have said proportional to the energy, because I forgot the factor $1/N$. – Matt L. Apr 16 '14 at 13:25