Winterwurst - 9 months ago 37

Java Question

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);`

`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:

Source (Stackoverflow)