- Posts: 4173
- Thank you received: 181
FFTW
- albertw
- Topic Author
- Offline
- IFAS Secretary
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
Albert White MSc FRAS
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: 9640
- Thank you received: 547
[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.
Midlands Astronomy Club.
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
[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!
Albert White MSc FRAS
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: 9640
- Thank you received: 547
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)))
Midlands Astronomy Club.
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
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
Albert White MSc FRAS
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: 9640
- Thank you received: 547
Midlands Astronomy Club.
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.