FFTW

Hi,

I'm doing a project at the moment to work out the frequency of pulsars. Basically I take the raw data, put it into time bins and run a FFT over that, then look for peaks in the FFT data to get the frequency of the pulsar. Straight forward enough until you implement it.

My problem is with the FFT, in particular the fftw3 libraries. Well that and not knowing too much about fft's!

Anyway, I wrote this little test program just to get the basics right. Create data for a sine wave then get its frequency from the fft.

[code:1]#include <stdio.h>
#include <math.h>
#include <fftw3.h>

int main(){

int N=5000;
double in[N];
double out[6000];
fftw_plan p;
float x=0;
int y=0 ;
int z=0 ;
FILE *fp;

fp = fopen("datasource.data","w");

/*
* Fill in in[N] with a sine wave to simulate a very regular pulsar
*/
while ( x < 100 ){
y=(int)(100*(1+sin(x)));
in[z]=y;
fprintf(fp, "%f\n", in[z]);
x=x+0.02;
z=z++;
}

fclose(fp);

/*
* Perhaps I should be using complex numbers here, but the inputs are real
* and I expect I'm looking for real outputs...
* fftw_plan_r2r_1d and FFTW_DHT seem to be the necessary choice...
*/
p = fftw_plan_r2r_1d(N, in, out, FFTW_DHT, FFTW_ESTIMATE);

fftw_execute(p); /* repeat as needed */


z=0;
while (z !=5500){

printf("%f\n",out[z]);
z++;


}

fftw_destroy_plan(p);
/* fftw_free(in); */

fftw_free(out);

}[/code:1]

The sine wave that we generate looks like this:


And the plot of the fft output looks like:


So in principle this looks ok.

However if I change the amplitude of the sine wave to y=(int)(200*(1+sin(x)));

The fft changes to:


Given that a fft is meant to transform from the time domain to the frequency domain the amplitude of the signal should not matter. So I'm getting someting fairly fundamental wrong here!

Anyone able to help out with some basics on using fftw to get the frequency (ies) of an imput signal?

Cheers,
~Al

[code] while ( x < 100 ){
y=(int)(100*(1+sin(x)));
in[z]=y;
fprintf(fp, "%f\n", in[z]);

Al,
should that second last line not read as int[z]=y and same for last part? Is that denoting an integer of [z]? Could it be this typo that's a problem?

Seanie.

[code] while ( x < 100 ){
y=(int)(100*(1+sin(x)));
in[z]=y;
fprintf(fp, "%f\n", in[z]);

Al,
should that second last line not read as int[z]=y and same for last part? Is that denoting an integer of [z]? Could it be this typo that's a problem?

Seanie.

no its an array called 'in'. But thanks for at least reading through the post anyway!

It's the best I could come up with!

Is it because you are not getting a 'proper' wave pattern that is wrong? t
The fact you increased 100 to 200 shows that your FFT wave is from almost 50 to almost 100. What (pictorially) are you hoping it would have changed when you 'doubled' it with the y=(int)(200*(1+sin(x)))

It's the best I could come up with!

Is it because you are not getting a 'proper' wave pattern that is wrong? t
The fact you increased 100 to 200 shows that your FFT wave is from almost 50 to almost 100. What (pictorially) are you hoping it would have changed when you 'doubled' it with the y=(int)(200*(1+sin(x)))

Have a peek at the images here:
www.numerit.com/samples/fft/doc.htm

The 'normal' graph has amplitude on the Y axis and Time on the X axis; thats the normal wave graph.

Applying a fourier transform changes you from the time domain to the frequency domain so you get the frequency on the X axis. (yes I know its on the Y in my plot, thats part of the problem!)

So for figure 3 in that URL we have a wave that is made up of 2 sine waves of different frequency and amplitude.

Figure 4 shows what the frequency of the waves is.

So you feed in a signal and the output from the transform tells you the component wave that make it up. So for my pulsar I feed in the counts and the pulsar frequency should show as a spike in the transform.

I (vaguely) remember Fourier transforms in college being involved with imaginary numbers not real ones so I suspect that my cop out of using the real2real functions is whats causing the problem. More readoing to do; just not at 2am

Cheers,
~Al

Well I think I'm out of my league now! It sounds like an interesting project though.