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!

8)

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

8)

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.

Hmm .... FFT was never my strong point ..... code is .... so I was rather hopping this would be a code question but since it's a maths question that's me right out of my league ... sorry Al!

Bart.

Al --

A DFT is reversible; therefore it must also encode the amplitude of each frequency component (as well as its phase). I'm nearly a complete neophyte on FFTs, so you'll have to take this with a grain of salt -- but I think your results are the nature of the beast.

-- Jeff.

Hi,

ok I was a bit off

The inputs are real alright, however the output array is complex, a+ib. To plot the output array you plot the magnitude of the value sqrt(a*a+b*b).

Now wheres that pulsar?

I was just going to post that Al!

I was just going to post that Al!

Guess the telapathy worked better

Hi,

Just to let you know that there was a point to this topic!

Say you have some data from RXTE, which consists of the exact timing (to tens of microseconds) of when photons are detected from a pulsar (plus noise, ambient light etc.). Now split all this into 'bins' of 0.01s, so we are counting the number of photons that arrive in a sequence of 0.01 seconds. The plot of that looks like this:



So we get between 0 and 9 counts per sample. Great but it doesn't even look like there is a pattern in there. Now we put that through the fft transform which will show us what frequencies are in there:



Now something looks obvious! The peak is at 38245, which means that it has 38245 cycles during out timespan, which equates to a source with a frequency of 19.814792 Hz

By looking at later observations to see how the pulsar frequency has changed we can work out the age and magnetic field of the object.

Cheers,
~Al

Nice report, Al.

-- Jeff.