#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define		pi		3.14159265359
#define		datasize        1024

double COSINE[datasize], SINE[datasize];       /* Twiddle tables           */
/* RData[] is the real component, CData[] is the imaginary component */
double RData[datasize], CData[datasize];

/********************************************************/
/* MAKE_TWIDDLE_TABLE                                   */
/* - Generates a FFT twiddle table.                     */
/********************************************************/
void make_twiddle_table(int Inputs)
{
int    i;
int    size = (Inputs / 2);
double arg1, arg2;
   arg1 = (2.0 * pi / (double)Inputs);
   for (i = 1; i <= Inputs; i++) {
	   arg2 = (double)(i - size) * arg1;
	   COSINE[i-1] = cos(arg2);
	   SINE[i-1]   = sin(arg2);
   }
}

/****************************************/
/* FFT                                  */
/****************************************/
void fft(double RData[datasize], double CData[datasize], int Inputs)
{
int    i, z, zz, j, j2, j3, k, r, d, st;
int    nd2, nml, arg, twf;
double tempr, tempc, c, s;
   for (i = 0; i < Inputs; i+=2) {
	   RData[i] = -1.0 * RData[i];
	   CData[i] = -1.0 * CData[i];
   }
   if (Inputs ==  16) st = 4;
   if (Inputs ==  32) st = 5;
   if (Inputs ==  64) st = 6;
   if (Inputs == 128) st = 7;
   if (Inputs == 256) st = 8;
   if (Inputs == 512) st = 9;
   if (Inputs == 1024) st = 10;
   j = 1;
   for (z = 1; z < Inputs; z++) {
	   if (z < j) {
		  tempr      = RData[j-1];
		  RData[j-1] = RData[z-1];
		  RData[z-1] = tempr;
		  tempc      = CData[j-1];
		  CData[j-1] = CData[z-1];
		  CData[z-1] = tempc;
       }
	   k = (Inputs / 2);
	   while (k < j) {
		  j -= k;
		  k  = (k / 2);
       }
	   j += k;
   }
   for (z = 1; z <= st; z++) {
	   r   = 1;
	   for (i = 1; i <= z; i++) r = (r * 2);
	   d   = (r / 2);
	   arg = (Inputs / r);
	   for (zz = 1; zz <= d; zz++) {
		   twf = (Inputs/4-1) + (arg * (zz-1));
		   c   = COSINE[twf];
	   s   = SINE[twf];
		   k   = zz;
	   while (k <= Inputs) {
			  j2 = (k + d) - 1;
			  j3 = (k - 1);
			  tempr     = (RData[j2] * c) + (CData[j2] * s);
			  tempc     = (CData[j2] * c) - (RData[j2] * s);
			  RData[j2] = (RData[j3] - tempr);
	      CData[j2] = (CData[j3] - tempc);
			  RData[j3] = (RData[j3] + tempr);
	      CData[j3] = (CData[j3] + tempc);
	      k += r;
		   }
       }
   }
}

/****************************************/
/* MODULUS                              */
/* - Calculates modulus of data.        */
/****************************************/
void modulus(double rdata[datasize], double cdata[datasize], int Inputs)
{
int    i;
double temp;
	for (i = 0; i < Inputs; i++) {
		temp = sqrt(rdata[i] * rdata[i] + cdata[i] * cdata[i]);
		rdata[i] = temp;
		cdata[i] = 0.0;
	}
}

/********************************************************/
/* MAIN                                                 */
/********************************************************/
void main(void)
{
/* Pixels can be 32, 64, 128, 256, or 512 depends on data set size */

int k;
double fk;
int n;
int N = 32;
float kf = 2;
float f;
float fmax;
float fa;
float W;

/********************************************************/
printf ("\nFast-Fourier-Transformation eines Sinussignals.");
printf ("\nAnzahl der St�tzstellen: ");
scanf ("%d",&N);
printf ("Frequenz des Sinussignals [Hz]: ");
scanf ("%f",&f);
printf ("Maximale Frequenz des Nutzsignals [Hz] : ");
scanf ("%f",&fmax);
while (1) {
printf ("Abtastfrequenz [Hz]: ");
scanf ("%f",&fa);
if (fa < 2*fmax)
	{ printf ("Verletzung des Abtasttheorems nach Nyquist.\nBitte Eingabe wiederholen\n"); continue ;}
	break;
}
kf = fa / fmax;
W = 2.0 * pi * f;

/********************************************************/

for (n=0; n < datasize; n++)
  { RData[n] = (double) sin (W*n/fa);
    CData[n] = 0;
}
make_twiddle_table(N/2);

/* Call FFT routine with complex data */
	fft(RData, CData, N);

/* Call modulus routine - optional */
	modulus(RData, CData, N);

printf ("\n\n%d-Punkt FFT",N);
printf ("\nf = %g Hz\nH”chste Frequenz = %g Hz",f,fmax);
printf ("\nAbtastfrequenz = %g Hz\n\nFrequenzspektrum S(f):",fa);

for ( k=0; k <= N/2; k++)
    {  printf ( "\n%g Hz = %f",(double)k*fa/(N),2*RData[k]/N);


}
printf ("\n\nEnde");

}