File:FFTexample16T.png: Difference between revisions

of a self–Fourier function ${\displaystyle A(x)=\exp(-x^{2}/2)}$, shown with black dots, to the result of the numerical evaluation of the the Fourier operator of array ${\displaystyle A}$, shown with blue. The discrete representation is performed with number of nodes ${\displaystyle n\!=\!16}$.

The function ${\displaystyle A}$ is represented at the mesh with nodes ${\displaystyle x_{n}={\sqrt {\frac {\pi }{8}}}(n\!-\!8)~,~n=0..15}$ in such a way that ${\displaystyle A_{n}=A(x_{n})}$

The FFT of array ${\displaystyle A}$ is performed with the routine

void fft(z_type *a, int N, int o) // Arry is FIRST!
{z_type u,w,t;  int i,j,k,l,e=1,L,p,q,m;
q=N/2;  p=2; for(m=1;p<N;m++) p*=2; 
p=N-1;  z_type y=z_type(0.,M_PI*o); j=0;
for(i=0;i<p;i++){ if(i<j) { t=a[j]; a[j]=a[i]; a[i]=t;}
                  k=q; while(k<=j){ j-=k; k/=2;} 
                  j+=k; }
for(l=0;l<m;l++)
{ L=e; e*=2; u=1.; w=exp(y/((double)L));
 for(j=0;j<L;j++)
 {  for(i=j;i<N;i+=e){k=i+L; t=a[k]*u; a[k]=a[i]-t; a[i]+=t;}
    u*=w;
} }
}

However, the transformation with the routine fafo.cin would return practically the same array, while the routine fafo is the discrete analogy of the Fourier operator and function ${\displaystyle A}$ is self Fourier.

The figure shows difference between FFT and the discrete implementation of the Fourier operator. For the self-Fourier Gaussian, the Fourier operator returns the same array, while the FFT returns the saw-like structure. For the physical applications, fafo.cin seems to be more suitable.

C++ generator of curves

// The picture is generated with the //File fafo.cin should be loaded to the working directory for the compilation of the code below.

#include<math.h>
#include<stdio.h>
#include <stdlib.h>
#include <complex>
using namespace std;
#define z_type complex<double>
#define Re(x)  x.real()
#define Im(x)  x.imag()
#define RI(x)  x.real(),x.imag()
#define DB double
#define DO(x,y) for(x=0;x<y;x++)

void ado(FILE *O, int X, int Y)
{       fprintf(O,"%c!PS-Adobe-2.0 EPSF-2.0\n",'%');
       fprintf(O,"%c%cBoundingBox: 0 0 %d %d\n",'%','%',X,Y);
       fprintf(O,"/M {moveto} bind def\n");
       fprintf(O,"/L {lineto} bind def\n");
       fprintf(O,"/S {stroke} bind def\n");
       fprintf(O,"/s {show newpath} bind def\n");
       fprintf(O,"/C {closepath} bind def\n");
       fprintf(O,"/F {fill} bind def\n");
       fprintf(O,"/o {.025 0 360 arc C F} bind def\n");
       fprintf(O,"/times-Roman findfont 20 scalefont setfont\n");
       fprintf(O,"/W {setlinewidth} bind def\n");
       fprintf(O,"/RGB {setrgbcolor} bind def\n");}
// #include "ado.cin"

#include"fafo.cin"

// DB F(DB x){DB u=x*x; return u*(-3.+u)*exp(-x*x/2.);}
DB F(DB x){ return exp(-x*x/2.);}

main(){z_type * a, *b, c; int j,m,n, N=16; FILE *o;
       double step=sqrt(2*M_PI/N),x,y,u;
       a=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
       b=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
//for(j=0;j<N;j++) { x=step*(j-N/2); u=x*x; a[j]=b[j]=(3.+u*(-6.+u))*exp(-x*x/2); }
for(j=0;j<N;j++) { x=step*(j-N/2); u=x*x; a[j]=b[j]=F(x); }
fft(b,N,1);
for(j=0;j<N;j++) printf("%2d %18.15f %18.15f  %18.15f %18.15f\n", j, RI(a[j]), RI(b[j])  );
o=fopen("FFTexample16.eps","w"); ado(o,1024,780); 
#define M(x,y) fprintf(o,"%6.4f %6.4f M\n",0.+x,0.+y);
#define L(x,y) fprintf(o,"%6.4f %6.4f L\n",0.+x,0.+y);
#define o(x,y) fprintf(o,"%6.4f %6.4f o\n",0.+x,0.+y);
fprintf(o,"522 340 translate 100 100 scale\n");
// M(-5,0) L(5,0) M(0,0) L(0,1) fprintf(o,".01 W S\n");
// M(-5,1) L(5,1) M(-5,-1) L(5,-1)
for(m=-5;m<6;m++) {M(m,-3) L(m,3)} fprintf(o,".004 W S\n");
for(m=-3;m<4;m++) {M(-5,m) L(5,m)} fprintf(o,".004 W S\n");
fprintf(o,"1 setlinejoin 1 setlinecap\n");
DB *X; X=(DB *) malloc((size_t)((N+1)*sizeof(DB))); DO(j,N){ x=step*(j-N/2); X[j]=x; }
DO(j,N){x=X[j]; M(x,0)L(x,.15)} fprintf(o,".01 W S\n");
DO(j,N){x=X[j];y=Re(a[j]); if(j==0)M(x,y)else L(x,y);} fprintf(o,".04 W 0 .4 1 RGB S\n");
DO(j,N){x=X[j];y=Re(b[j]); if(j==0)M(x,y)else L(x,y);} fprintf(o,".04 W 1 0 0 RGB S\n");
// DO(j,N){x=X[j];y=Im(b[j]); if(j==0)M(x,y)else L(x,y);} fprintf(o,".04 W 1 0 0 RGB S\n");
// DO(j,N){x=X[j];y=100.*(Re(b[j])-F(x)); if(j==0)M(x,y)else L(x,y);} fprintf(o,"0.007 W 0 0 .3 RGB S\n");
printf("X[0]=%9.5f step=%9.6f\n",X[0],step);
// DO(m,101){x=-5.+.1*m; y=F(x); if(m/2*2==m)M(x,y)else L(x,y);} fprintf(o,".01 W 0 0 0 RGB S\n");
fprintf(o,".01 W 0 0 0 RGB S\n");
DO(m,101){x=-5.+.1*m; y=F(x); o(x,y)}
fprintf(o,"showpage\n%cTrailer",'%'); fclose(o);
 system("epstopdf FFTexample16.eps");
 system(    "open FFTexample16.pdf"); //these 2 commands may be specific for macintosh
 getchar(); system("killall Preview");// if run at another operational system, may need to modify
 free(a);
 free(b);
 free(X);
}

The image is generated in the following way.

The lines are drawn in the EPS format by the C++ code below. The result is concerted to PDF format.

The labels are added in the latex document below.

The result is concerted to the PNG format with default reaolution.

{{Image|FFTexample16T.png|right|350px|Add image caption here.}}