FFT algorithm analysis and design based on embedded system

First, the derivation process of the real number FFT algorithm is analyzed, and then a C language program that implements the FFT algorithm is given, which can be directly applied to embedded systems such as microcontrollers or DSPs that require FFT operations.

The FFT algorithm of decimation in time (DIT) usually stores the original data in reverse order and finally outputs the results X(0), X(1), ..., X(k), ... in normal order. Assume that at the beginning, the data is in the array float dataR[128]. We represent the subscript i as (b6b5b4b3b2b1b0)b. The reverse order storage is to store the original element at the i-th position at the (b0b1b2b3b4b5b6)b position. Since the C language has a strong bit operation capability, we can extract b6, b5, b4, b3, b2, b1, b0 respectively and then reassemble them into b0, b1, b2, b3, b4, b5, b6, which are the reverse order positions. The program segment is as follows (assuming 128-point FFT):
/* i is the original storage position, and the final invert_pos is the inverted storage position*/
int b0=b1=b2=b3=b4=b5=6=0;
b0=i&0x01; b1=(i/2)&0x01; b2=(i/4)&0x01;
b3=(i/8)&0x01; b4=(i/16)&0x01; b5=(i/32)&0x01;
b6=(i/64)&0x01; /*The above statements extract the 0 and 1 values ​​of each bit*/
invert_pos=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;

You can compare it with the bit reversal program in the textbook and you will find that this algorithm makes full use of the bit operation capabilities of the C language. It is very easy to understand and the bit operation speed is very fast.

2 Derivation of the real number butterfly operation algorithm

Let’s first look at the butterfly diagram shown in Figure 1.

It should be noted that: XR(K) and XR'(K) have the same storage location, so after (1) and (2), the value at that location has changed. X'(K) is needed to calculate X(K+B) below, so when programming, be sure to save XR'(K) and XI'(K) into the two temporary variables TR and TI.

Similarly: XR(K+B)+jXI(K+B)= XR'(K)+jXI'(K)- [ XR'(K+B)+jXI'(K+B)] *[ cos（2πP/N）-jsin（2πP/N）]Continue to decompose and get the following two formulas:
XR(K+B)= XR'(K)-XR'(K+B) cos（2πP/N）- XI'(K+B) sin (2πP/N) （3）
XI(K+B)= XI'(K)+ XR'(K+B) sin（2πP/N）- XI'(K+B) cos (2πP/N) （4）
Note:
① When programming, XR'(K) and XI'(K) in formulas (3) and (4) are replaced by TR and TI respectively.

② After formula (3), the value of XR(K+B) has changed, and the previous value at this position is used in formula (4), so the previous value XR'(K+B) must be saved before executing formula (3).

③ When programming, XR(K) and XR'(K), XI(K) and XI'(K) use the same variable.
Through the above analysis, we only need to convert equations (1), (2), (3), and (4) into C language statements. Pay attention to the intermediate storage of variables, as shown in the following program segment.

/* Butterfly operation program segment, dataR[] stores the real part, dataI[] stores the imaginary part*/
/* Make a table of cos and sin functions, and directly look up the table to speed up the operation*/
TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];/*Save variables for use in the following statements*/
dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];

3 Analysis of the basic idea of ​​DIT FFT algorithm

We know that the N-point FFT operation can be divided into LOGN2 levels, each level has N/2 discs. The basic idea of ​​DIT FFT is to use 3 layers of loops to complete the entire operation (N-point FFT).

First-level loop: Since N=2m requires m levels of calculation, the first-level loop controls the number of operations.

Second-level loop: Since there are 2L-1 butterfly factors (multipliers) at the Lth level, the second-level loop is controlled according to the multipliers to ensure that the third-level loop is executed once for each butterfly factor. In this way, the third-level loop is controlled by the second-level loop, and each level performs 2L-1 loop calculations.

The third loop: Since there are N/2L groups in the Lth level, and the multipliers of different groups in the same level have the same distribution, after the second loop determines a multiplier, the third loop will calculate the butterfly with this multiplier once for each group in this level, that is, N/2L butterfly calculations will be performed each time the third loop is executed.

We can conclude that in each level, the third loop completes N/2L disk calculations; the second loop makes the third loop perform 2L-1 times, so when the second loop is completed, a total of 2L-1 *N/2L=N/2 disk calculations are performed. In essence, the second and third loops complete the calculation of level L.

A few data to note:

① In the Lth level, the two input ends of each disc are b=2L-1 points apart.

② The same multiplier corresponds to N/2L disks with adjacent points spaced 2L apart.

③ The P in the 2L-1 dishing factors WPN of the Lth level can be expressed as p = j*2m-L, where j = 0, 1, 2, ..., (2L-1-1).

The above is an analysis and study of the FFT algorithm in embedded systems. Readers can directly apply the algorithm to their own systems. Welcome to write to discuss. (Email: xiaowanang@163.net)

Attached is the 128-point DIT FFT function: