Winterwurst Winterwurst - 2 months ago 11
Java Question

Understanding the FFT output

I'm currently occupied in a practicum and my boss wants me to have a FFT ready by the end of the week programmed in java.

Now I already got a code for the FFT from Princeton University: http://introcs.cs.princeton.edu/java/97data/FFT.java

I implemented and extended this code in my project with which I'm now being able to first read the binary input of a signal, then process the FFT on this sample values and then provide the magnitude.

Now I come to my problem.
I inputted the following values which I generated with

final double d = Math.sin(i);
and looped 8 times (this is just for testing purpose, next week I'll have to input real values).

0.0
0.8414709848078965
0.9092974268256817
0.1411200080598672
-0.7568024953079282
-0.9589242746631385
-0.27941549819892586
0.6569865987187891


so those values come from a pure sine (I don't know the right English word but with pure sine I mean a sine with exactly one frequency for example 50 Hz).

The output is now

0.553732750242242
2.3946469565385193 - 2.0970118573701813i
-1.386684423934684 + 0.9155598966338983i
-0.8810419659226628 + 0.28041399267903344i
-0.8075738836045867
-0.8810419659226628 - 0.28041399267903366i
-1.386684423934684 - 0.9155598966338983i
2.394646956538519 + 2.0970118573701817i


And the magnitudes of the output

0.553732750242242
3.183047718211326
1.6616689248786416
0.9245901540720989
0.8075738836045867
0.924590154072099
1.6616689248786416
3.183047718211326


Now I actually expected the output values to be 0 at each frequency sample point up until I reach the frequency the pure sine is dominated by, where the output should be >0 (for example at 50 Hz). At least that's what my boss expected when he gave me this task.

Summary:
So this is what I'm struggling with. I've read another thread asking about a similar issue, but there are still unanswered questions for me. What am I supposed to do with the given output data? How do I find the most occurring frequency?

I really could need some help or explanation of where my thinking is wrong.

Thanks for listening...

Answer

Computing a 512-point fourier transform after applying a simple window function:

w(i)= ((float)i-(float)(n-1f)/2f)

it gives peak at i=25 (max magnitude on result array).

Input was also added with more info such as frequency of sine wave generator (50 Hz) and sampling rate (1kHz or 0.001 seconds per sample) and adding 2PI constant:

initialization now looks like this(as a sin(2xPIxFxi) notation):

for (int i = 0; i < n; i++)
{ 
         v1[i] = ((float)i-(float)(n-1f)/2f)*
             (float)Math.Sin(0.001f* 50f*2f*Math.PI*(float)i);
                               ^      ^  ^                 ^
                               |      |  |                 |
                               |      F  |                 |
                           sampling      2xPi constant   sample bin
                           rate(simulation)
}

and result peak looks like:

v2[22] 2145,21852033773  
v2[23] 3283,36245333956  
v2[24] 6368,06249969329  
v2[25] 28160,6579468591    <-- peak
v2[26] 23231,0481898687  
v2[27] 1503,8455705291  
v2[28] 1708,68502071037  

so we have

25

now these steps are in frequency space and input was at 1kHz rate so you can perceive a maximum of 500 Hz harmonic signal(half of sampling rate).

25*500 = 12500

also result range is 0 to N/2 with other half mirrored, and dividing perceivable frequency range(500) to results range(256 for N=512) gives

48.83 Hz

big portion of the error must have been the window function used at the beginning but v2[26] has higher value than v2[24] so actual pick is somewhere closer to v2[26] and smoothed graph of these points should show 50 Hz.

Ignore first element of result array as it is about constant signal level or infinity wavelength or zero frequency.

Here are dft computing codes just to be sure if fft is returning right results:

//a---->b Fourier Transformation brute-force
__kernel void dft(__global float *aRe,
                  __global float *aIm,
                  __global float *bRe,
                  __global float *bIm)
{
    int id=get_global_id(0); // thread id
    int s=get_global_size(0); // total threads = 512
    double cRe=0.0f;
    double cIm=0.0f;
    double fid=(double)id;
    double fmpi2n=(-2.0*M_PI)*fid/(double)s;
    for(int i=0;i<s;i++)
    {
            double fi=(float)i;
            double re=cos(fmpi2n*fi);
            double im=sin(fmpi2n*fi);

            cRe+=aRe[i]*re-aIm[i]*im;
            cIm+=aRe[i]*im+aIm[i]*re;
    }

    bRe[id]=cRe;
    bIm[id]=cIm;
}

and to be sure, testing result against inverse-transformation to check if original input signal is achieved again:

// a--->b inverse Fourier Transformation brute force
__kernel void idft(__global float *aRe,
                   __global float *aIm,
                   __global float *bRe,
                   __global float *bIm)
{
    int id=get_global_id(0); // thread id
    int s=get_global_size(0); // total threads = 512
    double cRe=0.0f;
    double cIm=0.0f;
    for(int i=0;i<s;i++)
    {
        double re=cos(2.0*M_PI*((double)id)*((double)i)/(double)s);
        double im=sin(2.0*M_PI*((double)id)*((double)i)/(double)s);

        cRe+=aRe[i]*re-aIm[i]*im;
        cIm+=aRe[i]*im+aIm[i]*re;
    }
    cRe/=(double)s;
    cIm/=(double)s;
    bRe[id]=cRe;
    bIm[id]=cIm;
}

I know it's bad to run slow code on fast machines but this looks so much simpler to try and is scalable to many cores(2.4ms for a 320core-gpu including array copies). Maybe fft document tab will be commited by at least 5 users and you can check detailed info from there:

http://stackoverflow.com/documentation/fft/commit

Comments