168 lines
4.4 KiB
C
168 lines
4.4 KiB
C
#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");
|
|
|
|
}
|
|
|