[C++] Fast Fourier Transform algorithm does not work properly
-
Hello, there is a pseudocode for recursive radix 2 DIT Fast Fourier Transform(page 11,12) at this link. I have written a code in c++ which does the same, but it doesn't work properly. For example when I take 1024 datapoints of sin(0.1*x) where x={0,1,...1023}, I get this. Reverse transform inserted for this gives initial function multiplied by 1024, so it works. What could be wrong? I tried with value dir and dir/2 inside "Shift" function too. Here is the code: Notation: E2-even part, E3-odd part I work without complex numbers library, so there are both "_re" and "_im" tables.
void FFT(double *tab_re, double *tab_im, double *results_re, double *results_im, int n, int dir) //n-data size
{
if(n==1)
{results_re[0]=tab_re[0]; //real values
results_im[0]=tab_im[0]; //imaginary
return;}int nh=n/2;
double *E2_re=new double [nh];
double *E3_re=new double [nh];
double *E2_im=new double [nh];
double *E3_im=new double [nh];double *F2_re=new double [nh];
double *F3_re=new double [nh];
double *F2_im=new double [nh];
double *F3_im=new double [nh];for(int k=0;k
-
Hello, there is a pseudocode for recursive radix 2 DIT Fast Fourier Transform(page 11,12) at this link. I have written a code in c++ which does the same, but it doesn't work properly. For example when I take 1024 datapoints of sin(0.1*x) where x={0,1,...1023}, I get this. Reverse transform inserted for this gives initial function multiplied by 1024, so it works. What could be wrong? I tried with value dir and dir/2 inside "Shift" function too. Here is the code: Notation: E2-even part, E3-odd part I work without complex numbers library, so there are both "_re" and "_im" tables.
void FFT(double *tab_re, double *tab_im, double *results_re, double *results_im, int n, int dir) //n-data size
{
if(n==1)
{results_re[0]=tab_re[0]; //real values
results_im[0]=tab_im[0]; //imaginary
return;}int nh=n/2;
double *E2_re=new double [nh];
double *E3_re=new double [nh];
double *E2_im=new double [nh];
double *E3_im=new double [nh];double *F2_re=new double [nh];
double *F3_re=new double [nh];
double *F2_im=new double [nh];
double *F3_im=new double [nh];for(int k=0;k
If you get the original function back after the reverse, you have probably coded the transform correctly. The problem is likely to be that your sampling rate is not a multiple of the basic frequency of sin(0.1*x). This means that you get a large FFT component close to the frequency, but you also get "noise" throughout the range, due to the fact that the residues do not cancel out. I would try taking an FFT of sin(x), sin(2*x), etc., and see if these functions give the expected values.
Freedom is the freedom to say that two plus two make four. If that is granted, all else follows. -- 6079 Smith W.
-
If you get the original function back after the reverse, you have probably coded the transform correctly. The problem is likely to be that your sampling rate is not a multiple of the basic frequency of sin(0.1*x). This means that you get a large FFT component close to the frequency, but you also get "noise" throughout the range, due to the fact that the residues do not cancel out. I would try taking an FFT of sin(x), sin(2*x), etc., and see if these functions give the expected values.
Freedom is the freedom to say that two plus two make four. If that is granted, all else follows. -- 6079 Smith W.
Thank You for the hint. I have tried with sin(nx), sin(2pi*x/n) and sin(2pi*x*n) too, but results were similar. The most promising result was for the option with sin(2pi*x*n): here is a result with n=10. Values for x>270 are 0. I used decimation in frequency FFT instead of decimation in time version, but it still doesn't show good values.
-
Hello, there is a pseudocode for recursive radix 2 DIT Fast Fourier Transform(page 11,12) at this link. I have written a code in c++ which does the same, but it doesn't work properly. For example when I take 1024 datapoints of sin(0.1*x) where x={0,1,...1023}, I get this. Reverse transform inserted for this gives initial function multiplied by 1024, so it works. What could be wrong? I tried with value dir and dir/2 inside "Shift" function too. Here is the code: Notation: E2-even part, E3-odd part I work without complex numbers library, so there are both "_re" and "_im" tables.
void FFT(double *tab_re, double *tab_im, double *results_re, double *results_im, int n, int dir) //n-data size
{
if(n==1)
{results_re[0]=tab_re[0]; //real values
results_im[0]=tab_im[0]; //imaginary
return;}int nh=n/2;
double *E2_re=new double [nh];
double *E3_re=new double [nh];
double *E2_im=new double [nh];
double *E3_im=new double [nh];double *F2_re=new double [nh];
double *F3_re=new double [nh];
double *F2_im=new double [nh];
double *F3_im=new double [nh];for(int k=0;k
You lost me with your translation of the code .. page 11 & 12 have the exp function which doesn't appear in your code?
In vino veritas
-
You lost me with your translation of the code .. page 11 & 12 have the exp function which doesn't appear in your code?
In vino veritas
The exp function is inside function called "Shift". As I wrote earlier, I don't use complex numbers library, so there are both tables for real and imagery parts. I used this relation: exp(i*phi)=cos(phi)+i*sin(phi).
-
Hello, there is a pseudocode for recursive radix 2 DIT Fast Fourier Transform(page 11,12) at this link. I have written a code in c++ which does the same, but it doesn't work properly. For example when I take 1024 datapoints of sin(0.1*x) where x={0,1,...1023}, I get this. Reverse transform inserted for this gives initial function multiplied by 1024, so it works. What could be wrong? I tried with value dir and dir/2 inside "Shift" function too. Here is the code: Notation: E2-even part, E3-odd part I work without complex numbers library, so there are both "_re" and "_im" tables.
void FFT(double *tab_re, double *tab_im, double *results_re, double *results_im, int n, int dir) //n-data size
{
if(n==1)
{results_re[0]=tab_re[0]; //real values
results_im[0]=tab_im[0]; //imaginary
return;}int nh=n/2;
double *E2_re=new double [nh];
double *E3_re=new double [nh];
double *E2_im=new double [nh];
double *E3_im=new double [nh];double *F2_re=new double [nh];
double *F3_re=new double [nh];
double *F2_im=new double [nh];
double *F3_im=new double [nh];for(int k=0;k
Bonjour, //hi essayez avec : //try with arg=v*pi*i/n; k++); //two time i++);