A thorough study of C-language algorithm for fast Fourier transform FFT

A thorough study of C-language algorithm for fast Fourier transform FFT

LED music spectrum display of the core algorithm is fast Fourier transform, FFT understanding and programming is more difficult, specially write this article to share the research results.

A thorough understanding of the Fourier transform

Rapid Fourier transform (Fast Fourier Transform) is a fast algorithm for discrete Fourier transform (FFT), which can transform a signal from time domain to frequency domain by FFT.

The process of analog signals being converted into digital signals via A/D conversion is called sampling. To ensure that the spectral shape of the signal after sampling is not distorted, the sampling frequency must be greater than twice times the maximum frequency component of the signal, which is called the sampling theorem.

Assuming the sampling frequency is FS, the sampling point is N, then the FFT result is a complex number of n points, each one corresponding to a frequency point, a point n (n starting from 1) The frequency expressed: fn= (n-1) *fs/n.

Example: 1kHz sampling frequency sampling 128 points, the FFT results of 128 data is the corresponding frequency point is 0,1k/128,2k/128,3k/128,...,127k/128 Hz.

The amplitude of this frequency point is: The modulus of the complex number divided by N/2 (N=1 is a DC component, whose amplitude is the modulus of the point divided by N).

The C language programming of Fourier transform

1, for fast Fourier transform FFT, the first problem to be solved is the code bit in reverse order.

Assuming an n-point input sequence, its ordinal binary digits are t=log2n.

Code bit reverse order to solve two problems: ① the T-bit binary number in reverse, ② the reverse of the two storage units to exchange.

① the t=3 bit binary number in reverse

② swap two units in reverse

If the input sequence of the natural order number I with the binary number, for example, if the maximum number is 15, that is, with 4 bits can represent n3n2n1n0, then the reverse after the J corresponding binary number is n0n1n2n3, then how to achieve reverse. Use the shift function of C language.

The procedure is as follows, I do not say much, can not understand the IQ must be below 180.

Complex number type definition and its operations

#define N 64//64 point

#define LOG2N 6//log2n=6

/* Plural type * *

typedef struct

{

float Real;

float img;

}complex;

Complex XData x[n]; Input sequence

/* plural addition * *

Complex Add (complex A,complex B)

{

Complex C;

C.real=a.real+b.real;

c.img=a.img+b.img;

return C;

}

/* Plural subtraction * *

Complex Sub (complex A,complex b)

{

Complex C;

C.real=a.real-b.real;

c.img=a.img-b.img;

return C;

}

/* plural multiplication * *

Complex Mul (complex A,complex B)

{

Complex C;

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

C.img=a.real*b.img + a.img*b.real;

return C;

}

/*** code bit reverse function ***/

void Reverse (void)

{

unsigned int i,j,k;

unsigned int t;

Complex temp;//Temporary Exchange variable

for (i=0;i<n;i++)//from No. 0 ordinal to N-1 ordinal

{

k=i;//Current First ordinal number

j=0;//stores the ordinal number after the reverse sequence, first initialized to 0

for (t=0;t<log2n;t++)//Shift T times, where log2n is a prior macro definition.

{

j<<=1;

J|= (k&1);//j left one and then add K to the lowest bit

K>>=1;//k to the right one, the secondary low to the lowest position

}

if (j>i)//If the reverse is greater than the original ordinal, the two storage units are exchanged (J>i is to prevent duplicate exchange)

{

Temp=x;

X=X[J];

X[j]=temp;

}

}

}

2, the second problem to be solved is the butterfly-shaped operation

① level 1th (1th column) each butterfly-shaped two-node "distance" is 1, the 2nd stage each Butterfly shape two node "distance" is 2, the 3rd level each butterfly shape two node "distance" is 4, the 4th level each Butterfly shape two node "distance" is 8. This pushes it,

The first M-class butterfly operation, each butterfly-shaped two-node "distance" l=2m-1.

② for the 16-point FFT, the 1th level has 16 groups of butterflies, each group has 1 butterflies, the 2nd level has 4 groups of butterflies, each group has 2 butterflies, the 3rd level has 2 groups of butterflies, each group has 4 butterflies, and the 4th level has 1 butterfly shapes, each group has 8 butterfly shape. This can be introduced,

For the N-point FFT, the M-class has n/2l group butterfly, each group has l=2m-1 butterfly shape.

Determination of ③ rotation factor

Taking 16-point FFT as an example, the first m-level K-rotation factor is, wherein k=0~2m-1-1, that is, the M-level common 2m-1 rotation factor, according to the rotation factor's irreducible sex, therefore the first M class K the rotation factor is, among them k=0~2m-1-1.

In order to improve the operation speed of FFT, we can set up a rotational factor array in advance and then realize it by look-up table method.

Complex code wn[n]=//rotation factor array
{//To save CPU calculation time, rotation factor using look-up table processing
★ According to the actual FFT point N, the table data needs to be modified by itself
The following results are automatically generated by Excel
Wn[k].real=cos (2*PI/N*K);
Wn[k].img=-sin (2*PI/N*K);

{1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
{0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
{0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
{0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
{0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
{ -0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
{ -0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
{ -0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
{ -1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
{ -0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
{ -0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
{ -0.38268,0.92388},{-0