Digital signal Processing--fft and butterfly-shaped algorithm

In digital signal processing, it is often necessary to use discrete Fourier transform (DFT) to obtain the frequency domain characteristics of the Signal. Although the traditional DFT algorithm can obtain the signal frequency domain characteristics, The algorithm is computationally large and time consuming, which is not conducive to the real-time processing of the Signal. therefore, since the DfT has been discovered, it can not be applied to the actual project for a long time, until a fast discrete Fourier calculation method is--fft, It is found that the discrete Fourier transform is widely used in the actual engineering. It should be emphasized that FFT is not a new method of frequency domain feature acquisition, but a fast algorithm for DFT. In this paper, the principle of FFT and the specific implementation process are explained in Detail.

DFT Calculation formula

In this paper, the calculation formula of DFT is given directly without deduction:

where x (n) represents the input discrete digital signal sequence, WN is the rotation factor, and x (k) is the relative amplitude of n discrete frequency points corresponding to the input sequence x (n). In general, assuming that x (n) is from Low-pass sampling and the sampling frequency is fs, then x (k) represents the relative amplitude of THE-FS/2 rate starting at the frequency interval of fs/n to the N frequency points fs/2-fs/n up to. Because a set of discrete frequency amplitude values obtained by DFT is actually changed from cycle to frequency axis, i.e. X (k+n) =x (k). therefore, any successive N points can be expressed as a DFT calculation effect, the negative frequency composition is abstract, difficult to understand, according to the X (k) periodic characteristics, So we can think that x (k) is the beginning of the zero frequency, frequency interval of fs/n, to fs-fs/n up to the N frequency point relative amplitude.

Calculation of N-point DFT

According to the DfT formula given in (1), we can know that each calculation of a frequency point x (k) requires n-times complex multiplication and N-1 complex addition, and the x (k) of the n points requires n^2 complex multiplication and n (N-1) complex Addition. When x (n) is a real number, the DfT that computes N points requires a 2*n^2 real-time multiplication, 2*n* (N-1) second real number Addition.

Characteristics of the rotational factor WN

Symmetry of the 1.WN

Periodicity of 2.WN

3.WN of availability

Based on these properties, we can get a series of useful results from formula (5)

Derivation of base-2 FFT algorithm

Suppose that the sample sequence number is n=2^l,l as an integer, and if this condition is not met, you can artificially add several 0 to make the sampling sequence points meet this Requirement. first, we divide the sequence x (n) into two groups according to the odd couple, as Follows:

So according to the DfT formula (1) Is:

thus, we convert an N-point DFT into the form of the formula (7), at which time the value of K is 0 to N-1, now divided into two sections to discuss, when K is 0~n/2-1, because X1 (r), x2 (r) is the sequence of N/2 points, therefore, this time formula (7) can be written as:

When k value is n/2~n-1, K is substituted with K ' +n/2, K ' value is 0~n/2-1. The formula (7) simplification can be:

Based on the above derivation, we can get the following conclusion: the DfT transformation process of n-point can be represented by the DfT transformation process of two N/2 points, and its concrete formula (10) shows the iterative formula of DFT fast Algorithm:

In the above formula, X ' (k ') is a discrete Fourier transform of even-numbered branches, and X ' (k ') is a discrete Fourier transform of the odd Branch. The calculation process of formula (10) can be visualized by using the butterfly algorithm flow graph of Figure 1.

Fig. 1 The time extraction method of butterfly operation flow diagram

In Figure 1, the DFT output of two N/2 points is the DFT result of an n point, and the input and output points are the Same. Using this notation, the 8-point DFT can be represented by the Diagram:

Figure 2 4-point decomposition of 8-point DFT

According to the formula (10), the DfT of an n-point can be composed of two n/2-point dft, and the first decomposition of the 8-point DFT of Figure 2 can be obtained by combining the butterfly signal flow diagram of Figure 1. The decomposition can be described in the following steps:

1. Divide the input sequence of n points into 2 groups of sequences of N/2 points, respectively, by parity

2. DfT transforms each set of sequences in 1 to obtain the DfT transformation values of two sets of points n/2 X1 and X2

3. Combining 2 results into one n-point DFT transformation results according to the butterfly signal flow graph

According to Formula (10) we can further decompose the 4-point DFT in Figure 2, and get the result of Figure 3, the decomposition step is consistent with the Previous.

Figure 3 2-point decomposition of 8-point DFT

finally, the 2-point DfT is further decomposed to obtain the final 8-point FFT signal calculation flow graph:

Figure 4 Total decomposition of the 8-point DFT

From the process of Fig. 2 to figure 4, It is necessary to explain the changing law of rotation coefficient. Seems to push forward a level, in the odd part of the rotation factor factor increment seems to be larger, in fact, it is Not. In fact, the rotation factor index of the odd grouping part is fixed at 1 per increment, only because each forward time, the number of data in the grouping sequence is less, in order to unify the use of the original data n-based rotation factor is caused by the Transformation. The coefficient of the odd part of each packet is wn, where n is the number of sequential points before this packet. The above side of the 8-point DFT for example, The first packet n=8, the second group N is 4, in order to unify the formula (4) for the transformation of N into 8, but the corresponding need for the exponent multiplied by 2.

Computational capacity of N-point base-2 FFT algorithm

From Figure 4, we can see that the FFT transform of N-point DfT can be converted to log2 (n) cascade Butterfly operations, Each of which contains n/2-times butterfly Calculation. Each butterfly operation consists of 1 complex multiplication and 2 complex addition. So the total computational amount of the N-point FFT calculation is: complex multiplication--n/2xlog2 (n) complex addition--nxlog2 (n). Assuming that the sampled sequence is a sequence of real numbers, then only the first level is calculated as a mixed calculation of real and complex numbers, and after one iteration the subsequent computations are made into complex numbers, which are inconsistent with the direct DFT calculations. therefore, it has little effect on the efficiency of the FFT algorithm for the input sequence to be complex or real. A complex multiplication consists of 4 real-time multiplication, 2 real-time addition, and one-time complex addition containing 2 complex additions. therefore, for the FFT calculation of N points, The total number of real multiplication is: 2xnxlog2 (n); The total number of complex additions is: 2xnxlog2 (n).

Implementation method of N-point base-2 FFT algorithm

From Figure 4 We can summarize the flow of a fast DfT calculation method for points N=2^l:

1. Reverse-order transformation for input data series.

The purpose of this transformation is to enable the output to get the sequential sequence of X (0) ~x (N-1). Also in the case of the 8-point dft, the transformation converts sequential input sequences x (0) ~x (7) to 4 x (0), x (4), x (2), x (6), x (1), x (5), x (3), x (7) sequence. It is implemented by assuming sequential input sequence once the village is in an array element of a (0) ~a (N-1), we first make the array subscript binary (example: for a sequence with a number of points of 8 only LOG2 (8) = 3 bit binary sequence representation, The ordinal 6 is represented as 110). Binary after binary is the binary sequence is inverted, the process of the inverted bit is the original sequence from right to left to compose a new sequence, such as the number of 6 binary representation of 110, inverted to 011, even if the decimal 3. The third step is to swap the data corresponding to the ordinal of the inverted position before and after the Reversal. Still taking the sequence number 6 as an example, the interchange process is as Follows:

temp = A (6); A (6) = a (3); A (3) = a (6);

In fact, given the efficiency of execution, it is a waste of time to use this process for every data Input. We can use pointers to pointers to implement this procedure, or use a pointer array to implement the PROCESS.

2. The circular structure of the butterfly-shaped Operation.
From Figure 4 we can see that for the number of n = 2^l FFT operation, can be decomposed into the l-order butterfly graph cascade, Each order of the butterfly shape is divided into m butterfly group, each butterfly-shaped group contains k-shaped butterfly. Based on this, we can construct a Third-order loop to achieve the butterfly-shaped Operation. The programming process needs to be aware that the rotation factor is associated with the number of butterflies in the butterfly Order and in the butterfly Group.

3. Key issues to be aware of floating point-to-point conversions
The above analysis is based on floating point arithmetic to get the conclusion, in fact, most embedded systems support the floating point operation is very little, so in the embedded system, the discrete Fourier transform should be used in fixed-point mode. For simple DFT operations, It is very easy to get from floating point to fixed point. According to the formula (1), assuming that the input x (n) is an ad-sampled number sequence, the ad bit number is 12 bits, the input signal range is 0~4096. In order to perform the fixed-point operation, we will enlarge the imaginary part of the real part of the rotation factor 2^12 times, and take the integer part to represent the rotation factor. After that, we can calculate according to the formula (1), the result is proportional to the original result, The new result is 2^12 times than the original Result. however, for FFT using the butterfly operation, We cannot use this simple amplification rotation factor to convert to integer calculation. Because the FFT is an asymmetric iterative process, Suppose we enlarge the rotation factor, according to the butterfly flow graph we can find that the final result is that the different inputs are magnified by different multiples, for the first input x (0) will never Enlarge. To give a more vivid example, take Figure 4 for Example. You can see that the X (0) on the right can be expressed directly in the following Formula:

From the above we can see the rotation factor multiplied by the different input itemsnumber(note here is the number, even if it is wn^0, is also considered in, because wn^0 equals 1 when no magnification, all rotation factor index modulus is not 1, so need to consider). This results in an unbalanced input and incorrect Operation. After consulting the relevant data, it is more appropriate to first of all rotation factors are magnified 2^q times, Q must be greater than equal to l, in order to ensure the difference of different rotation factors. Rotation factor amplification, in order to ensure that the modulus of 1, in each of the product operation of the butterfly operation we need to shift the result to the right Q bit to counteract this amplification, so as to obtain the correct Result. The reason why the use of magnification must be the power of the whole number of 2 is why we can then cancel the previous amplification by simple right shift operation, and the right shift replaces the division operation, which saves time greatly.

4. Overflow issues during the calculation
One of the last issues to be aware of is the overflow problem in the computational process. In practical applications, Although the ad has a 12 bit width, the sampled signal may be small, for example, may fluctuate between 0~8, that is, it may actually only be 3 bits in the Case. In this case, in order to not lose information in the calculation process, it is generally necessary to first shift the input data to the left P-bit amplification, the data amplification may lead to overflow, so that the calculation error, and the limit of overflow is the case: assume that we have a data bit width of d-bit (not including the sign bit), AD sample bits B , the rotation is therefore magnified by the q-bit, and the Fft-level transport calculates the maximum number of double L Bits. We get:

Assuming the ad bit width 12, the data bit width 32, The sign bit 1 bits, so the effective bit width 31 bits, The sampling number n, then we can obtain LOG2 (N) +p+q<=19, assumes the point 128, and q>=l can obtain the magnification p<=5.

FFT code example
Based on the above analysis, I have compiled a 128-point FFT code as follows, the code in the Keil simulation run, no problem found.

#define N #define POWER 6//this value represents the number of times the input data is magnified first, the magnification is 2^power #define P_COEF 8//the value represents the number of times the rotation factor is magnified, and the magnification is 2^power #if (n = = 4) #define L 2//L is defined to satisfy L = log2 (n) #elif (n = = 8) #define L 3//L definition satisfies L = log2 (n) #eli F (n = =) #define L 4//L definition satisfies L = log2 (n) #elif (n = =) #define L 5//L definition satisfies L = log2 (n) #elif (n = = () #define L 6//L definition satisfies L = log2 (n) #elif (n = =) #define L 7//L is defined to meet L = log2 (n) #endif//ad sample number of bits is 1 2-bit, The data space can be defined using S16 x[n], but in order to conserve storage space, the FFT results will also be stored in that variable space. Because THE//FFT affected variables need to be redefined, the way they are defined and the specific reasons are as Follows: The//fft process is multiplied by a larger rotation factor, so a 32-bit definition//fft process produces negative numbers, so a symbolic definition is required//fft the complex number is present, so you need to define a two-dimensional array [][0]  Store the real part, [][1] store imaginary part S32 x[n][2] = {};  *p[n] is defined so that after the first pointer is initialized, the array pointer points to the data variable x//p[i][0] in reverse order, representing the real part pointing to the data, p[i][1] represents the imaginary part of the data S32 *p[n]; Defines the rotation factor matrix//rotation factor matrix that stores wn^0,wn^1,wn^2...wn^ (n/2-1), which n/2 a complex rotation factor s16 w[n>>1][2] = {256,0,256,-13,255,-25,253,- 38,251,-50,248,-62,245,-74,241,-86,237,-98,231,-109,226,
-121,220,-132,213,-142,206,-152,198,-162,190,-172,181,-181,172,-190,162,-198,152,
-206,142,-213,132,-220,121,-226,109,-231,98,-237,86,-241,74,-245,62,-248,50,-251,38,
-253,25,-255,13,-256,0,-256,-13,-256,-25,-255,-38,-253,-50,-251,-62,-248,-74,-245,-86,
-241,-98,-237,-109,-231,-121,-226,-132,-220,-142,-213,-152,-206,-162,-198,-172,-190,-181,
-181,-190,-172,-198,-162,-206,-152,-213,-142,-220,-132,-226,-121,-231,-109,-237,-98,-241,
                  -86,-245  , -74,-248,-62,-251,-50,-253,-38,-255,-25,-256,-13};    void Init_pointer (void) {unsigned char i = 0;    unsigned char j = 0;       unsigned char k = 0;       For (i = 0; i < N; i++) {j = 0;       For (k = 0; k < L; k++) {j |= ((((i >> k) &0x01) << (l-k-1));    } p[i] = &x[j][0]; }}///*description: the base 2fft main function, which uses the pointer array p to perform an FFT calculation of the global variable array x and stores the result in the array x *global var:rw->x; R->p,w.   (r means read, w means write, rw means read/write)/void fft2 (void) {u8 i;//i is used to denote the order of the Butterflies Cascade U8 j;//represents The sequence of the starting point of the butterfly grouping, The butterfly grouping span is 2^i U8 k;//k indicates the number of even branch sequences within the butterfly Group. U8 gp_distance = 1;//the Butterfly grouping span U8 temp;//temp is used to record half the distance between the current group, and also represents the distance between the same disc two branch U8 GP_HF = 0;//record The middle point ordinal of the current group U8 Delta = n;//    The rotation factor subscript increment, Originally the subscript initial value should be n/2, But.    S16 *PW = & (w[0][0]);       int tp_result[2];      Product//input signal sequence for temporary storage of rotation factor and odd packet amplification for (i = 0; i < N; i++) {x[i][0] <<= POWER; x[i][1] <<= POWER;      } for (i = 0; i < L; i++) {temp = gp_distance;      Gp_distance <<= 1;         For (j = 0; J < N; j+=gp_distance) {GP_HF = temp + j;         PW = & (w[0][0]);  For (k = j; K < gp_hf; K++)//complete a set of complex multiplication processes within a group of all butterfly operations {//butterfly Operations tp_result[0] = pw[0] * (p[k+temp][0])-            pw[1] * (p[k+temp][1]);            tp_result[1] = pw[0] * (p[k+temp][1]) + pw[1] * (p[k+temp][0]);            tp_result[0] >>= p_coef;            tp_result[1] >>= p_coef;            2 groups of complex addition processes in butterfly operations p[k+temp][0] = p[k][0]-tp_result[0];            p[k+temp][1] = p[k][1]-tp_result[1];            p[k][0] + = tp_result[0];                         p[k][1] + = tp_result[1];         PW + = delta;   }} Delta >>= 1; } }

Digital signal processing--fft and butterfly algorithm