Fast Fourier Transform (FFT) Algorithm in C: Code & Explanation

Posted on Mar 22, 2025 in Personal Image

Fast Fourier Transform (FFT) Algorithm in C

This document provides a C implementation of the Fast Fourier Transform (FFT) algorithm.

C Code Implementation

The following C code demonstrates the FFT algorithm:

#define PI 3.1415926535897
typedef struct
{
    float real;
    float imag;
} COMPLEX;

int main(void)
{
    int n, N;
    void fft(COMPLEX *, COMPLEX *, int);
    FILE *fp1, *fp2;
    COMPLEX *x, *X;

    fp1 = fopen("input.dat","r");
    fscanf(fp1, "%d", &N);
    x = (COMPLEX *) malloc(N * sizeof(COMPLEX));
    X = (COMPLEX *) malloc(N * sizeof(COMPLEX));
    for(n = 0; n < N; n++){
        fscanf(fp1, "%f", &x[n].real);
        x[n].imag = 0.0;
    }

    fft(x, X, N);

    fp2 = fopen("output.dat","w");
    fprintf(fp2, "%d\n", N);
    for(n = 0; n < N; n++){
        fprintf(fp2, "%f + j %f\n", X[n].real, X[n].imag);
    }
    fclose(fp2);

    return 0;
}

void fft(COMPLEX *x, COMPLEX *X, int N)
{
    int n, k;
    COMPLEX *f1, *f2, *F1, *F2;
    f1 = (COMPLEX *) malloc(N/2 * sizeof(COMPLEX));
    f2 = (COMPLEX *) malloc(N/2 * sizeof(COMPLEX));
    F1 = (COMPLEX *) malloc(N/2 * sizeof(COMPLEX));
    F2 = (COMPLEX *) malloc(N/2 * sizeof(COMPLEX));

    if(N > 2)
    {
        for(n = 0; n < N/2; n++)
        {
            f1[n].real = x[2 * n].real;
            f1[n].imag = x[2 * n].imag;
            f2[n].real = x[2 * n + 1].real;
            f2[n].imag = x[2 * n + 1].imag;
        }

        fft(f1, F1, N/2);
        fft(f2, F2, N/2);

        for(k = 0; k < N/2; k++)
        {
            X[k].real = F1[k].real + (F2[k].real * cos(2 * PI / N * k) + F2[k].imag * sin(2 * PI / N * k));
            X[k].imag = F1[k].imag + (F2[k].imag * cos(2 * PI / N * k) - F2[k].real * sin(2 * PI / N * k));
            X[k + N/2].real = F1[k].real - (F2[k].real * cos(2 * PI / N * k) + F2[k].imag * sin(2 * PI / N * k));
            X[k + N/2].imag = F1[k].imag - (F2[k].imag * cos(2 * PI / N * k) - F2[k].real * sin(2 * PI / N * k));
        }
        free(f1);
        free(f2);
        free(F1);
        free(F2);
    }

    else
    {
        X[0].real = x[0].real + x[1].real;
        X[0].imag = x[0].imag + x[1].imag;
        X[1].real = x[0].real - x[1].real;
        X[1].imag = x[0].imag - x[1].imag;
    }
}

Explanation

The code implements a recursive Cooley-Tukey FFT algorithm. It takes an array of complex numbers as input and returns its Discrete Fourier Transform. The algorithm works by recursively breaking down a DFT of size N into two DFTs of size N/2.

Key Components:

• COMPLEX struct: Defines a complex number with real and imaginary parts.
• fft function: The recursive function that performs the FFT.
• Input/Output: Reads input from “input.dat” and writes output to “output.dat”.