Power Spectral Density from jTransforms DoubleFFT_1D

Question

I'm using Jtransforms java library to perform analysis on a given dataset.

An example of the data is as follows:

980,988,1160,1080,928,1068,1156,1152,1176,1264

I'm using the DoubleFFT_1D function in jTransforms. The data output is as follows:

10952, -152, 80.052, 379.936, -307.691, 12.734, -224.052, 427.607, -48.308, 81.472

I'm having trouble interpreting the output. I understand that the first element in the output array is the total of the 10 inputs (10952). It's

the other elements of the output array that i don't understand. Ultimately, I want to plot the Power Spectral Density of the input data on a graph and find amounts between 0 and .5 Hz.

The documentation for the jTransform functions states (where a is the data set):

public void realForward(double[] a) computes 1D forward DFT of real data leaving the result in a . The physical layout of the output data is as follows:

if n is even then

a[2*k] = Re[k], 0 <= k < n / 2
a[2*k+1] = Im[k], 0 < k < n / 2
a[1] = Re[n/2]

if n is odd then

a[2*k] = Re[k], 0 <= k < (n+1)/2
a[2*k+1] = Im[k], 0 < k< (n-1)/2
a[1] = Im[(n-1)/2]

This method computes only half of the elements of the real transform. The other half satisfies the symmetry condition. If you want the full real forward transform, use realForwardFull. To get back the original data, use realInverse on the output of this method.

Parameters: a - data to transform

Now using the methods above: (since the length of my data array is 10, the "n is even" methods are used)

Re[0] = 10952
Re[1] = 80.052
Re[2] = -307.691
Re[3] = -224.052
Re[4] = -48.308
Re[5] = 12.734

Im[0] = -152
Im[1] = 379.936
Im[2] = 12.734
Im[3] = 427.607
Im[4] = 81.472

So some questions: Does this output look correct? It seems to me that Re[0] should not be 10952 which is the sum of all elements in the original array.

Seems like the output should be slightly corrected: (am I wrong?)

Re[0] = 80.052
Re[1] = -307.691
Re[2] = -224.052
Re[3] = -48.308
Re[4] = -152

Im[0] = 379.936
Im[1] = 12.734
Im[2] = 427.607
Im[3] = 81.472

Now using the following method posted in the forum:

To get the magnitude of bin k you need to calculate sqrt(re * re + im * im), where re, im are the real and imaginary components in the FFT output for bin k.

For your particular FFT re[k] = a[2*k] and im[k] = a[2*k+1]. Therefore to calculate the power spectrum:

for k in 0 to N/2 - 1
{
    spectrum[k] = sqrt(sqr(a[2*k]) + sqr(a[2*k+1]))
}

Thus:

spectrum[0] = 388.278
spectrum[1] = 307.955
spectrum[2] = 482.75
spectrum[3] = 94.717

Some questions. Does this data look correct? Am I on the right track? Would this spectrum data then plot out something like this:

388.278 at .125 Hz
307.955 at .25 Hz
482.75 at .375 Hz
94.717 at .5 Hz

Am I way off? My goal is to produce a Power Spectral Density bar chart from 0 to .5Hz

Messinger

Isn't this a duplicate of your question from yesterday: #4997247 ? – Pyridine

yes. with some follow up. After I logged in, i was unable to manipulate the old thread. None of my old questions (including the one referenced) showed up in my account. So I moved it and added some followup data. Sorry. I don't want to add confusion. Just had trouble with stackoverflow going from guest to logged in user. I really appreciate the help. I feel like I'm making progress with your help. – Messinger

ask the mods to fix account problems for you. Also how about marking Paul R's answer as accepted and upvoting it, looks like it was pretty useful to you. – Orvil

I would love to upvote the answer but i was unable to on the original thread. That's one reason I moved it. I'm not sure who to ask to fix it but now that I'm logged in I figure it shouldn't be a problem again. Paul R's answer did help. Would love to give him credti but I couldn't figure out how on the original thread. Again I appologize. – Messinger

Merged accounts and merged the questions. All should be well now. (In the future, you can flag for moderator attention to get this sort of thing fixed up. Someone else flagged it this time.) – Penza

Answer 1

I think you need to interpret the output data as follows:

10952       Re[0] = sum of all inputs = DC component
 -152       Re[5] - see note about a[1] being special - there is no Im[0]
   80.052   Re[1]
  379.936   Im[1]
 -307.691   Re[2]
   12.734   Im[2]
 -224.052   Re[3]
  427.607   Im[3]
  -48.308   Re[4]
   81.472   Im[4]

The magnitudes are therefore:

spectrum[0] = 10952
spectrum[1] = sqrt(80.052^2 + 379.936^2) = 388.278
spectrum[2] = sqrt(-307.691^2 + 12.734^2) = 307.427
spectrum[3] = sqrt(-224.052^2 + 427.607^2) = 482.749
spectrum[4] = sqrt(-48.308^2 + 81.472^2) = 94.717

[Sorry about there being two separate answers from me now - I think two related questions got merged while I was working on the new answer]

Pyridine

Thanks Paul. That clears it up greatly! Forgive my lack, but how do the spectrum results[0-4] relate to frequencies between 0 and .5 Hz? – Messinger

for the example above, I think Fs = 1Hz (that way the math works out). – Messinger

@Damon: sorry, yes, I meant Fs = 1 Hz. I'll repost a corrected comment. – Pyridine

@Damon: the FFT doesn't know or care about the sampling rate. You don't say what your sample rate is in the question - let's call it Fs. The centre frequency for bin k is k * Fs / N, where N is the number of points. So e.g. if Fs = 1.0 Hz and N = 10, then bin 0 is DC (0 Hz), bin 1 is 0.1 Hz, bin 2 is 0.2 Hz, bin 3 is 0.3 Hz... – Pyridine

Answer 2

Each entry in the transform represent the (complex) magnitude of the frequency in the sample.

the power density in a given frequency is just the magnitude of the complex amplitude of the transform in that frequency. the magnitude of a complex number is computed from its components and you should not have a problem obtaining this

Each column represents amplitudes for increasing frequencies, starting from 0 (the first entry), then 2 Pi/T (where T is the length of your sample), until the last sample 2*Pi*N /T (where N is the number of samples)

there are other conventions where the transform is returned for the -Pi * N /T frequency up to Pi * N / T, and the 0 frequency component is in the middle of the array

hope this helps

Subshrub

Answer 3

To get the magnitude of bin k you need to calculate sqrt(re * re + im * im), wheer re, im are the real and imaginary components in the FFT output for bin k.

For your particular FFT re[k] = a[2*k] and im[k] = a[2*k+1]. Therefore to calculate the power spectrum:

for k in 0 to N/2 - 1
  spectrum[k] = sqrt(sqr(a[2*k]) + sqr(a[2*k+1]))

Pyridine

Thanks for the help. This is the best info I've come by so far. I've used your spectrum equation to calculate values for the following data set: <br> 1068,1116,1140,1144,1160,1136,1100,1080,1200,1212,1092,952,1060,1068 – Messinger

@Damon, if this helped you , please click right to the answer to make the checkmark become green – Subshrub