Implementation of fast Fourier transform of C data structure algorithm

Implementation of fast Fourier transform of C data structure algorithm

This example will realize the two-dimensional fast Fourier transform, and meanwhile, it will also learn the basic operation of matrix and the basic operation of complex number in c language, and review the dynamic memory allocation, file operation, function call of structure pointer and so on.

For a long time, 1 direct Fourier transform is many fields, such as a linear system, optical, probability theory, quantum physics, antenna, digital image processing and signal analysis, such as a basic analysis tool, but even at a staggering rate of computer use to calculate the discrete Fourier transform the time also is often difficult to accept, thus led to the fast Fourier transform (FFT).

This example will carry out the forward and reverse fast Fourier transform on a 2 - dimensional array. In the case of positive Fourier transform, dfft() function first calls fft() to transform the array by rows, and then calls fft() to transform the result by columns. At this point, fast Fourier transform is completed, and then calls rdfft() to perform inverse Fourier transform. If the program is designed correctly, the result of the transformation should be the same as the original array.

Example code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
 
struct COMPLEX
{
  float re;
  float im;
} cplx , * Hfield , * S , * R , * w;
 
int n , m;
int ln , lm;
 
void initiate ();
void dfft ();
void rdfft ();
void showresult ();
 
void fft (int l , int k);
int reverse (int t , int k);
void W (int l);
int loop (int l);
void conjugate ();
 
void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);
void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);
void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);
struct COMPLEX * Hread(int i , int j);
void Hwrite (int i , int j , struct COMPLEX x);
 
void main ()
{
  initiate ();
  printf("\n The original data :\n");
  showresult();
  getchar ();
  dfft ();
  printf("\n The result of rapid compound leaf transformation :\n");
  showresult ();
  getchar ();
  rdfft ();
  printf("\n The inverse transformation of a fast compound leaf :\n");
  showresult ();
  getchar ();
  free (Hfield);
}
 
void initiate ()
{// Program initialization operations, including allocation of memory, reading in the data to be processed, display, and so on 
  FILE * df;
   
  df = fopen ("data.txt" , "r");
  fscanf (df , "%5d" , &n);
  fscanf (df , "%5d" , &m);
  if ((ln = loop (n)) == -1)
  {
    printf ("  The number of columns is not 2 Omega to an integer power  ");
    exit (1);
  }
  if ((lm = loop (m)) == -1)
  {
    printf ("  Number of lines is not 2 Omega to an integer power  ");
    exit (1);
  }
  Hfield = (struct COMPLEX *) malloc (n * m * sizeof (cplx));
  if (fread (Hfield , sizeof (cplx) , m * n , df) != (unsigned) (m * n))
  {
    if (feof (df)) printf (" Premature end of file ");
    else printf (" File read error ");
  }
  fclose (df);
}
 
void dfft ()
{// for 2 D rapid compound leaf transformation 
  int i , j; 
  int l , k;
   
  l = n;
  k = ln;
  w = (struct COMPLEX *) calloc (l , sizeof (cplx));
  R = (struct COMPLEX *) calloc (l , sizeof (cplx));
  S = (struct COMPLEX *) calloc (l , sizeof(cplx));
  W (l);
  for ( i = 0 ; i < m ; i++ )
  {// Fast compound leaf transformation by row 
    for (j = 0 ; j < n ; j++)
    {      
      S[j].re = Hread (i , j)->re;
      S[j].im = Hread (i , j)->im;
    }
    fft(l , k);
    for (j = 0 ; j < n ; j++)
      Hwrite (i , j , R[j]);
  }
  free (R);
  free (S);
  free (w);
   
  l = m;
  k = lm;
  w = (struct COMPLEX *) calloc (l , sizeof (cplx));
  R = (struct COMPLEX *) calloc (l , sizeof (cplx));
  S = (struct COMPLEX *) calloc (l , sizeof (cplx));
  W (l);
  for (i = 0 ; i < n ; i++)
  {// Fast compound leaf transformation by column 
    for(j = 0 ; j < m ; j++)
    {
      S[j].re = Hread(j , i)->re;
      S[j].im = Hread(j , i)->im;
    }
    fft(l , k);
    for (j = 0 ; j < m ; j++)
      Hwrite (j , i , R[j]);
  }
  free (R);
  free (S);
  free (w);
}
 
void rdfft ()
{
  conjugate ();
  dfft ();
  conjugate ();
}
 
void showresult ()
{
  int i , j;
  for (i = 0 ; i < m ; i++)
  {
    printf ( " \n The first %d line \n " , i);
    for (j = 0 ; j < n ; j++)
    {
      if (j % 4 == 0) printf (" \n ");
      printf(" (%5.2f,%5.2fi) " , Hread (i , j)->re , Hread (i , j)->im);
    }
  }
}
 
void fft (int l , int k)
{
  int i , j , s , nv , t;
  float c;
  struct COMPLEX mp , r;
  nv = l;
  c = (float) l;
  c = pow (c , 0.5);
  for (i = 0 ; i < k ; i++)
  {
    for (t = 0 ; t < l ; t += nv)
    {
      for (j = 0 ; j < nv / 2 ; j++)
      {
        s = (t + j) >> (k - i -1);
        s = reverse(s , k);
        r.re = S[t + j].re;
        r.im = S[t + j].im;
        mul (&w[s] , &S[t + j + nv / 2] , &mp);///////// Explains the difference between passing a pointer to a structure and the structure itself 
        add (&r , &mp , &S[t + j]);
        sub (&r , &mp , &S[t + j + nv / 2]);        
      }
    }
    nv = nv >> 1;   
  }
 
  for (i = 0 ; i < l ; i++)
  {
    j = reverse(i , k);
    R[j].re = S[i].re / c;
    R[j].im = S[i].im / c;
  }
}
 
int reverse (int t , int k)
{
  int i , x , y;
  y = 0;
  for (i = 0 ; i < k ; i++)
  {
    x = t & 1;
    t = t >> 1;
    y = (y << 1) + x;   
  }
  return y;
}
 
void W (int l)
{
  int i;
  float c , a;
  c = (float) l;
  c = 2 * PI / c;
  for (i = 0 ; i < l ; i++)
  {    
    a = (float) i;
    w[i].re = (float) cos(a * c);
     
    w[i].im = -(float) sin(a * c);
  }
}
 
int loop (int l)
{// Verify that the input data is 2 To an integer power of 2 Base is the number of digits in time 
  int i , m;
  if (l != 0)
  {
    for (i = 1 ; i < 32 ; i++)
    {
      m = l >> i;
      if (m == 0)
        break;
    }
    if (l == (1 << (i - 1)))
      return i - 1;
  }
  return -1;
}
 
void conjugate ()
{// Find the conjugate of the complex matrix 
  int i , j;
  for (i = 0 ; i < m ; i++)
  {
    for (j = 0 ; j < n ; j++)
    {
      Hread (i , j)->im *= -1;
    }
  }
}
 
struct COMPLEX * Hread (int i , int j)
{// Returns as a read matrix Hfield A pointer to the specified location in 
  return (Hfield + i * n + j);
}
 
void Hwrite (int i , int j , struct COMPLEX x)
{// Construct complex Numbers as a matrix x Write to specified Hfield position 
  (Hfield + i * n + j)->re = x.re;
  (Hfield + i * n + j)->im = x.im;
}
 
void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z)
{// Define complex addition 
  z->re = x->re + y->re;
  z->im = x->im + y->im; 
}
 
void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z)
{// Define the subtraction of complex Numbers 
  z->re = x->re - y->re;
  z->im = x->im - y->im;
}
 
void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z)
{// Define complex multiplication 
  z->re = (x->re) * (y->re) - (x->im) * (y->im);
  z->im = (x->im) * (y->re) + (x->re) * (y->im);
}

Thank you for reading, I hope to help you, thank you for your support of this site!