Objective: C program to implement Fourier transformation
Concept:A fast Fourier transform (FFT) is an efficient algorithm to compute the discrete Fourier transform (DFT) and its inverse. There are many distinct FFT algorithms involving a wide range of mathematics, from simple complex-number arithmetic to group theory and number theory.A DFT decomposes a sequence of values into components of different frequencies. This operation is useful in many fields but computing it directly
from the definition is often too slow to be practical. An FFT is a way to compute the same result more quickly: computing a DFT of N points in the obvious way, using the definition, takes O(N 2) arithmetical operations, while an FFT can compute the same result in only O(N log N) operations. The difference in speed can be substantial, especially for long data sets where N may be in the thousands or millions—in practice, the computation time can be reduced by several orders of magnitude in such cases, and the improvement is roughly proportional to N/log(N). This huge improvement made many DFT-based algorithms practical; FFTs are of great importance to a wide variety of applications, from digital signal processing and solving partial differential equations to algorithms for quick multiplication of large integers. The most well known FFT algorithms depend upon the factorization of N, but (contrary to popular misconception) there are FFTs with O(N log N) complexity for all N, even for prime N. Many FFT algorithms only depend on the fact that is an Nth primitive root of unity, and thus can be applied to analogous transforms over any finite field, such as number-theoretic transforms. Since the inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it
/*----------------------------------------------------------------------------
fft.c - fast Fourier transform and its inverse (both recursively)
----------------------------------------------------------------------------*/
/******************************************************************************
* This file defines a C function fft that, by calling another function *
* fft_rec (also defined), calculates an FFT recursively. Usage: *
* fft(N, x, X); *
* Parameters: *
* N: number of points in FFT (must equal 2^n for some integer n >= 1) *
* x: pointer to N time-domain samples given in rectangular form (Re x, *
* Im x) *
* X: pointer to N frequency-domain samples calculated in rectangular form *
* (Re X, Im X) *
* Similarly, a function ifft with the same parameters is defined that *
* calculates an inverse FFT (IFFT) recursively. Usage: *
* ifft(N, x, X); *
* Here, N and X are given, and x is calculated. *
******************************************************************************/
#include <stdlib.h>
#include <math.h>
/* macros */
#define TWO_PI (6.2831853071795864769252867665590057683943L)
/* function prototypes */
void fft(int N, double (*x)[2], double (*X)[2]);
void fft_rec(int N, int offset, int delta,
double (*x)[2], double (*X)[2], double (*XX)[2]);
void ifft(int N, double (*x)[2], double (*X)[2]);
/* FFT */
void fft(int N, double (*x)[2], double (*X)[2])
{
/* Declare a pointer to scratch space. */
double (*XX)[2] = malloc(2 * N * sizeof(double));
/* Calculate FFT by a recursion. */
fft_rec(N, 0, 1, x, X, XX);
/* Free memory. */
free(XX);
}
/* FFT recursion */
void fft_rec(int N, int offset, int delta,
double (*x)[2], double (*X)[2], double (*XX)[2])
{
int N2 = N/2; /* half the number of points in FFT */
int k; /* generic index */
double cs, sn; /* cosine and sine */
int k00, k01, k10, k11; /* indices for butterflies */
double tmp0, tmp1; /* temporary storage */
if(N != 2) /* Perform recursive step. */
{
/* Calculate two (N/2)-point DFT's. */
fft_rec(N2, offset, 2*delta, x, XX, X);
fft_rec(N2, offset+delta, 2*delta, x, XX, X);
/* Combine the two (N/2)-point DFT's into one N-point DFT. */
for(k=0; k
{
k00 = offset + k*delta; k01 = k00 + N2*delta;
k10 = offset + 2*k*delta; k11 = k10 + delta;
cs = cos(TWO_PI*k/(double)N); sn = sin(TWO_PI*k/(double)N);
tmp0 = cs * XX[k11][0] + sn * XX[k11][1];
tmp1 = cs * XX[k11][1] - sn * XX[k11][0];
X[k01][0] = XX[k10][0] - tmp0;
X[k01][1] = XX[k10][1] - tmp1;
X[k00][0] = XX[k10][0] + tmp0;
X[k00][1] = XX[k10][1] + tmp1;
}
}
else /* Perform 2-point DFT. */
{
k00 = offset; k01 = k00 + delta;
X[k01][0] = x[k00][0] - x[k01][0];
X[k01][1] = x[k00][1] - x[k01][1];
X[k00][0] = x[k00][0] + x[k01][0];
X[k00][1] = x[k00][1] + x[k01][1];
}
}
/* IFFT */
void ifft(int N, double (*x)[2], double (*X)[2])
{
int N2 = N/2; /* half the number of points in IFFT */
int i; /* generic index */
double tmp0, tmp1; /* temporary storage */
/* Calculate IFFT via reciprocity property of DFT. */
fft(N, X, x);
x[0][0] = x[0][0]/N; x[0][1] = x[0][1]/N;
x[N2][0] = x[N2][0]/N; x[N2][1] = x[N2][1]/N;
for(i=1; i
{
tmp0 = x[i][0]/N; tmp1 = x[i][1]/N;
x[i][0] = x[N-i][0]/N; x[i][1] = x[N-i][1]/N;
x[N-i][0] = tmp0; x[N-i][1] = tmp1;
}
}
/******************************************************************************
* This program demonstrates how to use the file fft.c to calculate an FFT *
* of given time-domain samples, as well as to calculate an inverse FFT *
* (IFFT) of given frequency-domain samples. First, N complex-valued time- *
* domain samples x, in rectangular form (Re x, Im x), are read from a *
* specified file; the 2N values are assumed to be separated by whitespace. *
* Then, an N-point FFT of these samples is found by calling the function *
* fft, thereby yielding N complex-valued frequency-domain samples X in *
* rectangular form (Re X, Im X). Next, an N-point IFFT of these samples is *
* is found by calling the function ifft, thereby recovering the original *
* samples x. Finally, the calculated samples X are saved to a specified *
* file, if desired. *
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "fft.c"
int main()
{
int i; /* generic index */
char file[FILENAME_MAX]; /* name of data file */
int N; /* number of points in FFT */
double (*x)[2]; /* pointer to time-domain samples */
double (*X)[2]; /* pointer to frequency-domain samples */
double dummy; /* scratch variable */
FILE *fp; /* file pointer */
/* Get name of input file of time-domain samples x. */
printf("Input file for time-domain samples x(n)? ");
scanf("%s", file);
/* Read through entire file to get number N of points in FFT. */
if(!(fp = fopen(file, "r")))
{
printf(" File \'%s\' could not be opened!", file);
exit(EXIT_FAILURE);
}
N=0;
while(fscanf(fp, "%lg%lg", &dummy, &dummy) == 2) N++;
printf("N = %d", N);
/* Check that N = 2^n for some integer n >= 1. */
if(N >= 2)
{
i = N;
while(i==2*(i/2)) i = i/2; /* While i is even, factor out a 2. */
} /* For N >=2, we now have N = 2^n iff i = 1. */
if(N < 2 || i != 1)
{
printf(", which does not equal 2^n for an integer n >= 1.");
exit(EXIT_FAILURE);
}
/* Allocate time- and frequency-domain memory. */
x = malloc(2 * N * sizeof(double));
X = malloc(2 * N * sizeof(double));
/* Get time-domain samples. */
rewind(fp);
for(i=0; i<N; i++) fscanf(fp, "%lg%lg", &x[i][0], &x[i][1]);
fclose(fp);
/* Calculate FFT. */
fft(N, x, X);
/* Print time-domain samples and resulting frequency-domain samples. */
printf("\nx(n):");
for(i=0; i<N; i++) printf("\n n=%d: %12f %12f", i, x[i][0], x[i][1]);
printf("\nX(k):");
for(i=0; i<N; i++) printf("\n k=%d: %12f %12f", i, X[i][0], X[i][1]);
/* Clear time-domain samples and calculate IFFT. */
for(i=0; i<N; i++) x[i][0] = x[i][1] = 0;
ifft(N, x, X);
/* Print recovered time-domain samples. */
printf("\nx(n):");
for(i=0; i<N; i++) printf("\n n=%d: %12f %12f", i, x[i][0], x[i][1]);
/* Write frequency-domain samples X to a file, if desired. */
printf("\nOutput file for frequency-domain samples X(k)?"
"\n (if none, abort program): ");
scanf("%s", file);
if(!(fp = fopen(file, "w")))
{
printf(" File \'%s\' could not be opened!", file);
exit(EXIT_FAILURE);
}
for(i=0; i<N; i++) fprintf(fp, "%23.15e %23.15e\n", X[i][0], X[i][1]);
fclose(fp);
printf("Samples X(k) were written to file %s.", file);
/* Free memory. */
free(x);
free(X);
return 0;
}
/*=========================================================================== === * Program output (example) *=========================================================================== === * Input file for time-domain samples x(n)? data.txt * N = 8 * x(n): * n=0: 3.600000 2.600000 * n=1: 2.900000 6.300000 * n=2: 5.600000 4.000000 * n=3: 4.800000 9.100000 * n=4: 3.300000 0.400000 * n=5: 5.900000 4.800000 * n=6: 5.000000 2.600000 * n=7: 4.300000 4.100000 * X(k): * k=0: 35.400000 33.900000 * k=1: 3.821320 0.892893 * k=2: -5.800000 -3.300000 * k=3: 5.971068 7.042641 * k=4: -0.400000 -14.700000 * k=5: -0.421320 2.307107 * k=6: -1.600000 -3.900000 * k=7: -8.171068 -1.442641 * x(n): * n=0: 3.600000 2.600000 * n=1: 2.900000 6.300000 * n=2: 5.600000 4.000000 * n=3: 4.800000 9.100000 * n=4: 3.300000 0.400000 * n=5: 5.900000 4.800000 * n=6: 5.000000 2.600000 * n=7: 4.300000 4.100000 * Output file for frequency-domain samples X(k)? * (if none, abort program): X.txt * Samples X(k) were written to file X.txt. */