1D real FFT help

[fruki]

Hi everyone,

I've been trying to build a simple FFT program to analyze different frequencies in a simple wave with no luck. It's tough because I'm fairly new to DSP and I haven't been able to find sample code to do it with IPP. Anyway, this is my shot at it. The problem is that the spikes of the FFT show up in wrong frequencies. In the example I constructed a signal in matlab like so:

>> Fs = 1000; % Sampling frequency

>> T = 1/Fs; % Sample time

>> L = 1000; % Length of signal

>> t = (0:L-1)*T; % Time vector

>> x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

Basically a two sine waves added together at 50Hz and 120Hz. I then copy the output and place it in my code.

The problem is I get the spikes at 52Hz and 124Hz. It's pretty frustrating as I need it as accurate as I can.

Maybe I can't interpret the output, I don't know.

My attempt at doing an FFT goes like this:

#define mag_sqrd(re,im) (re*re+im*im)

const int size = 1000;

double vector[size] = { ... }; //1000 points from matlab

//crapy function that tells me what power of 2 to use

int nextpow2(int a)

{

for(int i = 0 ; i < INT_MAX; i++)

if( pow(2.0,i) >= abs(a) )

return i;

//big error

return -1;

}

int main()

{

IppsFFTSpec_R_32f* pFFTSpec;

IppStatus status;

int FFT_LEN;

Ipp32f *fin;

FFT_LEN = pow(2.0,nextpow2(size)); //FFT_LEN is power of 2

fin = ippsMalloc_32f(FFT_LEN);

for(int i = 0; i < size; i++) //copy first part of vector

fin = vector;

for(int i = size; i < FFT_LEN; i++) //zero-padding the rest

fin = 0;

fclose(f);

//Perform IPP Real FFT

//1. Initializing

status = ippsFFTInitAlloc_R_32f(&pFFTSpec, nextpow2(FFT_LEN), IPP_FFT_DIV_FWD_BY_N, ippAlgHintAccurate);

if( status != ippStsNoErr )

return false;

//output vector

Ipp32f *e = ippsMalloc_32f(FFT_LEN);

//2. real FFT

status = ippsFFTFwd_RToPack_32f( fin, e, pFFTSpec, NULL );

if( status != ippStsNoErr )

return false;

//3. Real to complex

Ipp32fc *c = ippsMalloc_32fc(FFT_LEN);

ippsConjPack_32fc((const Ipp32f*) e, c, FFT_LEN);

//4. Obtain frequency spectrum

double re, im;

for(int i=0;i < FFT_LEN/2;i++)

{

re = c.re;

im = c.im;

fprintf(freal, "%f\\n", re);

fprintf(fimg, "%f\\n", im);

fprintf(fft, "%f\\n", (mag_sqrd(re,im)));

fflush(freal);

}

/// free FFT tables

ippsFree(c);

ippsFree(e);

ippsFree(fin);

ippsFFTFree_R_32f( pFFTSpec );

fclose(freal);

fclose(fimg);

system("pause");

return 0;

}

Any help/advice will be greatly appreciated. (Even if you see other errors! I'm willing to learn)

Thank you very much.

[rodney_thomson · Accepted solution]

Hi fruki,

You are very nearly there, but I think you may have just had a slight problem in determining the frequency value out of your complex array.

I assume that when you say your spikes are at '52Hz' and '124Hz' that you mean the spikes were at the 52nd and 124th indexes into your complex array?

This is correct as the frequency resolution of your FFT is equivalent to sample frequency divided by fft size:

freq_res = Fs / FFT_SIZE;

or approximately 0.9766 in this case (1000/1024).

The first value in your array is the DC component, the second value represents the frequency 0.9766 Hz, the third represents 2*0.9766 etc.

So for the 52nd and 124th:

52nd = 51*0.9766 = 49.8047 Hz
124th = 123*0.9766 = 120.1172 Hz

Which is what you were expecting.

Cheers

Rod