## FFTW

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

*FFTW* was created by *albertw*

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

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- Seanie_Morris
- Offline
- Administrator

- Posts: 9637
- Thank you received: 544

### Replied by *Seanie_Morris* on topic *Re: FFTW*

[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

*in*and same for last part? Is that denoting an integer of [z]? Could it be this typo that's a problem?

**t**[z]=ySeanie.

Radio Presenter (Midlands 103), Space Enthusiast, Astronomy Outreach Co-ordinator.

Former IFAS Chairperson and Secretary.

Please Log in or Create an account to join the conversation.

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

### Replied by *albertw* on topic *Re: FFTW*

[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 asinand same for last part? Is that denoting an integer of [z]? Could it be this typo that's a problem?t[z]=y

Seanie.

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

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- Seanie_Morris
- Offline
- Administrator

- Posts: 9637
- Thank you received: 544

### Replied by *Seanie_Morris* on topic *Re: FFTW*

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

Radio Presenter (Midlands 103), Space Enthusiast, Astronomy Outreach Co-ordinator.

Former IFAS Chairperson and Secretary.

Please Log in or Create an account to join the conversation.

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

### Replied by *albertw* on topic *Re: FFTW*

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

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- Seanie_Morris
- Offline
- Administrator

- Posts: 9637
- Thank you received: 544

### Replied by *Seanie_Morris* on topic *Re: FFTW*

Radio Presenter (Midlands 103), Space Enthusiast, Astronomy Outreach Co-ordinator.

Former IFAS Chairperson and Secretary.

Please Log in or Create an account to join the conversation.

- voyager
- Offline
- Super Giant

- Posts: 3663
- Thank you received: 2

### Replied by *voyager* on topic *Re: FFTW*

Bart.

Please Log in or Create an account to join the conversation.

- jeyjey
- Offline
- Red Giant

- Posts: 757
- Thank you received: 10

### Replied by *jeyjey* on topic *Re: FFTW*

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.

Solarscope SF70 / TV Pronto / AP400QMD Coronado SolarMax40 DS / Bogen 055+3130

APM MC1610 / Tak FC-125 / AP1200GTO Tak Mewlon 250 / AP600EGTO

Please Log in or Create an account to join the conversation.

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

### Replied by *albertw* on topic *Re: FFTW*

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?

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- dmcdona
- Offline
- Administrator

- Posts: 4557
- Thank you received: 0

### Replied by *dmcdona* on topic *Re: FFTW*

Please Log in or Create an account to join the conversation.

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

### Replied by *albertw* on topic *Re: FFTW*

I was just going to post that Al!

Guess the telapathy worked better

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- albertw
- Topic Author
- Offline
- IFAS Secretary

- Posts: 4173
- Thank you received: 181

### Replied by *albertw* on topic *Re: FFTW*

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

Chairperson, International Dark Sky Association - Irish Section

www.darksky.ie/

Please Log in or Create an account to join the conversation.

- jeyjey
- Offline
- Red Giant

- Posts: 757
- Thank you received: 10

### Replied by *jeyjey* on topic *Re: FFTW*

-- Jeff.

Solarscope SF70 / TV Pronto / AP400QMD Coronado SolarMax40 DS / Bogen 055+3130

APM MC1610 / Tak FC-125 / AP1200GTO Tak Mewlon 250 / AP600EGTO

Please Log in or Create an account to join the conversation.