Fast Fourier transform method and device for improving utilization rate of storage space

Abstract

The invention relates to a fast Fourier transform method and device for improving the utilization rate of a storage space. The fast Fourier transform method comprises the following steps of: 1, writing input data into the storage space in an address close arrangement form; 2, calculating a twiddle factor, and storing the twiddle factor into a storage space; 3, FFT input data processing is carried out, specifically, step 31, FFT input data is subjected to data inverted-order storage; step 32, performing data operation, specifically, step 321, expanding the twiddle factors obtained in the step 2 to obtain a twiddle factor complete set; and step 322, performing butterfly operation, and obtaining output Fourier transform data after continuous traversal is completed. According to the method, the storage resources required by the input data and the output data are reduced by converting the floating point data format of the input data and adapting the adder and the multiplier of the computing resources; and the aim of saving memory resources is fulfilled on the premise of ensuring the calculation efficiency.

Description

Technical Field

[0001] The present invention belongs to the field of software technology, relates to digital signal processing technology, and particularly relates to a fast Fourier transform method and device for improving storage space utilization rate. Background Art

[0002] With the rapid development of electronic technology and integrated circuit technology, digital signal processing has been widely applied in fields such as automatic control. Digital signal processing is basically carried out through data in the time domain and frequency domain. In this process, the fast Fourier transform (FFT) of the discrete Fourier transform is one of the basic technologies in the field of digital signal processing.

[0003] In the field of speech recognition, the fast Fourier transform is used for spectral analysis of audio signals, which can convert time-domain signals into frequency-domain signals. It can also be used in multiple key steps such as signal feature extraction, noise suppression, model training, and real-time speech processing.

[0004] In the prior art, the CPU can be called by software to calculate the input data of the FFT. However, the input data of the FFT is large and the calculation process is complex, so it takes a long time. And the voice interaction technology has high requirements for real-time performance, so the method of implementing FFT by the CPU takes a long time.

[0005] In the prior art, there are also dedicated FFT processors that accelerate the fast Fourier transform process through hardware implementation. The most common method is to store data in 32-bit floating-point format. For each complex number, its real part and imaginary part need to be stored. For example, in the FFT algorithm with 512 input points, storing 512 input data requires 4KB of storage space, storing the results requires 4KB of storage space, and storing the twiddle factors also requires 4KB of storage space, totaling 12KB of storage space, which occupies a large amount of storage resources. In order to save the storage space of the twiddle factors, some hardware implementations also use SIN (sine) and COS (cosine) function calculation units to implement real-time calculation of the twiddle factors. However, the SIN and COS calculation units themselves have a large area and complex logic, a long combinational logic chain, and poor timing, resulting in low calculation efficiency.

[0006] In summary, in the implementation of the common FFT algorithm, it is difficult to solve the conflict between the use of memory resources and the area efficiency. Summary of the Invention

[0007] To overcome the technical defects existing in the prior art, the present invention discloses a fast Fourier transform method and device for improving storage space utilization rate.

[0008] The fast Fourier transform method for improving storage space utilization rate according to the present invention includes the following steps: Step 1. Write the FFT input data into the storage space in the form of address-packed arrangement; Step 2. Calculate the rotation factor W(j) = (-sin(2π*j / m), cos(2π*j / m)), where the rotation factor serial number j = (0, 1, 2…m / 8), and m is the FFT input point number representing the number of FFT input data; The rotation factor obtained in this step is called the basic rotation factor, and the basic rotation factor is stored in the storage space in the same format as the FFT input data in Step 1; Step 3. Perform FFT input data processing, specifically including: Step 31. For the FFT input data, perform data bit-reversal storage; Step 32. Perform data operations, specifically: Step 321. Expand all the basic rotation factors with the rotation factor serial number j = (0, 1, 2…m / 8) obtained in Step 2 to obtain the complete set of rotation factors; The expansion method is to set three interval separation points as T1 = m / 4, T2 = 3m / 8, and T3 = m / 2 respectively; When the interval of the rotation factor serial number j is [(m / 8)+1, T1], W(j) = [-W(T1 - j) r , -W(T1 - j) i , The previous number in the square brackets on the right side of the above equation represents the real part of W(j), and the latter number represents the imaginary part of W(j); When the interval of the rotation factor serial number j is [T1 + 1, T2], W(j) = [W(j - T1) r , -W(j - T1) i ; When the interval of the rotation factor serial number j is [T2 + 1, T3 - 1], W(j) = [W(T3 - j) i , -W(T3 - j) r ; W(T1 - j) i , W(T1 - j) r respectively represent the real part value and the imaginary part value of the rotation factor with the rotation factor serial number T1 - j, and so on; Step 322. Perform butterfly operations, and the butterfly operations need to be traversed log2(m) times; First traversal: Divide the input data into m / 2 pairs in ascending order of the storage address number. Each pair consists of 2 adjacent points in address. Input the data X0 and X1 of each pair and the rotation factor W(0) of the first traversal into the FFT complex operation module, calculate using formula 1, and store the obtained calculation results Y0 and Y1 into the addresses corresponding to the data X0 and X1 of this pair respectively as the input data for the next traversal. After calculating all pairs, enter the next traversal. Formula 1 is: Y0 = X0 + X1 * W Y1 = X0 - X1 * W; In formula 1, X0 and X1 are the input data of formula 1 respectively, Y0 and Y1 are the output data of formula 1 respectively, and W is the rotation factor used in formula 1. For the first traversal, W = W(0). Starting from the second time, the kth traversal is specifically as follows: Divide the m data obtained from the previous traversal into m / 2 k groups in ascending order of the storage address number. Each group consists of 2 k adjacent points in address. For the 2 k data in the same group, divide them into 2 k-1 pairs. Each pair has two data. The data within each pair are X n and X n+kn , where kn = 2 k-1 , and n is the number of data in the group in ascending order of the storage address number, n = 0, 1... 2 k-1 -1; For the data pair X n , X n+kn , use the rotation factor W(n * m / 2 k ), input into formula 1, calculate the output data Y0 and Y1, and store them into the addresses corresponding to X n , X n+kn respectively; Continue traversing until it ends after traversing to the log2(m)th time, and obtain the output Fourier transform data.

[0009] Preferably, in steps 1 and 2, the FFT input data and the basic rotation factor are both stored in the form of word granularity, and the word granularity includes 16-bit complex imaginary part data and 16-bit complex real part data.

[0010] Preferably, step 31 is specifically as follows: Set the first register and the second register, traverse all input data in ascending order of the storage address number, reverse the address number of each input data, and then judge whether the reversed address number is greater than the address number before reversal; If so, perform a reverse operation, specifically read the data at the current address sequence number addr1 and the address sequence number addr2 after reverse, store them in the first register and the second register respectively, then store the data in the first register at the address sequence number addr2 after reverse, and store the data in the second register at the current address sequence number addr1; Otherwise, skip the current address and continue to perform the above judgment on the address sequence number of the next input data until all input data are traversed.

[0011] The present invention discloses a fast Fourier transform device for improving the storage space utilization rate, which includes an FFT control unit module, a memory connected thereto, and an FFT complex operation module. The memory includes an input data storage memory and a rotation factor storage memory. The FFT complex operation module is used to perform complex operations in the fast Fourier transform method, and the FFT control unit module is used to implement the fast Fourier transform method.

[0012] Preferably, the memory further includes a first register and a second register.

[0013] The present invention reduces the storage resources required for input data and output data by converting the floating-point data format of the input data and adapting the adders and multipliers of the computing resources; and achieves the purpose of saving memory resources on the premise of ensuring the computing efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 is a schematic flowchart of a specific implementation manner of the fast Fourier transform method for improving the storage space utilization rate of the present invention; Figure 2 is a schematic diagram of a specific distribution manner of the rotation factor in the complex coordinate space of the present invention; Figure 3 is a schematic diagram of a specific implementation manner of the reverse operation of the present invention; Figure 4 is a schematic diagram of a specific implementation manner of the fast Fourier transform device for improving the storage space utilization rate of the present invention; Figure 5 is a schematic diagram of the input-output data relationship of Formula 1 of the present invention; Figure 6 is a schematic diagram of a specific implementation manner of the butterfly operation of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0015] The method for generating a clock reset network based on multi-constraint drive and self-adaptation described in the present invention, as Figure 1 shown, is implemented according to the following steps: The following combines the drawings to further elaborate on the specific implementation manners of the present invention: Before performing the transformation method of the present invention, divide the size of the storage resources required for the present invention. Define m as the number of FFT input points, then the specific storage space required is as follows: Data storage size = m * 2 * 2Byte Twiddle factor storage size = (m / 8 + 1) * 2 * 2Byte For example, when the FFT input point number m = 512, the memory size required for data storage is 2048Byte, and the memory size required for twiddle factor storage is 256Byte. Byte represents byte, and each byte is equal to 8bit. The FFT input point number m is the number of FFT input data.

[0016] Step 1. Import the FFT input data into the storage space. Specifically: Write 512 FFT input data into the storage space in an address-packed form at the word granularity. Each word granularity contains 16bit complex imaginary part data and 16bit complex real part data. The address-packed form means that the addresses are continuous without interruption.

[0017] Step 2. Calculate the twiddle factor W(j) = (-sin(2π*j / m), cos(2π*j / m)), where the twiddle factor serial number j = (0, 1, 2…m / 8), When the FFT input point number m = 512, it can be known that there are 65 twiddle factor serial numbers j, and then 65 complex value form twiddle factors W(j) can be obtained; In this step, the reason why the number of twiddle factors only takes m / 8 is that all twiddle factors are symmetrically distributed on the same circumference in the complex two-dimensional coordinate space. Only 1 / 8 of the twiddle factors need to be calculated. According to symmetry, the real part and imaginary part values of all the remaining twiddle factors can be obtained. For example Figure 2 As shown, when m = 16, each twiddle factor is centrosymmetrically distributed around the origin on the circumference. In this step, only m / 8 twiddle factors need to be temporarily stored. When performing operations later, all twiddle factors can be extended.

[0018] The twiddle factors obtained in Step 2 are called basic twiddle factors. Convert the imaginary part and real part of all basic twiddle factors into 16bit floating-point format data, that is, write them into the storage space at the word granularity for storage, and obtain a data form with the same format as the FFT input data in Step 1.

[0019] Step 3. FFT input data input Send a start flag signal to the FFT control unit module to perform FFT algorithm calculation. Step 31. Store the FFT input data in reverse order; The fast Fourier transform requires reversing the data. By reversing, it means swapping the stored data corresponding to the current address sequence number and the address sequence number after reversing the current address sequence number.

[0020] Figure 3 A specific implementation of the above reversal is given. Figure 3 In the upper part, the 16 input data from 0 to 15, with addresses from 0000 to 1111, are reversed. For example, Figure 3 In the second number, the address sequence number is 0001, and the address sequence number after reversing the current address sequence number is 1000. Then, swap the data 1 and 8 corresponding to the two address sequence numbers 0001 and 1000. After all reversals are completed, Figure 3 The reversed data in the lower part is obtained.

[0021] In the present invention, the following method can be used to reverse the FFT input data. Specifically: Set the first register and the second register. Traverse all input data in ascending order of the address sequence number where the data is stored. For each input data, judge the address sequence number. If the address sequence number after reversal is greater than the address sequence number before reversal, it means that the reversal has not occurred yet, and then perform the reversal operation. Specifically, read the data of the current address sequence number addr1 and the address sequence number addr2 after reversal, and store them in the first register and the second register respectively. Then, store the data in the first register into the address sequence number addr2 after reversal. e， Store the data in the second register into the current address sequence number addr1; Otherwise, skip the current address and continue to perform the above judgment on the address sequence number of the next input data until all input data are traversed; If it is equal to the sequence number after reversal, it means that the current position is the same as the position after reversal, and no reversal is needed; if the sequence number is less than the sequence number after reversal, it means that these two points have been swapped before, and no reversal is needed either.

[0022] The specific principle implemented by Step 31 is: Since the address sequence numbers are traversed in ascending order, except for individual address sequence numbers that do not change before and after reversal, such as special address sequence numbers like 0000, 0110, 1111, 1001, etc., most of the address sequence numbers in the first half of the data sequence are less than the address sequence numbers after reversal, so the reversal operation is performed; when traversing to the second half, since most of the address sequence numbers in the second half, except for the special address sequence numbers, are greater than the address sequence numbers before reversal, and the operations in the first half have already achieved reversal, no reversal is needed for the second half.

[0023] Step 32. Perform data operations. Specifically: Step 321. Expand all the basic rotation factors with rotation factor serial number j = (0, 1, 2…m / 8) obtained in Step 2 to obtain the complete set of rotation factors.

[0024] In the complete set of rotation factors, the rotation factor serial number j = (0, 1, 2…m / 2), that is, it expands from the basic rotation factors which are roughly only one-eighth of the FFT input point number m to the complete set of rotation factors including m / 2 rotation factors; The expansion method is to set the interval separation points as T1 = m / 4, T2 = 3m / 8, and T3 = m / 2 respectively; When the interval of the rotation factor serial number j is [(m / 8)+1, T1], W(j)=[-W(T1-j) r , -W(T1-j) i , It means that the first number in the square brackets of W(j) represents the real part of W(j), and the second number represents the imaginary part of W(j); When the interval of the rotation factor serial number j is [T1+1, T2], W(j)= [W(j-T1) r ,-W(j-T1) i ; When the interval of the rotation factor serial number j is [T2+1, T3-1], W(j)= [W(T3-j) i ,-W(T3- j) r ; W(T1-j) i , W(T1-j) r respectively represent the real part value and the imaginary part value of the rotation factor with rotation factor serial number T1-j, and so on for the rest.

[0025] After the expansion, m / 2 rotation factors are obtained, and the serial numbers are from 0 to (m / 2)-1 respectively.

[0026] The above expansion principle is based on the symmetric distribution of the rotation factors. As Figure 2 shown, taking m = 16 as an example, the rotation factors are symmetrically distributed on the circumference. When calculating three basic rotation factors with j = 0, 1, 2 in the fourth quadrant through Step 2, for j = 3, from Figure 2It can be seen that at this time, the rotation factors W(3) and W(1) are centrosymmetric in the fourth quadrant, that is, the imaginary part of the rotation factor with j = 3 is equal to the real part of the rotation factor with j = 1, and the real part of the rotation factor with j = 3 is equal to the imaginary part of the rotation factor with j = 1; through similar analysis, it can be obtained that the imaginary part of the rotation factor with j = 4 is equal to the real part of the rotation factor with j = 0, and the real part of the rotation factor with j = 4 is equal to the imaginary part of the rotation factor with j = 0; similarly, for the other rotation factors distributed in the other three quadrants, similar conclusions can be drawn, that is, the real and imaginary parts of each of the remaining rotation factors can be represented by exchanging the positions and adding or not adding a negative sign to the real and imaginary parts of the three basic rotation factors, and finally all 16 rotation factors are obtained. The reason for only calculating m / 2 rotation factors in step 321 is that in the subsequent butterfly operations, since the two rotation factors in each butterfly operation are symmetric about the origin, only the first m / 2 rotation factors are needed.

[0027] Step 322. Perform butterfly operations, which need to be traversed log2(m) times; for example, when m = 512, the butterfly operations need to be traversed 9 times; the butterfly operation is the minimum unit calculation method of FFT, which is used to combine the calculation results of two discrete Fourier transforms (DFT). Each butterfly operation is a complex number calculation, using the rotation factor to independently process two input data, which is suitable for hardware parallel acceleration.

[0028] The first traversal: Divide the input data into m / 2 pairs in ascending order of the storage address number. Each pair consists of 2 adjacent points in address. Input the data X0 and X1 of each pair and the rotation factor W(0) of the first traversal into the FFT complex operation module, and calculate using formula 1. Store the obtained calculation results Y0 and Y1 into the addresses corresponding to the data X0 and X1 of this pair respectively as the input data for the next traversal. After calculating all pairs, enter the next traversal.

[0029] Formula 1 is Y0 = X0 + X1 * W Y1 = X0 - X1 * W; The input and output data processing methods in formula 1 are as Figure 5 shown. X0 and X1 are the input data of formula 1 respectively, Y0 and Y1 are the output data of formula 1 respectively, and W is the rotation factor used in formula 1. For the first traversal, W = W(0); each parameter such as X0, X1, and W is in complex number form, and the calculation is a complex number operation method. Starting from the second time, the kth traversal is specifically: Divide the m data obtained from the previous traversal into m / 2 k groups in ascending order of the storage address number. Each group consists of 2 k adjacent ones in address. For the 2 in the same groupk A data is divided into 2 k-1 pairs, and the data within each pair is X n and X n+kn , where kn = 2 k-1 , n is the number assigned to the data within the group in ascending order of the storage address sequence number, n = 0, 1…2 k-1 -1; For the data pair X n , X n+kn , the rotation factor W(n * m / 2 k ) is adopted, and substituting into Formula 1, the output data Y0 and Y1 are calculated and respectively stored into the addresses corresponding to X n , X n+kn ; As Figure 6 shown, taking m = 512, and k = 2 and 3 for the second and third traversals as examples: In the second traversal, the input data is divided into 512 / 4 = 128 groups, and the 4 data in each group are X0, X1, X2, X3 in ascending order of the storage address sequence number. They are divided into two pairs. The first pair is X0, X2; the second pair is X1, X3; for the first pair, n = 0, and the rotation factor W(0) is adopted. For the second pair, n = 1, and the rotation factor W(1 * 512 / 4) = W(128) is adopted, where W(128) can be obtained after the rotation factor expansion according to Step 321.

[0030] Substituting into Formula 1 for calculation, for the first pair, Y0 = X0 + X2 * W(0) Y1 = X0 - X2 * W(0); For the second pair, Y0 = X1 + X3 * W(128) Y1 = X1 - X3 * W(128).

[0031] Four output values are calculated and respectively stored into the addresses corresponding to X0, X1, X2, X3.

[0032] When k = 3 in the third traversal; In the third traversal, the input data is divided into 512 / 8 = 64 groups, and the 8 data in each group are X0, X1, X2…X7 in ascending order of the storage address sequence number. They are divided into four pairs. The first pair is X0, X4; the second pair is X1, X5; the third pair is X2, X6; the fourth pair is X3, X7; For the first pair, n = 0, the rotation factor W(0) is adopted; for the second pair, n = 1, the rotation factor W(1 * 512 / 8) = W(64) is adopted; for the third pair, n = 2, the rotation factor W(2 * 512 / 8) = W(128) is adopted; for the fourth pair, n = 3, the rotation factor W(3 * 512 / 4) = W(192) is adopted, where W(128) and W(192) can be obtained after the rotation factor expansion according to step 321.

[0033] Perform continuous traversal until it ends after traversing log2(m) times, and the output Fourier transform data is obtained.

[0034] The fast Fourier transform method described in the present invention can be implemented based on the fast Fourier transform device as Figure 4 shown, which includes an FFT control unit module, a memory connected thereto, and an FFT complex arithmetic module. The memory includes an input data storage memory and a rotation factor storage memory. The FFT complex arithmetic module is used to perform complex arithmetic in the fast Fourier transform method, and the FFT control unit module is used to implement the fast Fourier transform method. Figure 4 In the specific embodiment shown, the memory further includes a first register and a second register for implementing the bit-reversal operation.

[0035] The foregoing are the preferred embodiments of the present invention. If the preferred embodiments in each preferred embodiment are not obviously self-contradictory or premised on a certain preferred embodiment, each preferred embodiment can be arbitrarily superimposed and combined for use. The embodiments and the specific parameters in the embodiments are only for clearly expressing the inventor's invention verification process and are not used to limit the patent protection scope of the present invention. The patent protection scope of the present invention still depends on its claims. All equivalent structural changes made by using the content of the specification and drawings of the present invention should, by the same token, be included in the protection scope of the present invention.

Claims

1. A fast Fourier transform method for improving the utilization rate of storage space, characterized in that, It includes the following steps: Step 1. Write the FFT input data into the storage space in the form of address-packed arrangement; Step 2. Calculate the rotation factor \(W(j)=(-\sin(2\pi*j / m),\cos(2\pi*j / m))\), where the rotation factor serial number \(j = (0,1,2\cdots m / 8)\), and \(m\) is the FFT input point number representing the number of FFT input data; The rotation factor obtained in this step is called the basic rotation factor, and the basic rotation factor is stored into the storage space in the same format as the FFT input data in Step 1; Step 3. Perform FFT input data processing, specifically including: Step 31. For the FFT input data, perform data reverse storage; Step 32. Perform data operations, specifically: Step 321. Expand all the basic rotation factors with the rotation factor serial number \(j=(0,1,2\cdots m / 8)\) obtained in Step 2 to obtain the complete set of rotation factors; The expansion method is to set three interval separation points as \(T1 = m / 4\), \(T2 = 3m / 8\), \(T3 = m / 2\) respectively; When the interval of the rotation factor serial number \(j\) is \([(m / 8)+1,T1]\), W(j)=[-W(T1 - j) r , -W(T1 - j) i ​ The number in the front of the square brackets on the right side of the above formula represents the real part of \(W(j)\), and the number in the back represents the imaginary part of \(W(j)\); When the interval of the rotation factor serial number \(j\) is \([T1 + 1,T2]\), W(j)= [W(j - T1) r , -W(j - T1) i ; When the interval of the rotation factor serial number \(j\) is \([T2 + 1,T3 - 1]\), W(j)= [W(T3 - j) i , -W(T3 - j) r ; W(T1 - j) i , W(T1 - j) r respectively represent the real part value and the imaginary part value of the rotation factor with the rotation factor serial number T1 - j, and so on for the rest; Step 322. Perform butterfly operations, and the butterfly operations need to be traversed \(\log_2(m)\) times; The first traversal: Divide the input data into \(m / 2\) pairs evenly according to the storage address serial number from low to high. Each pair is 2 points with adjacent addresses. Input the data \(X0\) and \(X1\) of each pair and the rotation factor \(W(0)\) of the first traversal into the FFT complex operation module, calculate using Formula 1, and store the obtained calculation results \(Y0\) and \(Y1\) into the addresses corresponding to the data \(X0\) and \(X1\) of this pair respectively as the input data for the next traversal; After calculating all pairs, enter the next traversal; Formula 1 is: \(Y0 = X0+X1*W\) \(Y1 = X0 - X1*W\); In Formula 1, \(X0\) and \(X1\) are the input data of Formula 1 respectively, \(Y0\) and \(Y1\) are the output data of Formula 1 respectively, and \(W\) is the rotation factor adopted by Formula 1. For the first traversal, \(W = W(0)\); Starting from the second time, the \(k\)th traversal is specifically: Divide the m data obtained in the previous traversal into m / 2 groups in ascending order of storage address numbers. k Each group consists of 2 k adjacent addresses. For the 2 k data in the same group, divide them into 2 k-1 pairs. Each pair has two data, and the data within each pair are X n and X n+kn , where kn = 2 k-1 , n is the number of data in the group in ascending order of storage address number, n = 0, 1... 2 k-1 - 1; For data pair X n , X n+kn , the rotation factor W (n*m / 2 k ) is adopted, and formula 1 is input to calculate the output data Y0 and Y1, which are respectively stored in the addresses corresponding to X n , X n+kn ; Perform continuous traversals until the \(\log_2(m)\)th traversal is completed and then end to obtain the output Fourier transform data.

2. The fast Fourier transform method according to claim 1, wherein In Steps 1 and 2, both the FFT input data and the basic rotation factor are stored in the form of word granularity, and the word granularity includes 16-bit complex imaginary part data and 16-bit complex real part data.

3. The fast Fourier transform method according to claim 1, characterized in that, Step 31 is specifically to set the first register and the second register, traverse all input data in ascending order according to the data storage address serial number, reverse the address serial number of each input data, and then judge whether it is greater than the address serial number before the reverse; If so, perform a reverse operation. Specifically, read the data of the current address number addr1 and the address number addr2 after reverse, and store them in the first register and the second register respectively. Then store the data in the first register into the address number addr2 after reverse e， Store the data in the second register into the current address number addr1; Otherwise, skip the current address and continue to perform the above judgment on the address serial number of the next input data until all input data are traversed.

4. A fast Fourier transform device for improving the storage space utilization rate, characterized in that, It includes an FFT control unit module, a memory connected thereto, and an FFT complex arithmetic module. The memory includes an input data storage memory and a twiddle factor storage memory. The FFT complex arithmetic module is used to perform complex arithmetic in the fast Fourier transform method described in claim 1, and the FFT control unit module is used to implement the fast Fourier transform method described in claim 1.

5. The fast Fourier transform device according to claim 4, wherein The memory further includes a first register and a second register.

Images