FFT Fast Fourier Transform algorithm suitable for single-chip microcomputer, 51 single-chip microcomputer can be used

Source code

FFT.c

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

Fast Fourier Transform C package

Function Introduction: This package is a general fast Fourier transform C language function with strong portability. The following parts are not dependent on

This package uses a union to represent a complex number. The input is a complex number in natural order.

The input real number can make the complex imaginary part 0, and the output is the natural order of the FFT transformation.

Complex numbers. This package can call the create_sin_tab() function to create a sine function table during initialization.

In the future, the table lookup method can be used to calculate the sin and cos operations that take a long time to calculate, which can speed up the calculation speed.

Compared with Ver1.1, Ver1.2 only creates 1/4 of the sine wave sampling value when creating the sine table.

In comparison, it saves FFT_N/4 storage space

Instructions for use: To use this function, you only need to change the value of the macro definition FFT_N to change the number of points.

It should be 2 to the power of N. If this condition is not met, 0 should be added at the end.

cos value, you should call create_sin_tab() function to create a sine table before calling FFT function

Function call: FFT(Compx);

Author: Ji Shuaihu

Date: 2010-2-20

Version: Ver1.2

references:

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

#include

#include "FFT.h"

struct compx Compx[FFT_N] = {0}; //FFT input and output: start from Compx[0] and define it according to the size

double SIN_TAB[FFT_N / 4 + 1]; //Define the storage space for the sine table

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

Function prototype: struct compx EE (struct compx b1, struct compx b2)

Function: Multiply two complex numbers

Input parameters: two complex numbers a and b defined as a union

Output parameter: the product of a and b, output in the form of a union

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

struct compx EE (struct compx a, struct compx b)

{

struct compx c;

c.real = a.real*b.real - a.imag*b.imag;

c.imag = a.real*b.imag + a.imag*b.real;

return(c);

}

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

Function prototype: void create_sin_tab(double *sin_t)

Function: Create a sine sampling table with the same number of sampling points as the Foley transform points

Input parameter: *sin_t array pointer storing sine table

Output parameters: None

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

void create_sin_tab(double *sin_t)

{

int i;

for (i = 0; i <= FFT_N / 4; i++)

sin_t[i] = sin(2 * PI*i / FFT_N);

}

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

Function prototype: void sin_tab(double pi)

Function: Calculate the sine value of a number using a table lookup method

Input parameter: pi is the radian value of the sine value to be calculated, ranging from 0 to 2*PI, and needs to be converted if it does not meet the requirements

Output parameter: sine value of input value pi

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

double sin_tab(double pi)

{

int n;

double a = 0;

n = (int)(pi*FFT_N / 2 / PI);

if (n >= 0 && n <= FFT_N / 4)

a = SIN_TAB[n];

else if (n>FFT_N / 4 && n {

n -= FFT_N / 4;

a = SIN_TAB[FFT_N / 4 - n];

}

else if (n >= FFT_N / 2 && n<3 * FFT_N / 4)

{

n -= FFT_N / 2;

a = -SIN_TAB[n];

}

else if (n >= 3 * FFT_N / 4 && n<3 * FFT_N)

{

n = FFT_N - n;

a = -SIN_TAB[n];

}

return a;

}

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

Function prototype: void cos_tab(double pi)

Function: Calculate the cosine value of a number using a table lookup method

Input parameter: pi is the radian value of the cosine value to be calculated, ranging from 0 to 2*PI, and needs to be converted if it does not meet the requirements

Output parameter: cosine value of input value pi

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

double cos_tab(double pi)

{

double a, pi2;

pi2 = pi + PI / 2;

if (pi2>2 * PI)

pi2 -= 2 * PI;

a = sin_tab(pi2);

return a;

}

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

Function prototype: void FFT(struct compx *xin)

Function: Perform fast Fourier transform (FFT) on the input complex array

Input parameter: *xin first address pointer of complex structure group, struct type

Output parameters: None

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

void FFT(struct compx *xin)

{

register int f, m, nv2, nm1, i, k, l, j = 0;

struct compx u, w, t;

nv2 = FFT_N / 2; //Address operation, that is, changing the natural order into the reverse order, using the Reid algorithm

nm1 = FFT_N - 1;

for (i = 0; i < nm1; ++i)

{

if (i < j) //If i {

t = xin[j];

xin[j] = xin[i];

xin[i] = t;

}

k = nv2; //Find the next inversion sequence of j

while (k <= j) //If k<=j, it means the highest bit of j is 1

{

j = j - k; //Change the highest bit to 0

k = k / 2; //k/2, compare the second highest bit, and so on, compare one by one until a bit is 0

}

j = j + k; //Change 0 to 1

}

{

int le, lei, ip; //FFT operation core, use butterfly operation to complete FFT operation

f = FFT_N;

for (l = 1; (f = f / 2) != 1; ++l); //Calculate the value of l, that is, calculate the butterfly series

for (m = 1; m <= l; m++) // Control the number of butterfly knots

{   

//m represents the mth butterfly, l is the total number of butterfly levels l=log（2）N

le = 2 << (m - 1); //le butterfly knot distance, that is, the distance between the butterfly knot of the mth level butterfly and the point le

lei = le / 2; //The distance between two points in the same butterfly knot

u.real = 1.0; //u is the butterfly knot operation coefficient, the initial value is 1

u.imag = 0.0;

w.real = cos_tab(PI / lei); //w is the coefficient quotient, that is, the quotient of the current coefficient and the previous coefficient

w.imag = -sin_tab(PI / lei);

for (j = 0; j <= lei - 1; j++) //Control the calculation of different types of butterfly knots, that is, the calculation of butterfly knots with different coefficients

{

for (i = j; i <= FFT_N - 1; i = i + le) //Control the same butterfly knot operation, that is, calculate the butterfly knot with the same coefficient

{

ip = i + lei; //i, ip represent the two nodes participating in the butterfly operation

t = EE(xin[ip], u); //Butterfly operation, see the formula for details

xin[ip].real = xin[i].real - t.real;

xin[ip].imag = xin[i].imag - t.imag;

xin[i].real = xin[i].real + t.real;

xin[i].imag = xin[i].imag + t.imag;

}

u = EE(u, w); //Change the coefficient and perform the next butterfly operation

}

}

}

}

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

Function prototype: void Get_Result(struct compx *xin, double sample_frequency)

Function: Calculate the modulus of the transformed result, store the real part of the complex number, store the frequency in the imaginary part of the complex number, and the valid data is the first FFT_N/2 numbers

Input parameters: *xin first address pointer of complex structure group, struct type, sample_frequency: sampling frequency

Output parameters: None

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

void Get_Result(struct compx *xin, double sample_frequency)

{

int i = 0;

for (i = 0; i < FFT_N / 2; ++i)

{          

//Find the modulus of the transformed result and store it in the real part of the complex number

xin[i].real = sqrt(xin[i].real*xin[i].real + xin[i].imag*xin[i].imag) / (FFT_N >> (i != 0));

xin[i].imag = i * sample_frequency / FFT_N;

}

}

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

Function prototype: void Refresh_Data(struct compx *xin, double wave_data)

Function: Update data

Input parameters: *xin first address pointer of complex structure group, struct type, id: label, wave_data: value of a point

Output parameters: None

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

void Refresh_Data(struct compx *xin, int id, double wave_data)

{

xin[id].real = wave_data;

xin[id].imag = 0;

}

FFT.h

#ifndef FFT_H

#define FFT_H

#define FFT_N 16 //Define the number of points for Fourier transform

#define PI 3.14159265358979323846264338327950288419717 //define the value of pi

struct compx { double real, imag; }; //define a complex number structure

extern struct compx Compx[]; //FFT input and output: start from Compx[0] and define it according to the size

extern double SIN_TAB[]; //sine signal table

extern void Refresh_Data(struct compx *xin, int id, double wave_data);

extern void create_sin_tab(double *sin_t);

extern void FFT(struct compx *xin);

extern void Get_Result(struct compx *xin, double sample_frequency);

#endif

Instructions

Modify FFT_N 16 in FFT.h to define the number of Fourier transform points. After the FFT transform is normalized, the entire result table is symmetrical except for the DC component of 0Hz, that is, if the number of points is 16, we only need to look at the first 8 points. Therefore, the number of points of this Fourier transform can be determined according to the number of pixels in the long direction of your screen. For example, a 128x64 screen can consider using a 256-point FFT. Here I use an 8x8 LED dot matrix screen for display, so I use a 16-point FFT.

Before the operation, you need to call create_sin_tab(SIN_TAB) once to create a sine signal table. Then you will use the table lookup method to calculate the sine value to speed up the calculation.

After filling in data using the Refresh_Data (Compx, i, wave_data) function

Call FFT(Compx) to complete the transformation

Use the Get_Result (Compx, Sample_Frequency) function to find the modulus of the transformed result, store the real part of the complex number, and store the frequency in the imaginary part of the complex number. The valid data is the first FFT_N/2 numbers.