CN101937424A - FPGA (Field Programmable Gate Array) based method for realizing high-speed FFT (Fast Fourier Transform) processing - Google Patents

Abstract

The invention discloses a FPGA (Field Programmable Gate Array) based method for realizing high-speed FFT (Fast Fourier Transform) processing. The method comprises the following steps of: combining a discrete fourier transform X(k) and an input signal x(n) into a discrete fourier transform pair, wherein the discrete fourier transform X(k) and the input signal x(n) are complex sequences with the same length of N, decomposing the discrete fourier transform X(k) and the input signal x(n) into short sequences, wherein the even number items are in one group and the odd number items are in another group, and thus, the subsequences of two N/2 points are obtained; and so forth, obtaining the rapid algorithm FFT of the discrete fourier transform. By improving the FFT algorithm structure, the invention obtains high-speed parallel FFT so that the FFT operation speed is greatly improved compared with the traditional method.

Description

Realize the method that high speed FFT handles based on FPGA

Technical field

The invention discloses a kind of method that fpga chip realizes that high speed FFT handles of using, it relates to digital processing field.

Background technology

DFT (Discrete Fourier Transfer) and fast algorithm FFT thereof are the core ingredients of digital processing field, and fft algorithm is varied, can be divided into base 2, base 4 etc. again by data pick-up mode difference.Each algorithm relative merits is looked different restraining factors and difference.The FFT implementation method is also varied, can realize with software, also can realize with hardware.Realize that with software then computing velocity is very slow on PC or workstation, the concrete system of general many combinations realizes with hardware, for example realize with single-chip microcomputer or DSP (DigitalSignal Process) processor, but speed still can not satisfy the demand of some application, is difficult to and A/D device matching fast.Pursue arithmetic speed in applications such as radar, electronic countermeasures, require very high real-time.FPGA is the abbreviation of English Field Programmable Gate Array, it is field programmable gate array, it is the product that further develops on the basis of PAL (Programmable Array Logic), GAL (GateArray Logic), EPLD programming devices such as (Embedded Programmable Logic Device), has higher integrated level, stronger logic realization ability and better design flexibility.FPGA is made up of many independently programmed logical module, the user can couple together these modules by programming and realize different designs, it is as special IC (ASIC, Application Specific Integrated Chip) a kind of semi-custom circuit in the field occurs, both solve the deficiency of custom circuit, overcome the limited shortcoming of original programming device gate circuit number again.In the last few years, the FPGA technical development was rapid.On the one hand, the programmable logic device (PLD) of various high capacity, high-performance, low-power consumption is constantly released.On the other hand, many FPGA design aids occurred, these instruments have improved the design efficiency based on the new type integrated circuit of FPGA greatly, make more low-cost, more short-period complex digital system development become possibility.

In the digital signal processing, establishing input signal X (n) is that length is the complex sequences of N, and its DFT (DiscreteFourier Transfer) is defined as:

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{\mathrm{nk}},(0\≤k\≤N-1)---\left(1\right)$

Its IDFT (Inverse Discrete Fourier Transfer) is defined as:

$x\left(n\right)=\mathrm{IDFT}\left[X\left(k\right)\right]=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{-\mathrm{nk}},(0\≤k\≤N-1)---\left(2\right)$

Wherein,

Formula (1) is called the discrete fourier direct transform, and formula (2) is called the discrete fourier inverse transformation, and it is right that x (n) and X (k) have constituted discrete Fourier transformation.

According to above-mentioned formula, calculate an X (k), need N complex multiplication and N-1 complex addition, and calculate whole X (k) (0≤k≤N-1), need N altogether ²The inferior complex addition of inferior complex multiplication and N (N-1), directly calculating whole X (k) needs 4N altogether ²The inferior real number addition of inferior real multiplications and 2N (2N-1).As seen workload and N ²Be directly proportional, N is big more, and operand will significantly increase.

List of references:

[1]Cooley?J?W，Tukey.An?algorithm?for?the?machine?calculationof?complex?Fourier?series[J].Math.Compute，1965，Vol.19，No.4，p297～301

[2] Ren Shuyan. the fft algorithm of using VHDL language is realized [J], Harbin University of Science and Technology's journal, Vol.8, No.6,2003.12, P.24-26

[3] Liu Chaohui, Han Yueqiu. realize the research [J] of FFT with FPGA. Beijing Institute of Technology's journal [J], Vol.19, No.2,1999.4, P.234-238

[4] Wang Xudong, Liu Yu. full parallel structure for FFT and FPGA thereof realize [J]. Nanjing Aero-Space University's journal, Vol.38 No.1 Feb.2006, p96-100

[5]Lo?Sing?Cheng，Ali?Miri，Tet?Hin?Yeap.Efficient?FPGAimplementation?of?FFT?based?Multipliers[C]，Proc.IEEE?CCECE/CCGEI，Saskatoon，May?2005，pp.1300-1303

[6]Ayan?Banerjee，Anindya?Sundar?Dhar，Swapna?Banerjee.FPGArealization?of?a?CORDIC?based?FFT?processor?for?biomedical?signalprocessing[J]，Microprocessor?and?Microsystems?2001，25(1)，pp.131-142

[7]PEREZ-PASCUAL，A.，SANSALONI，T.，and?VALLS，J.On-lineradix-2?butterflies?on?FPGA[C].Proc.WSES?International?Conferenceon?Speech，signal?and?image?processing，2002，pp.2201-2205

[8]PEREZ-PASCUAL，A.，SANSALONI，T.，and?VALLS，J.FPGA?basedradix-4?butterflies?for?HiperLAN2[C].Proc.IEEE?InternationalSymposium?on?Circuits?and?Systems(ISCAS?2002)，2002，pp.III.277-III.280

[9]Z.Szadkowski.16-point?Discrete?Fourier?transform?based?onthe?Radix-2?FFT?algorithm?implemented?into?Cyclone?FPGA?as?the?UHERCtrigger?for?horizontal?air?showers[C]，Proc.of?the?29th?ICRC，Pune，2005.

[10]Wen-Chang?Yeh，Chein-Wei?Jen.High-Speed?and?Low-PowerSplit-Radix?FFT[J].IEEE?Tans.Signal?Processing，2003，Vol.51，No.3，p864-874.

[11]Yutai?Ma，Lars?Wanhammar.A?Hardware?Efficient?Control?ofMemory?Addressing?for?High-Performance?FFT?Proces?sors[J].IEEE?Tans.Signal?Processing，2000，Vol.48，No.3，p917-920.

Summary of the invention

Technical matters: the present invention seeks to improve arithmetic speed into reducing operand, the defective that exists at background technology provides a kind of method that realizes high speed FFT processing based on FPGA.

Technical scheme: the present invention adopts following technical scheme for achieving the above object:

The present invention is based on the method that FPGA realizes that high speed FFT handles, it is right to it is characterized in that discrete Fourier transformation X (k) and input signal x (n) have constituted discrete Fourier transformation, all is that length is that the complex sequences of N, sequence length N are 2 integer power power, i.e. N=2 ^M, wherein M is a positive integer, expression power number of times, and wherein k, n represent the sequence of discrete Fourier transformation X (k) and input signal x (n) respectively;

At first discrete Fourier transformation x (n) is decomposed into two groups, even number Xiang Weiyi group, odd term is one group, obtains the subsequence that two N/2 are ordered, that is:

x ₁(r)＝x(2r)，x ₂(r)＝x(2r+1)，(0≤r≤N/2-1) (3)

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left(n\right)W_N^{\mathrm{kn}}=\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x\left(2r\right){\mathrm{<figure-callout\; id="2"\; label="W\; N"\; filenames="HSA00000266303200011.png,HSA00000266303200012.png"\; state="\{\{state\}\}">W</figure-callout>}}_N^{}+\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x(2r+1)W_N^{k(2r+1)}---\left(4\right)$

Utilize

Formula (4) can be write as:

$\left\{\begin{array}{c}X\left(k\right)=X_1\left(k\right)+W_N^k{X}_2\left(k\right)\\ X(k+N/2)=X_1\left(k\right)-W_N^k{X}_2\left(k\right)\end{array}\right.,(0\&le;k\&le;N/2-1)---\left(5\right)$

X wherein ₁(k) and X ₂(k) be respectively x ₁(r) and x ₂(r) DFT.

Beneficial effect:

The operation time of its 32 FFT IP kernels of fft algorithm of altera corp's employing traditional structure is greater than 1.0us.Also want more than the 1.0us with 32 FFT times that DSP does.The present invention utilizes the abundant resources such as logic, RAM (Read Access Memory), ROM (Read Only Memory) and DSP of FPGA device, adopts a kind of new construction fft algorithm, makes the FFT arithmetic speed than classic method large increase arranged.Along with improving constantly of chip integration, this its superiority of high speed FFT disposal route will be more and more obvious.And this method is easily expanded.

Description of drawings

Fig. 1: butterfly computation figure;

Fig. 2: decimation in time: 8 DFT are decomposed into 42 DFT;

Fig. 3: the parallel full fft algorithm structural drawing of basic 2 decimations in time;

Fig. 4: FFT structural drawing after the step exchange;

Fig. 5: high-speed parallel FFT Processing Structure;

The FPGA of Fig. 6: high speed FFT realizes block diagram;

Fig. 7: emulation sequential chart;

Fig. 8: FPGA operation result;

Fig. 9: theoretical output result.

Embodiment

The present invention adopts a kind of base-2 new structure fft algorithm, uses the advanced FPGA device Stratix of altera corp series EP1S25F780-5 and does simulation hardware.The Stratix device is the FPGA device of a employing high performance structures system.It combines powerful core performance, big memory bandwidth, digital signal processing (DSP) function, the FPGA of High Speed I/O performance and modular design and one.Its embedded DSP module has very high multiplying speed, when programming, can specify with device inside DSP module and generate multiplier with the MegaWizard method with VHDL, do butterfly with this multiplier, constitute FFT computing level with a plurality of butterflies, can realize FFT computing parallelization by water operation.The present invention also will adopt altera corp's Quartus II software to do logical design and time series analysis, and utilize powerful simulation hardware of Quartus II software and logic analysis function, the VHDL hardware description language is designed compiling, comprehensive, layout, wiring, download at last and generated special-purpose FFT process chip in the fpga chip.

Calculate in the coefficient that needs the computing finished in the DFT process, have considerable symmetry.By studying this symmetry, can simplify the computing in the computation process, calculate the required time of DFT thereby reduce.

Utilize

Periodicity, symmetry, reducibility three big characteristics, x (n) or X (k) sequence can be resolved into short sequence according to certain rules and carry out computing, can avoid a large amount of repetitive operations like this, improve the arithmetic speed of calculating DFT.Algorithm pattern can be divided into two big class, i.e. decimation in time fft algorithm and decimation in frequency fft algorithms.

This patent adopts base-2 decimation in time fft algorithms.If sequence length N is 2 integer power power, i.e. N=2 ^M, wherein M is a positive integer, expression power number of times.

At first sequence x (n) is decomposed into two groups, even number Xiang Weiyi group, odd term is one group, obtains the subsequence that two N/2 are ordered, that is:

x ₁(r)＝x(2r)，x ₂(r)＝x(2r+1)，(0≤r≤N/2-1) (3)

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left(n\right)W_N^{\mathrm{kn}}=\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x\left(2r\right)W_N^{2\mathrm{kr}}+\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x(2r+1)W_N^{k(2r+1)}---\left(4\right)$

X wherein ₁(k) and X ₂(k) be respectively x ₁(r) and x ₂(r) DFT.Utilize Formula (4) can be write as:

$\left\{\begin{array}{c}X\left(k\right)=X_1\left(k\right)+W_N^k{X}_2\left(k\right)\\ X(k+N/2)=X_1\left(k\right)-W_N^k{X}_2\left(k\right)\end{array}\right.,(0\&le;k\&le;N/2-1)---\left(5\right)$

Represent as shown in Figure 1 with patterned mode.Each butterfly computation needs a complex multiplication and twice complex addition computing.Adopt this method for expressing, the procedure graph of above-mentioned decomposition operation as shown in Figure 2.After decomposing, each N/2 point DFT needs N ²/ 4 complex multiplication, two N/2 point DFT need N altogether ²/ 2 complex multiplication, combinatorial operation need N/2 butterfly computation altogether, need N/2 complex multiplication, therefore need N altogether ²/ 2+N/2 ≈ N ²/ 2 complex multiplication are compared the closely operand of half of saving with direct computing.

Because N/2=2 ^M-1Still be integer, this sequence can be decomposed down always, to the last be 2 DFT till.For one 8 FFT, can obtain the computing flow graph of base-2FFT of a complete N=8 according to above-mentioned algorithm, as shown in Figure 3.

According to above-mentioned fft algorithm principle and Fig. 3, can summarize the certain law and the characteristics of base-2FFT algorithm:

1, whole FFT flow graph is made up of butterfly computation entirely, so butterfly computation is the core of FFT computing, and the specific implementation of this algorithm is exactly how to calculate successively in certain sequence to finish whole butterfly computations.Each butterfly computation needs once the addition of multiplying each other of plural number and twice plural number.For a N=2 ^MSequence, can progressively decompose at last is 2 DFT computing entirely.The FFT computing that whole N are ordered is total

Individual butterfly computation needs complex multiplication altogether Inferior, complex addition a _f=NM=Nlog ₂N time, compare operand with direct calculating and significantly reduce.

2, the input and output amount of each butterfly is not heavy mutually mutually in the flow graph, two input quantities of any one butterfly are behind butterfly computation, just lost value, do not need to preserve, therefore can realize with the location computing, promptly can leave in the unit of the same address of original storage input data through the one-level calculated result.Therefore, only need the storage unit of N plural number, both can deposit the raw data of input, can deposit intermediate result again, but also can deposit last operation result.This original position compute mode has been saved a large amount of storage unit, is a big advantage of fft algorithm, thereby effectively reduces equipment cost.

3, for location computing structure, computing finish back output as a result X (k) still deposit by natural order.And list entries x (n) owing to press the decomposition result of idol, strange time sequencing extraction one by one, has rearranged the order of depositing of sequence data, so it is to discharge by the order of sign indicating number bit reversal.So-called sign indicating number bit reversal, be exactly with binary highest significant position to the bit-reversed of least significant bit (LSB) place and must binary number.For example, during N=8, the binary code position 001 of n=1 is 100 behind the sign indicating number bit reversal, and promptly corresponding binary number is 4.In the actual operation, directly x (n) is not exported according to bit reversed order, earlier it is imported storage unit according to natural order usually, and according to fft algorithm, the operation transform that indexes then obtains the arrangement of bit reversed order, as table 1.

4, the required coefficient of butterfly computation At different levels different.Observe from top to bottom for every grade, all be with Beginning increases progressively successively by geometric progression, and the cycle repeats.For example, the computing of m level, coefficient is

1=0,1 ..., 2 ^M-1-1, totally 2 ^M-1Individual coefficient, index 1 increases 1 one by one, and the cycle repeats 2 ^M-mInferior.Required coefficient can calculate good back in advance and exists in the numerical table during calculating, fast operation like this, but need the expense internal memory.Recursion calculating successively when needed can be saved internal memory like this, but will be increased the certain calculation workload.

Table 1 yard position ordering (N=8)

This full parallel organization characteristics shown in Figure 3, if employing synchronizing sequential circuit, can realize that then each timeticks exports one group of FFT result of calculation, thereby give full play to the parallel fast processing characteristics that add flowing structure, make the FFT arithmetic speed be greatly improved.

[high speed FFT Processing Algorithm]

Introduction by the front knows that 2,8 fft algorithms of tradition base as shown in Figure 3.Digitized representation twiddle factor on the arrow

In k.Input is by the order discharging of sign indicating number bit reversal, and output is natural order.The characteristics of this structure are that each butterfly output data still is placed in the original input data storage cell, so only need 2N storage unit (data are plural forms among the FFT, and two unit of every some needs are stored, and N is that the FFT computing is counted).Its shortcoming is that same position butterfly input data addressing not at the same level is unfixing, is difficult to realize cycle control.Be difficult to Parallel Implementation when programming with FPGA, data processing speed is slow.It is all the more so when FFT counts increase.

Can find by observing the traditional structure fft algorithm, if, then can obtain structure as shown in Figure 4 two butterfly exchanges in the middle of the first order.To the further conversion of this structure, with do not bother to see me out reposition but it is stored and deposits in order of second level output, then two and then transposings of butterflies in the middle of the third level, and input arranged in order, then can constitute high-speed parallel FFT Processing Structure as shown in Figure 5.In the butterfly conversion simultaneously, its twiddle factor is also followed conversion.

The input of this structure is an order, and output is the bit code inverted order, and every grade of twiddle factor all is stored in the FPGA sheet in the ROM, and wherein the addressing of twiddle factor is crucial.By the above FFT structural drawing conversion process every grade of twiddle factor addressing rule of can deriving, the FFT (N=2 of ordering for N ^k, k is a progression), twiddle factor has

N/2 altogether, they are lined up an array that contains N/2 element by bit code inverted order form count:

I level (i=0,1,2 then ..., it is with w (0) that twiddle factor k-1) puts in order, w (1), and w (2) ...., w (2 ⁱ) repeat 2 ^K-i-1Inferior obtain.This design feature is that every grade of input, output data order are constant, and therefore every grade of geometry is fixed, and is easy to realize parallel FFT hardware configuration with FPGA, thereby obviously accelerates the FFT arithmetic speed.

[realization of FPGA hardware]

The characteristics of FPGA device are that the available hardware descriptive language carries out flexible programming to it.But utilize the function of the software simulation hardware that FPGA manufacturer provides, make hardware design flexible as software design.Shortened the system research and development cycle.Utilize JTAG (Joint Test Action Group) interface to carry out the dirigibility that ISP (In System Programmable, in-system programming) has improved system to it.Along with the raising of chip integration, not only have in the monolithic FPGA a large amount of logical blocks and also can also integrated RAM, ROM, I/O and DSP piece etc., thereby SOC (System On a Chip, SOC (system on a chip)) is become a reality.That the present invention adopts is the Stratix family chip EP1S25F780-5 of altera corp, does simulation hardware and logic analysis with QuartusII software, and will export the result and Matlab K-theoretic operation K result compares.System architecture diagram as shown in Figure 6.The design feature of native system is:

1. for improving data precision, system all uses 16 bit wides.Use data-array, three arrays of write-array and fly-array have realized the parallel processing of kernel, can calculate 32 multiple FFT in 10 clock period.Clock period was 25 nanoseconds, and therefore 32 FFT only needed for 250 nanoseconds.

2. the flowing water input and output of data have been realized.Calculating i group data simultaneously, i-1 group data FFT result is just in serial output, and i+1 group data are then just imported in serial.Because that kernel calculate to be to walk abreast, speed is fast, input tape is roomy.Native system data input rate can reach 200MHz.For the correctness of checking design, below it is carried out emulation, input signal is as follows:

x(n)＝sin(2*n*pi/4.7)+sin(2*n*pi/16.3)+w(n)?(n＝0，1，2…31) (6)

Emulated data is made of two sine wave signal superimposed noises, and signal to noise ratio (S/N ratio) is 3dB, and sinusoidal signal frequency is respectively 0.213f _sWith 0.061f _s, f _s=1 is the normalization sample frequency, and w (n) is that average is that 0 variance is 1 Gaussian noise.The system emulation waveform as shown in Figure 7.

The FFT operation result of FPGA output and the FFT notional result of calculating with Matlab are compared, shown in Fig. 8,9.What simulate signal adopted is real signal, and its frequency spectrum has symmetry, and half of only having drawn 32 simulation results among the figure is preceding 16 output results.

Claims (1)

1. realize the method that high speed FFT handles based on FPGA for one kind, it is right to it is characterized in that discrete Fourier transformation X (k) and input signal x (n) have constituted discrete Fourier transformation, all is longly to be that the complex sequences of N, sequence length N are 2 integer power power, i.e. N=2 ^M, wherein M is a positive integer, wherein k, n represent the sequence of discrete Fourier transformation X (k) and input signal x (n) respectively;

At first discrete Fourier transformation x (n) is decomposed into two groups, even number Xiang Weiyi group, odd term is one group, obtains the subsequence that two N/2 are ordered, that is:

x ₁(r)＝x(2r)，x ₂(r)＝x(2r+1)，(0≤r≤N/2-1) (1)

$X\left(k\right)=\mathrm{DFT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left(n\right)W_N^{\mathrm{kn}}=\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x\left(2r\right)W_N^{2\mathrm{kr}}+\underset{r=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}x(2r+1)W_N^{k(2r+1)}---\left(2\right)$

Utilize

Formula (2) can be write as:

$\left\{\begin{array}{c}X\left(k\right)=X_1\left(k\right)+W_N^k{X}_2\left(k\right)\\ X(k+N/2)=X_1\left(k\right)-W_N^k{X}_2\left(k\right)\end{array}\right.,(0\&le;k\&le;N/2-1)---\left(3\right)$

X wherein ₁(k) and X ₂(k) be respectively x ₁(r) and x ₂(r) DFT.

Images