98

I need some help understanding the output of the DFT/FFT computation.

I'm an experienced software engineer and need to interpret some smartphone accelerometer readings, such as finding the principal frequencies. Unfortunately, I slept through most of my college EE classes fifteen years ago, but I've been reading up on DFT and FFT for the last several days (to little avail, apparently).

Please, no responses of "go take an EE class". I'm actually planning to do that if my employer will pay me. :)

So here is my problem:

I've captured a signal at 32 Hz. Here is a 1 second sample of 32 points, which I've charted in Excel.

enter image description here

I then got some FFT code written in Java from Columbia University (after following the suggestions in a post on "Reliable and fast FFT in Java").

The output of this program is as follows. I believe it is running an in-place FFT, so it re-uses the same buffer for both input and output.

Before: 

Re: [0.887  1.645  2.005  1.069  1.069  0.69  1.046  1.847  0.808  0.617  0.792  1.384  1.782  0.925  0.751  0.858  0.915  1.006  0.985  0.97  1.075  1.183  1.408  1.575  1.556  1.282  1.06  1.061  1.283  1.701  1.101  0.702  ]

Im: [0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ]

After: 

Re: [37.054  1.774  -1.075  1.451  -0.653  -0.253  -1.686  -3.602  0.226  0.374  -0.194  -0.312  -1.432  0.429  0.709  -0.085  0.0090  -0.085  0.709  0.429  -1.432  -0.312  -0.194  0.374  0.226  -3.602  -1.686  -0.253  -0.653  1.451  -1.075  1.774  ]

Im: [0.0  1.474  -0.238  -2.026  -0.22  -0.24  -5.009  -1.398  0.416  -1.251  -0.708  -0.713  0.851  1.882  0.379  0.021  0.0  -0.021  -0.379  -1.882  -0.851  0.713  0.708  1.251  -0.416  1.398  5.009  0.24  0.22  2.026  0.238  -1.474  ]

So, at this point, I can't make heads or tails of the output. I understand the DFT concepts, such as the real portion being the amplitudes of the component cosine waves and the imaginary portion being the amplitudes of the component sine waves. I can also follow this diagram from the great book "The Scientist and Engineer's Guide to Digital Signal Processing": enter image description here

So my specific questions are:

  1. From the output of the FFT, how do I find the "most occurring frequencies"? This is part of my analysis of my accelerometer data. Should I read the real (cosine) or imaginary (sine) arrays?

  2. I have a 32-point input in the time domain. Shouldn't the output of the FFT be a 16-element array for reals and a 16-element array for imaginary? Why does the program give me real and imaginary array outputs both of size 32?

  3. Related to the previous question, how do I parse the indexes in the output arrays? Given my input of 32 samples sampled at 32 Hz, my understanding is that a 16-element array output should have its index uniformly spread up to 1/2 the sampling rate (of 32 Hz), so am I correct in understanding that each element of the array represents (32 Hz * 1/2) / 16 = 1 Hz?

  4. Why does the FFT output have negative values? I thought the values represent amplitudes of a sinusoid. For example, the output of Real[ 3 ] = -1.075 should mean an amplitude of -1.075 for a cosine wave of frequency 3. Is that right? How can an amplitude be negative?

Cœur
  • 37,241
  • 25
  • 195
  • 267
stackoverflowuser2010
  • 38,621
  • 48
  • 169
  • 217
  • What would you like to compute from the accelerometer readings: velocity, distance? The noise of the accelerometer readings follows Gaussian distribution and I cannot see how fitting a sine wave would remedy that. – Ali Jul 19 '11 at 11:47
  • 3
    the java tag should be removed as it's more generic than to a specific language – user3791372 Mar 24 '16 at 23:43
  • Looking at the souce of Columbia University, it's not efficient at all. It's a plain, unoptimized implementation of Cooley-Tucky with butterfly lookup tables, and the bit-reversal is done manually instead of using existing library functions – Mark Jeronimus May 15 '20 at 14:54
  • @MarkJeronimus: Can you recommend an efficient FFT implementation in Java? If I recall correctly, the reason I went with the Columbia University code was the the FFTW library was too complex to run on an Android smartphone. – stackoverflowuser2010 May 15 '20 at 17:45
  • I found some scattered 'optimized' implementations, but they are basically one algorithm *per* N size, so if you need a range of sizes you need all those routines. In practice I've mainly used Intel Integrated Performance Primitives (yes, from Java, through JNA), but that's non-free. At home I use basically the same algorithm you linked, but written from scratch in 2005 using a textbook. Its just FFT (Fast Fourier Transform), nothing so 'Fast' about it to justify the name 'Fast FFT'. – Mark Jeronimus May 17 '20 at 16:58

4 Answers4

97

  1. You should neither look for the real or imaginative part of a complex number (that what's your real and imaginary array is). Instead you want to look for the magnitude of the frequency which is defined as sqrt (real * real + imag * imag). This number will always be positive. Now all you have to search is for the maximum value (ignore the first entry in your array. That is your DC offset and carries no frequency dependent information).

  2. You get 32 real and 32 imaginary outputs because you are using a complex to complex FFT. Remember that you've converted your 32 samples into 64 values (or 32 complex values) by extending it with zero imaginary parts. This results in a symetric FFT output where the frequency result occurs twice. Once ready to use in the outputs 0 to N/2, and once mirrored in the outputs N/2 to N. In your case it's easiest to simply ignore the outputs N/2 to N. You don't need them, they are just an artifact on how you calculate your FFT.

  3. The frequency to fft-bin equation is (bin_id * freq/2) / (N/2) where freq is your sample-frequency (aka 32 Hz, and N is the size of your FFT). In your case this simplifies to 1 Hz per bin. The bins N/2 to N represent negative frequencies (strange concept, I know). For your case they don't contain any significant information because they are just a mirror of the first N/2 frequencies.

  4. Your real and imaginary parts of each bin form a complex number. It is okay if real and imaginary parts are negative while the magnitude of the frequency itself is positive (see my answer to question 1). I suggest that you read up on complex numbers. Explaining how they work (and why they are useful) exceeds what is possible to explain in a single stackoverflow-question.

Note: You may also want to read up what autocorrelation is, and how it is used to find the fundamental frequency of a signal. I have a feeling that this is what you really want.

André Chalella
  • 13,788
  • 10
  • 54
  • 62
Nils Pipenbrinck
  • 83,631
  • 31
  • 151
  • 221
  • 2
    Thanks. Regarding 1: I saw this Matlab page that shows a frequency spectrum (http://www.mathworks.com/help/techdoc/ref/fft.html). On that page is a plot with the title of "Single-Sided Amplitude Spectrum of y(t)". Is that plotting the magnitude of the frequency as you suggested, sqrt(real^2 + img^2)? Regarding 3: I still don't get the 2Hz per bin result. In my case, N=32 and freq=32, right? So there are N/2=32/2=16 bins, and the highest frequency (Nyquist) is freq/2=32/2=16 Hz, resulting in 16 Hz per 16 bins, giving 1 Hz per bin? – stackoverflowuser2010 Jul 19 '11 at 04:56
  • 1
    Yes, the plot shows the magnitude of the spectrum - |Y(f)|. The absolute-value bars means magnitude. Bin width = sample rate / FFT size. Your sample rate is 32 hz, your FFT size is 32. Yes, you are right about the bin width! – Matt Montag Jul 20 '11 at 07:51
13

You already have some good answers, but I'll just add that you really need to apply a window function to your time domain data prior to the FFT, otherwise you will get nasty artefacts in your spectrum, due to spectral leakage.

Paul R
  • 208,748
  • 37
  • 389
  • 560
6

1) Look for the indices in the real array with the highest values, besides the first one (that's the DC component). You'll probably need a sample rate considerably higher than 32 Hz, and a larger window size, to get much in the way of meaningful results.

2) The second half of both arrays is the mirror of the first half. For instance, note that the last element of the real array (1.774) is the same as the second element (1.774), and the last element of the imaginary array (1.474) is the negative of the second element.

3) The maximum frequency you can pick up at a sample rate of 32 Hz is 16 Hz (Nyquist limit), so each step is 2 Hz. As noted earlier, remember that the first element is 0 Hz (i.e, the DC offset).

4) Sure, a negative amplitude makes perfect sense. It just means that the signal is "flipped" -- a standard FFT is based on a cosine, which normally has value = 1 at t = 0, so a signal which had value = -1 at time = 0 would have a negative amplitude.

  • Thanks for the reply. (1) Do you mean that I can ignore the imaginary (sine) array, and if so, why? Surely the sine component must be important? (2) Why does this mirroring occur? Is it just the result of the FFT algorithm? Do most people ignore the mirrored half? (3) How did you calculate the steps of 2Hz? I understand the Nyquist limit of 16Hz, so if there are 16 (non-mirrored) array elements, each element must be 16 Hz/16 = 1 Hz each? (4) To find principal frequencies, do I just just take the absolute value of the amplitude values in the output arrays? – stackoverflowuser2010 Jul 19 '11 at 01:10
  • 1
    You shouldn't look in the real array for the highest value, and you can't ignore the sine/I array. Instead you want the magnitude of the composite complex vector. Mirroring occurs because half the input (the I array) is all zeros, so the result has half the degrees of freedom. You can ignore it if your data is strictly real. – hotpaw2 Jul 19 '11 at 22:05
  • (3), the value of step=2Hz remains implicit to me so far. We have 16 bins, represented by array of lenght=16. We need to describe all frequencies from 0Hz to 16Hz. I assume every bin describes a piece of that range, doesn't it? – krafter Jun 16 '13 at 17:56
  • @krafter I think it's halved because you can't deduce a frequency from a single value (as there is no repetition). – JVE999 Jun 06 '20 at 23:56
6

Note that the "most occurring frequency" can get splattered into multiple FFT bins, even with a window function. So you may have to use a longer window, multiple windows, or interpolation to better estimate the frequency of any spectral peaks.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153