Power Spectral Density from jTransforms DoubleFFT_

问题:

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

回答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]

回答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

回答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]))