1D real FFT help

fruki — Mon, 01 Aug 2011 03:45:32 GMT

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.

1D real FFT help

rodney_thomson — Mon, 01 Aug 2011 06:40:21 GMT

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

topic 1D real FFT help in Intel® Integrated Performance Primitives

1D real FFT help

1D real FFT help