Fast co-prime spectrum analysis system based on FPGA (Field Programmable Gate Array)

Abstract

The invention discloses a fast co-prime spectrum analysis system based on an FPGA (Field Programmable Gate Array), which comprises the following processing processes: initialization: setting parameters of spectrum detection, designing a convolution window filter bank according to the set parameters, performing fixed-point quantification and splicing on filter coefficients and twiddle factors required in calculation, and storing the filter coefficients and twiddle factors in a predetermined ROM (Read Only Memory); a multi-phase filtering stage; first-stage IDFT: performing first sequential transformation on the multi-phase filtering result, reading a twiddle factor, calling a complex multiplication accumulator group, completing first-stage IDFT calculation, and storing the result; a second stage of IDFT: performing second sequential transformation on the calculation result of the first stage of IDFT, reading the twiddle factor, calling the complex multiplication accumulator group, completing second-stage calculation of IDFT, performing third sequential transformation on the calculation result, and storing the result in the storage area; and in a cross-correlation stage, outputting a final power spectrum.

Description

Technical Field

[0001] The invention relates to the field of radio broadband spectrum perception, and in particular to a fast coprime spectrum analysis system based on FPGA. Background Art

[0002] Wideband spectrum sensing involves real-time, large-scale scanning and monitoring of the radio spectrum, which is crucial for improving spectrum utilization efficiency and optimizing the allocation of spectrum resources [1]. When hardware resources are limited, its detection range will be limited to a certain extent. Due to the expansion of bandwidth, if sampling is performed at the Nyquist frequency, wideband spectrum sensing will require an ADC with a high sampling rate, resulting in higher hardware costs and power consumption. In addition, Nyquist sampling of large bandwidths will result in a huge amount of sampled data, which will bring great difficulties to subsequent signal processing and will also lead to poor resource consumption and real-time performance of hardware implementation. In response to the shortcomings of the Nyquist sampling scheme, researchers have proposed a spectrum sensing scheme based on sub-Nyquist sampling. At present, the main sub-Nyquist sampling methods include compressed sensing [2] and coprime frequency undersampling [3].

[0003] At present, the main methods of compressed sensing sampling include analog information converter (AIC) [4], modulation wideband converter (MWC) [5], and multiple coset sampling [6]. Spectrum sensing methods based on compressed sensing get rid of the dependence on high-speed ADC and significantly reduce the amount of data to be processed. However, the MWC sampling method requires a high-speed pulse generation circuit and a multi-branch analog filter to complete the undersampling task, which has high hardware costs. The multiple coset sampling method requires the driving clock of the analog-to-digital converter of each branch to meet the sub-nanosecond delay relationship. At present, there are only simulation reports on the multiple coset method, and no reports on engineering implementation.

[0004] Therefore, realizing under-sampling broadband spectrum analysis with low hardware cost is an urgent problem to be solved in the radio field. The present invention intends to realize the under-sampling broadband spectrum analysis function on a general digital signal processing platform by taking the coprime under-sampling spectrum analysis method as a breakthrough.

[0005] The coprime spectrum analysis algorithm only needs to perform low-speed parallel digital signal processing on two coprime sampling channels to complete the broadband spectrum analysis task. The key operations on any channel include three basic digital signal processing operations: polyphase filtering, IFFT, and cross-correlation operations [7]. The core operations of these three types of operations are all multiplication and accumulation (MAC). Usually, the deployment platforms of high-speed signal processing algorithms include DSP processors and FPGAs. Since the DSP processor platform adopts a serial computing structure, running the coprime spectrum analysis algorithm on this type of platform cannot give full play to the parallel computing advantages of the algorithm. Compared with DSP, the current FPGA not only has more massive digital signal processing circuit resources for multiplication and accumulation calculations, but also has the advantage of high parallel processing that DSP does not have. Therefore, the present invention abandons DSP and uses FPGA hardware platform to implement

[0006] Given that the computing resources and storage space in FPGA are relatively limited, and the hardware resource requirements of the coprime spectrum analysis algorithm increase with the improvement of the engineering requirements of spectrum resolution, it is necessary to strike a balance between spectrum resolution, hardware resource consumption, and detection speed and to improve the performance of all aspects as much as possible.

[0007] Related Literature

[0008] [1] Xiao Dong, Zhang Zhengmin, Hou Yanxi, et al. Wideband spectrum sensing technology based on transposed sampling[J]. Information Countermeasures Technology, 2024, 3(06): 60-70.

[0009] [2]Donoho D L.Compressed sensing[J].IEEE Transactions on informationtheory,2006,52(4):1289-1306.

[0010] [3]Vaidyanathan PP,Pal P.Sparse sensing with co-prime samplers andarrays[J].IEEE Transactions on Signal Processing,2010,59(2):573-586.

[0011] [4]Haque T, Yazicigil RT, Pan KJL, et al. Theory and design of aquadrature analog-to-information converter for energy-efficient widebandspectrum sensing[J]. IEEE Transactions on Circuits and Systems I: Regular Papers, 2014, 62(2):527-535.

[0012] [5]Byambadorj Z, Asami K, Yamaguchi TJ, et al. High-precision sub-Nyquist sampling system based on modulated wideband converter for communication device testing [J]. IEEE Transactions on Circuits and Systems I: Regular Papers, 2021, 69 (1): 378-388.

[0013] [6]Song Z, Yang J, Zhang H, et al.

[0014] [7] Huang Xiangdong, Song Jinshui, Li Yanping. High-precision carrier frequency estimation of multi-band communication signals under undersampling[J]. Journal of Electronics & Information Technology, 2024, 46(03): 906-913. Summary of the invention

[0015] The purpose of the present invention is to design a parallel processing system based on FPGA to realize a fast coprime spectrum analysis method, and to perform multiphase filtering, IDFT and cross-correlation processing in sequence according to the coprime spectrum, to construct a multi-level high-efficiency MAC parallel digital processing system, and to adopt a built-in pipeline structure for modular design for complex multiplication and complex modulus basic operations, so that in the process of frequent multiplexing, not only the fast response of spectrum analysis is realized, but also hardware resources are saved. The technical solution is as follows:

[0016] A fast coprime spectrum analysis system based on FPGA includes the following processing steps:

[0017] Step 1, initialization: set the parameters of spectrum detection, design the convolution window filter group according to the set parameters, perform fixed-point quantization and splicing on the filter coefficients and rotation factors required in the calculation, and save them in the predetermined ROM;

[0018] Step 2, multi-phase filtering stage: write the sampled data into the multi-phase decomposition data cache, perform rate conversion on the upper and lower channel sampled data x(Nn) and x(Mn), and decompose them into M-way data streams and N-way data streams respectively; read out the multi-phase decomposition data and filter coefficients, call the complex multiplication accumulator group, complete the multi-phase filtering calculation, and save the results;

[0019] Step 3, IDFT stage 1: perform the first order transformation on the polyphase filtering result, read out the rotation factor at the same time, call the complex multiplication and accumulator group, complete the IDFT stage 1 calculation, and save the result;

[0020] Step 4, IDFT stage 2: Perform a second order transformation on the calculation results of the IDFT stage 1, read out the rotation factors at the same time, call the complex multiplication and accumulator group, complete the IDFT stage 2 calculation, perform a third order transformation on the calculation results, and save them to the storage area;

[0021] Step 5, cross-correlation stage: perform periodic extension on the IDFT calculation results of the upper and lower channels respectively, call the complex multiplier accumulator group to perform cross-correlation calculation on the extended upper and lower channel data, and save the results; when the cross-correlation calculation of all the data to be processed is completed, call the CORDIC module to perform modulus calculation and output the final power spectrum.

[0022] Further, the method of step one is as follows:

[0023] Set the equivalent Nyquist frequency of the spectrum detection upper limit, set the coprime integer pair M and N and the spectrum resolution, and design the convolution window filter group according to the parameters M and N; perform fixed-point quantization on the filter coefficients, and for the upper channel, concatenate every M data; for the lower channel, concatenate every N data, and generate the rotation factors of points M1, M2, N1, and N2 respectively;

[0024] Construct M and N complex multiplier accumulators in the upper and lower channels respectively to prepare for global calls at each stage; convert the polyphase filter coefficients H0(z),...,H M-1 (z), the lower channel polyphase filter coefficients G0(z), ..., G N-1 (z) are stored in a predetermined ROM respectively; the IDFT coefficients at point M1 and the IDFT coefficients at point M2 of the upper channel are stored in a predetermined ROM, and the IDFT coefficients at point N1 and the IDFT coefficients at point N2 of the lower channel are stored in a predetermined ROM.

[0025] Furthermore, the method of step 2 is as follows:

[0026] Multiphase decomposition data cache: In the upper channel, the coprime undersampled serial high-speed data stream x(Nn) is converted into M parallel low-speed data streams Cache into RAM; Similarly, in the lower channel, the mutually prime undersampled serial high-speed data stream x(Mn) is converted into N parallel low-speed data streams Cache to RAM;

[0027] Call the M complex multiplier accumulators on the upper channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients H m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 1, m=0,…,M-1;

[0028] Call the N complex multiplier accumulators of the lower channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients G m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 2, m=0,…,N-1.

[0029] Furthermore, the method of step three is as follows:

[0030] On the upper channel, read the filtering result data of storage area 1 It is further decomposed into M2 M1 point sub-data streams, and the M complex multiplier accumulators are correspondingly divided into M2 groups, each group containing M1 complex multiplier accumulators; for each group of complex multiplier accumulators in the M2 group, the IDFT coefficients of the corresponding M1 point are fed into, and after the multiplication and accumulation operation is performed, the corresponding transformation result is stored in the IDFT first stage output storage area 1;

[0031] In the lower channel, read the filtering result data of storage area 2 It is further decomposed into n2 n1-point sub-data streams, and the N complex multiplier accumulators are correspondingly divided into n2 groups, each group containing N1 complex multiplier accumulators; for each group of complex multiplier accumulators in the N2 groups, the corresponding N1-point IDFT coefficients are fed into, and after the multiplication and accumulation operations, the corresponding transformation results are stored in the IDFT first stage output storage area 2.

[0032] Furthermore, the method of step 4 is as follows:

[0033] In the upper channel, the transformed output data stored in the storage area 1 in the first stage of IDFT is read and further decomposed into M1 M2 point sub-data streams. Correspondingly, the M complex multiplier accumulators are divided into M1 groups, each group contains M2 complex multiplier accumulators, and for each group of complex multiplier accumulators in the M1 group, the IDFT coefficients of the corresponding M2 points are fed into the corresponding calculation results after multiplication and accumulation operations. Store in IDFT second stage output storage area 1;

[0034] Each time the data outputted from the second stage of IDFT to the storage area 1 is read out, the data inside the storage area 1 is cyclically shifted with a step length of N, and a data write operation is performed on its extension area; the extension area of ​​the storage area 1 includes N data registers, which are sequentially connected to the 1st to the Nth data registers of the storage area 1;

[0035] Whenever new data is written into IDFT phase 2 storage area 1 or an internal shift is performed, the new data in storage area 1 is synchronously written into the connected registers in the extension area; the output of IDFT phase 2 output storage area 1 after expansion is in in((.)) M is the modulo operation on M;

[0036] In the lower channel, the transformed output data of storage area 2 is read and further decomposed into N1 N2 point sub-data streams. Correspondingly, the N complex multiplier accumulators are divided into N1 groups, each group containing N2 complex multiplier accumulators; for each group of complex multiplier accumulators in the N1 group, the corresponding N2 point IDFT coefficients are fed into the corresponding N2 point IDFT coefficients, and after the multiplication and accumulation operation, the corresponding transformed result is Store in IDFT second stage output storage area 2;

[0037] After each reading of the data in the IDFT second-stage output storage area 2, the data inside the storage area 2 is cyclically shifted with a step length of M, and a data write operation is performed on its extension area; the extension area of ​​the storage area 2 includes M data registers, wherein the first N registers are sequentially connected to the 1st register to the Nth data register of the storage area 2, and the N+1th register to the Mth register in the extension area of ​​the storage area 2 are sequentially connected to the 1st register to the MNth register in the storage area 2;

[0038] Whenever new data is written into IDFT phase 2 storage area 2 or an internal shift is performed, the new data in storage area 2 is synchronously written into the connected registers in the extension area; the output of the IDFT phase 2 output storage area 2 after expansion is in in((.)) N It is the modulo operation of N.

[0039] Furthermore, the method of step five is as follows:

[0040] At the same time, the M complex multiplier accumulators of the upper channel and the N complex multiplier accumulators of the lower channel are called, and the M+N storage area data after the upper channel extension and the M+N storage area data after the lower channel extension are respectively fed into the above M+N complex multiplier accumulators, and the following cross-correlation operation is performed:

[0041]

[0042] The M+N cross-correlation operation results are stored in the cross-correlation storage area with a capacity of MN. When the spectrum number i has traversed all MN integer values, the M+N complex multiplication accumulators are continuously called. A total of Round parallel cross-correlation operation;

[0043] The data in the cross-correlation storage area with a cross-correlation output capacity of MN is read out in serial form, and the CORDIC module is called to perform a modulo operation on the cross-correlation results of MN complex values; first, the complex coordinates of the input complex data are judged. If the real input sign bit is 1 and the imaginary input sign bit is 0, it corresponds to the second quadrant, or the real input sign bit is 1 and the imaginary input sign bit is 1, it corresponds to the third quadrant. The real input and the imaginary input are inverted at the same time, and the complex coordinates are corrected from the second quadrant to the fourth quadrant or from the third quadrant to the first quadrant; if the complex coordinates of the input data are in the first quadrant or the fourth quadrant, no correction is made; set the CORDIC module The number of iterations is K, and the corrected data passes through K CORDIC units for coordinate transformation in succession; for the nth CORDIC unit, first determine the sign bit of the imaginary input. If the sign bit of the imaginary input is 1, the real output is the real input plus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input minus the imaginary input shifted right by n bits; if the sign bit of the imaginary input is 0, the real output is the real input minus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input plus the imaginary input shifted right by n bits; after 16 transformations, the real output of the last unit multiplied by the scale factor is regarded as the modulus result, and MN power spectrum values ​​|S are output one by one. xx (w i )|, i∈[0,MN-1].

[0044] The beneficial effects of the present invention are as follows:

[0045] (1) The present invention is based on FPGA, and through time division multiplexing and parallel structure design of the core computing module, the multiplier-accumulator group completes the polyphase filtering, IDFT, and cross-correlation calculations required in the system, which greatly reduces the hardware resources required for the algorithm while ensuring the calculation speed;

[0046] (2) The present invention introduces the Good-Thomas algorithm into the coprime spectrum analysis algorithm to solve the problem that the coprime spectrum analysis cannot apply the FFT based on butterfly operation in the upper and lower channels at the same time; at the same time, combined with the software and hardware design of FPGA, the calculation amount of IDFT is reduced and the calculation speed of IDFT is improved without the need for additional computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram of the coprime spectrum analysis algorithm flow chart

[0048] Figure 2 Schematic diagram of hardware structure design of fast interprime spectrum analyzer based on FPGA

[0049] Figure 3 Schematic diagram of data interaction between the IDFT second-stage output storage area and its extension area

[0050] Figure 4 Schematic diagram of the hardware structure design of the multiplication and accumulation unit

[0051] Figure 5 Schematic diagram of CORDIC module and CORDIC unit structure

[0052] Figure 6 System operation control state transition diagram

[0053] Figure 7 Schematic diagram of the panoramic spectrum generated in the simulation results and the system operation time

[0054] Figure 8 Frequency points and their position numbers in the simulation results

[0055] Fig. 9 Hardware implementation diagram of the present invention. DETAILED DESCRIPTION

[0056] The following first gives the overall architecture and calculation process of the invention, then introduces its technical principles and implementation process in detail, and finally gives a specific experiment to verify the process and performance.

[0057] The overall architecture and calculation process of the present invention are as follows:

[0058] The present invention is based on Figure 1 The coprime spectrum analysis signal processing flow shown in FIG. 1 (including a polyphase filtering module, an IDFT module and a cross-correlation module in sequence) is proposed as follows: Figure 2 The hardware structure system design based on FPGA is shown (let the prime number pair M=M1M2, N=N1N2) to realize full digital broadband spectrum analysis.

[0059] Figure 2 The FPGA hardware system, the processing process of each stage is described as follows:

[0060] Initialization: Construct M and N multiplication accumulators in the upper and lower channels respectively, so as to prepare for global calls in each stage. Set the polyphase filter coefficients H0(z),...,H M-1 (z), the lower channel polyphase filter coefficients G0(z), ..., G N-1 (z) are stored in a predetermined ROM respectively; the IDFT coefficients at point M1 and the IDFT coefficients at point M2 of the upper channel are stored in a predetermined ROM, and the IDFT coefficients at point N1 and the IDFT coefficients at point N2 of the lower channel are stored in a predetermined ROM.

[0061] Step 1: Polyphase filtering stage

[0062] (1) Multiphase decomposition data cache: In the upper channel, the coprime undersampled serial high-speed data stream x(Nn) is converted into M parallel low-speed data streams Cache into RAM; Similarly, in the lower channel, the mutually prime undersampled serial high-speed data stream x(Mn) is converted into N parallel low-speed data streams Cache to RAM.

[0063] (2) Call the M multiplication accumulators on the upper channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients H m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 1, m=0,…,M-1;

[0064] Similarly, the N multiplication accumulators of the lower channel are called to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients G m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 2, m=0,…,N-1;

[0065] Step 2: IDFT

[0066] (1) IDFT Phase 1

[0067] On the upper channel, read the filtering result data of storage area 1 It is further decomposed into M2 M1 point sub-data streams, and the M multiplier accumulators are correspondingly divided into M2 groups (each group contains M1 multiplier accumulators). For each multiplier accumulator, the corresponding M1 point IDFT coefficient is fed into it, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT first stage output storage area 1;

[0068] Similarly, in the lower channel, read the filtering result data of storage area 2 It is further decomposed into N2 N1-point sub-data streams, and the N multipliers and accumulators are correspondingly divided into N2 groups (each group contains N1 multipliers and accumulators). For each multiplier and accumulator, the corresponding N1-point IDFT coefficients are fed into each multiplier and accumulator, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT first-stage output storage area 2;

[0069] (2) IDFT Phase 2

[0070] In the upper channel, the transformed output data of storage area 1 is read and further decomposed into M1 M2 point sub-data streams. The M multiplier accumulators are divided into M1 groups (each group contains M2 multiplier accumulators). For each multiplier accumulator, the IDFT coefficients of the corresponding M2 points are fed into the multiplication and accumulation operations, and the corresponding calculation results are Stored in IDFT second stage output storage area 1. Figure 3As shown in (a), each time the data in the IDFT two-stage output storage area 1 is read out, the data inside the storage area 1 is cyclically shifted with a step length of N, and the data is written to its extension area at the same time. The extension area of ​​storage area 1 includes N data registers, which are sequentially connected to the 1st to Nth data registers of storage area 1. Whenever new data is written to the IDFT two-stage storage area 1 or an internal shift is performed, the new data in storage area 1 is synchronously written to the connected registers in the extension area. The output of the IDFT two-stage output storage area 1 after expansion is in in((.)) M It is the modulo operation of M.

[0071] Similarly, in the lower channel, the transformed output data of storage area 2 is read and further decomposed into N1 N2 point sub-data streams. The N multiplier accumulators are divided into N1 groups (each group contains N2 multiplier accumulators). For each multiplier accumulator, the corresponding N2 point IDFT coefficients are fed into the multiplication and accumulation operation, and the corresponding transformation result is Stored in IDFT second stage output storage area 2. Figure 3 As shown in (b), each time the data in the IDFT two-stage output storage area 2 is read out, the data inside the storage area 2 is cyclically shifted with a step size of M, and a data write operation is performed on its extension area. The extension area of ​​storage area 2 includes M data registers, of which the first N registers are sequentially connected to the 1st to Nth data registers in storage area 2, and the N+1th to Mth registers in the extension area are sequentially connected to the 1st to MNth registers in storage area 2. Whenever new data is written to the IDFT two-stage storage area 2 or an internal shift is performed, the new data in storage area 2 is synchronously written to the connected registers in the extension area. The output of the IDFT two-stage output storage area 2 after expansion is in in((.)) N It is the modulo operation of N.

[0072] Step 3 Cross-correlation stage

[0073] (1) Simultaneously call the M multiplier accumulators of the upper channel and the N multiplier accumulators of the lower channel, feed the M+N storage area data after the upper channel extension and the M+N storage area data after the lower channel extension into the above M+N multiplier accumulators, and perform the following cross-correlation operation:

[0074]

[0075] (2) The M+N cross-correlation operation results are stored in a cross-correlation storage area with a capacity of MN. When the spectrum number i has traversed all MN integer values, the M+N multiplication accumulators are continuously called. A total of Rounds of parallel cross-correlation operations.

[0076] (3) Read the data in the cross-correlation output storage area in serial form, and call the CORDIC module to perform a modulo operation on MN cross-correlation results (complex values). First, the complex coordinates of the input complex data are judged. If the real input sign bit is 1 and the imaginary input sign bit is 0 (corresponding to the second quadrant), or the real input sign bit is 1 and the imaginary input sign bit is 1 (corresponding to the third quadrant), the real input and the imaginary input are inverted at the same time, and the complex coordinates are corrected from the second quadrant to the fourth quadrant or from the third quadrant to the first quadrant. If the complex coordinates of the input data are in the first quadrant or the fourth quadrant, no correction is made. Considering the characteristic that the calculation accuracy of the CORDIC algorithm increases with the increase in the number of iterations, in order to ensure the calculation accuracy of the modulo operation, 16 iterations are taken here, and the corrected data passes through 16 CORDIC units for coordinate transformation. For the nth unit, first determine the sign bit of the imaginary input. If the sign bit of the imaginary input is 1, the real output is the real input plus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input minus the imaginary input shifted right by n bits; if the sign bit of the imaginary input is 0, the real output is the real input minus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input plus the imaginary input shifted right by n bits. After 16 transformations, the real output of the last unit multiplied by the proportional factor is the approximate modulus value result. The CORDIC module receives the serial complex cross-correlation results, and after calculating the delay, outputs MN power spectrum values ​​one by one |S xx (ω i )|,i∈[0,MN-1]

[0077] Next, the principles and processes of the coprime spectrum analysis algorithm and the Good-Thomas FFT algorithm involved in the present invention are introduced, and then the FPGA implementation process of the coprime spectrum analysis algorithm is proposed, and the calculation process of the calculation module of the specific module in the system is explained.

[0078] 1 Coprime Spectrum Analysis Algorithm

[0079] Set the spectrum detection upper limit equivalent Nyquist frequency F Nyq , set the coprime integer pair M and N, the spectrum resolution Δf = F Nyq / MN.

[0080] The convolution window filter group is designed according to the parameters M and N. A Hamming window of length M is convolved with a rectangular window of length M to obtain a prototype filter H(jω) with a length of 2M-1 and a passband range of [-π / M,π / M]. Using the same method, a prototype filter G(jω) with a length of 2N-1 and a passband range of [-π / N,π / N] is obtained.

[0081] Next, the calculation process of coprime spectrum analysis is carried out:

[0082] (1) Downsample the discrete signal x(n) to be analyzed by M times to obtain the signal x(Mn), and downsample it by N times to obtain the signal x(Nn); use the filter bank H (0) (jω) Perform multi-phase filtering on the signal x(Nn) to obtain the branches after filtering Using filter bank G (0) (jω) Perform multi-phase filtering on the signal x(Mn) to obtain the branches after filtering

[0083] (2) Yes Perform M-point IDFT to obtain the branches of the upper channel that have completed channelization right Perform N-point IDFT to obtain the branches of the lower channel that have completed channelization

[0084] (3) Yes and Perform cross-correlation calculation to obtain the final panoramic power spectrum S xx (w i ).

[0085] 2Good-Thomas IFFT calculation steps

[0086] Assume that the number of IFFT points is N, and decompose N to set N=N1·N2, where N1 and N2 are coprime integers. The process of N-point FFT is decomposed into first completing N2 N1-point FFTs, then completing N2 N1-point FFTs, and sequentially transforming the results to obtain the result of N-point FFT.

[0087] Determine the input index n based on N1 and N2, the rules are as follows:

[0088] n=(N2·n1+N1·n2)mod N (1)

[0089] Among them, 0≤n1≤N1-1, 0≤n2≤N2-1

[0090] Determine the middle index m based on N1 and N2, the rules are as follows:

[0091]

[0092] Determine the output index k based on N1 and N2, the rules are as follows:

[0093]

[0094] Among them, 0≤k1≤N1-1, 0≤k2≤N2-1

[0095] It means to find the value of Euler function for N2, then subtract 1 from the value of Euler function and find the remainder of it with respect to N1. Same reason.

[0096] When N=12, let N1=4, N2=3, the transform indexes of the three times in Good-Thomas FFT are shown in Table 1:

[0097] Original serial number Input index n Intermediate index m Output index k 0 0 0 0 1 3 4 4 2 6 8 8 3 9 1 9 4 4 5 1 5 7 9 5 6 10 2 6 7 1 6 10 8 8 10 2 9 11 3 3 10 2 7 7 11 5 11 11

[0098] Table 112-point Good-Thomas FFT transform index

[0099] By sequentially transforming the data using the index obtained above and numbering them in combination with the order and multiple relationship of IFFT and FFT, an N-point IFFT can be converted into a combination of N1 N2-point FFTs and N2 N1-point FFTs. The number of multiplications required is determined by N. 2 Reduced to N·(N1+N2).

[0100] 3. The calculation process of the designed complex multiplier accumulator

[0101] The designed complex multiplier accumulator structure is as follows Figure 4 The complex multiplier is a pipeline structure, which divides the complex multiplication into three stages of calculation. The three stages of calculation are performed synchronously, and registers are inserted after each stage to ensure data transmission.

[0102] The inputs of the complex multiplier are

[0103]

[0104] The output of the complex multiplier is

[0105] O=O r +O j ·j (5)

[0106] The first stage completes three addition operations, ab, cd, and c+d, and beats a, b, and c.

[0107] The second stage completes a total of three multiplication operations: (ab)d, (cd)a, and (c+d)b.

[0108] The third stage is completed r =(ab)d+(cd)a,O j =(ab)d+(c+d)b, a total of 2 addition operations, at this point the complex multiplication calculation is completed, and the real part and imaginary part of the complex multiplication result are obtained.

[0109] The workflow of the accumulator is as follows:

[0110] The value stored in the register in the accumulator is initialized to 0; then in each clock cycle, the accumulator adds the input to the data stored in the register, and then stores the result of the addition in the register instead of the previous data;

[0111] When the input data is valid, the count value is incremented. When the count value is equal to the maximum count value, the count value is cleared and the register is cleared at the same time.

[0112] 4 Calculation process of the designed CORDIC module

[0113] The designed CORDIC module structure is as follows Figure 5 The CORDIC module is a pipeline structure, including 16 serially connected CORDIC units. The input complex data undergoes 16 pseudo-rotation operations to obtain the modulus value of the input data.

[0114] The input of the CORDIC module is x = a + b·j, and the output is O real and O image The input to the nth CORDIC unit is x n =a n +b n j, the output is O n real and O n imag The CORDIC module first determines the quadrant of the complex coordinates of the signal: when the complex coordinates of the signal are in the second or third quadrant, the real and imaginary parts of the signal are inverted to make it fall in the first and fourth quadrants; if the complex coordinates of the signal are in the first or fourth quadrant, no correction is made. The data is then sent to the CORDIC unit for pipeline processing.

[0115] The nth CORDIC unit first processes the imaginary part b of the input data. n The positive or negative value is judged, and then calculation is performed and the result is output according to the judgment result.

[0116] When b n <0,

[0117]

[0118] When b n >0,

[0119]

[0120] The output of the last CORDIC unit, i.e., the 16th CORDIC unit, is multiplied by the scale factor to obtain the final modulo value result.

[0121] 5. FPGA Implementation Process of Coprime Spectrum Analysis Algorithm

[0122] The filter coefficients are quantized at fixed point, and converted from floating point numbers to fixed point signed numbers with a bit width of 16 bits. For the upper channel, every M data is concatenated; for the lower channel, every N data is concatenated. The coefficients are saved in the form of a coe file.

[0123] Generate the rotation factors of points M1, M2, N1, and N2 respectively, and convert the rotation factors from floating point numbers to fixed-point signed numbers with a bit width of 16 bits. For the upper channel, concatenate every M data. For the lower channel, concatenate every N data. Save the coefficients in the form of a coe file.

[0124] The signal to be analyzed is sampled and the collected data is saved in the data rate conversion module for subsequent processing. The clock frequency of all operations described below is 100 MHz.

[0125] The system works in the form of a state machine, and its state transition diagram is as follows Figure 6 As shown, it specifically includes the following five stages:

[0126] The system first enters the data multi-phase decomposition stage: the sampled data is stored, one data is stored in each clock cycle, and the data is stored in the multi-phase decomposition data cache after delay; the upper channel stores 2M clock cycles, and the lower channel stores 2N clock cycles, and the read enable is controlled by an independent counter. For the multi-phase decomposition data cache, within one cycle, the data is written to the register with the lowest address bit, and the value of the low address bit register is sent to the high address bit register. When writing, the column address has a higher priority than the row address; after completing all operations, the system switches to the next state.

[0127] The system enters the FIR calculation stage: first, the multi-phase decomposition data cache 1 and 2 are read, and the two modules read out one column of data each time, and the column address pointer is switched at the same time; the filter coefficient storage area 1 and 2 are read, and the data in one address in the ROM is read out each time, and the address pointer is switched at the same time; the output of the above modules is synchronously sent to the multiplication accumulator group for calculation after one clock cycle; after a calculation delay of several clock cycles, the multiplication accumulator group outputs the FIR calculation result and stores it in the filter result storage area; after completing all operations, the system switches to the next state; this process completes a total of 2 parallel complex multiplication and accumulation operations, and the number of parallel operations in one clock cycle is M+N.

[0128] The system enters the first stage of IFFT calculation: first, the filter result storage area is read, M2 and N2 data are read out from the filter result storage areas 1 and 2 respectively, and M1 times and N1 times of the data are copied and spliced ​​to form parallel data, and the address pointer is switched at the same time; the M1-IDFT coefficient storage module and the N1-IDFT coefficient storage module are read, and the data in one address in the ROM is read out each time, and the address pointer is switched at the same time; the output of the filter result storage area is sequentially transformed according to the Good-Thomas FFT input index; the interleaved data and the corresponding IDFT coefficients are sent to the multiplication and accumulator group for calculation; the multiplication and accumulator group outputs the IFFT first stage calculation result and stores it in the IFFT first stage output storage area; after completing all operations, the system switches to the next state; this process completes a total of M (M1 + N1) parallel complex multiplication and accumulation operations, and the number of parallel operations in one clock cycle is M + N.

[0129] The system enters the second stage of IFFT calculation: first, the IFFT first stage output storage area is read, and M1 and N1 data are read out from the IFFT first stage output storage area 1 and 2 respectively, and M2 times and N2 times of replication and splicing are performed to form parallel data, and the address pointer is switched when reading data; the M2-IDFT coefficient storage module and the N2-IDFT coefficient storage module are read, and the data in one address in the ROM is read out each time, and the address pointer is switched at the same time; the data in the IFFT first stage output storage area are sequentially transformed according to the Good-Thomas FFT intermediate index; the interleaved data and the corresponding IDFT coefficients are sent to the multiplier accumulator group for calculation; the multiplier accumulator group outputs the IFFT second stage calculation result, according to the Good-Thomas The FFT output index is transformed sequentially and stored in the IFFT second-stage output storage area; the IFFT second-stage output storage area 1 is extended by a period of length N, and the IFFT second-stage output storage area 2 is extended by a period of length M; after all operations, the system switches to the next state; this process completes a total of M (M2 + N2) parallel complex multiplication and accumulation operations, and the number of parallel operations in one clock cycle is M + N.

[0130] The system takes the cross-correlation calculation stage into account: read the IFFT second-stage output storage area 1, each time read out the data with a length of M+N after expansion, increment the address when reading the data, and perform a circular shift with a step length of N on the data in the register group; read the IFFT second-stage output storage area 2, each time read out the data with a length of M+N after expansion and take its conjugate value, increment the address when reading the data, and perform a circular shift with a step length of M on the data in the register group; delay 1 clock in the read command After the cycle, the output data of the expanded module 1 and module 2 are sent to the multiplier group for multiplication. After a delay of several clock cycles, the multiplier group outputs the multiplication result; the cross-correlation output storage area is read to read out the cumulative cross-correlation calculation result of the previous snapshot, and M+N data are read out each time. If the current snapshot is the first snapshot, this step is skipped; the multiplication result and the cumulative cross-correlation calculation result of the previous snapshot are sent to the adder for addition, and the addition result is output and written into the cross-correlation output storage area. This process completes a total of NM parallel complex multiplication and addition operations, and the number of parallel operations in one clock cycle is M+N.

[0131] After all snapshots are processed, the cross-correlation output storage area is read, and one complex data is read out each time. After the read command is delayed by one clock cycle, the read data is sent to the CORDIC for modulo operation. After a calculation delay of several clock cycles, the real power spectrum result is output.

[0132] 6 Validation Experiment

[0133] In order to verify the function and performance of the designed coprime spectrum analyzer, the spectrum analysis results and operation speed were experimented through logic function simulation at a clock frequency of 100Mhz.

[0134] First, the required convolution window filter group coefficients and rotation factors and other parameters are generated. The filter coefficients are quantized into 16-bit fixed-point signed numbers, and the rotation factors are quantized into 16-bit fixed-point signed numbers. The generated data is saved in the ROM in the format of a coe file.

[0135] The spectrum detection range is set to 0 to 6.006 GHz, the downsampling multiples are set to N = 44, N = 39, and the spectrum resolution Δf is set to 7 MHz. The obtained panoramic spectrum is the average result of the spectrum analysis of 20 snapshot data.

[0136] Generate test data, the signal model is:

[0137] x(n)=cos(2πf0n) (8)

[0138] in,

[0139] The signal is downsampled at M=44 points and N=39 points to obtain two undersampled signals:

[0140]

[0141] The above two signals are quantized into fixed-point signed numbers with a bit width of 16 bits and saved in ROM. When the system enters the data preprocessing stage, the input signals in the ROM are read.

[0142] The spectrum analysis results are as follows Figure 7 and Figure 8 As shown in the figure, where analyzer_data_out is the output panoramic spectrum result, and analyzer_addr_out is the position number synchronized with the panoramic spectrum. The time to complete a coprime spectrum calculation is about 33.255us, and the time required to process a single snapshot is about 1.66us. The position number of the peak in the obtained panoramic spectrum is i0=412. The frequency corresponding to the number is The error from the actual frequency is δ0=1Mhz.

[0143] Experiments have shown that the estimated frequency of the signal to be detected calculated by the designed system is close to the real frequency of the signal and has a faster calculation speed.

[0144] like Fig. 9As shown, the embodiment of the present invention implants the proposed FPGA-based fast coprime spectrum analysis system into the FPGA device, based on which the high-resolution panoramic spectrum generation of the signal to be tested is completed.

Claims

1. A fast coprime spectrum analysis system based on FPGA, including the following processing procedures: Step 1, initialization: set the parameters of spectrum detection, design the convolution window filter group according to the set parameters, perform fixed-point quantization and splicing on the filter coefficients and rotation factors required in the calculation, and save them in the predetermined ROM; Step 2, multi-phase filtering stage: write the sampled data into the multi-phase decomposition data buffer, perform rate conversion on the upper and lower channel sampled data x(Nn) and x(Mn), and decompose them into M-way data streams and N-way data streams respectively; Read out the polyphase decomposition data and filter coefficients, call the complex multiplier accumulator group, complete the polyphase filtering calculation, and save the results; Step 3, IDFT stage 1: perform the first order transformation on the polyphase filtering result, read out the rotation factor at the same time, call the complex multiplication and accumulator group, complete the IDFT stage 1 calculation, and save the result; Step 4, IDFT stage 2: Perform a second order transformation on the calculation results of the IDFT stage 1, read out the rotation factors at the same time, call the complex multiplication and accumulator group, complete the IDFT stage 2 calculation, perform a third order transformation on the calculation results, and save them to the storage area; Step 5, cross-correlation stage: perform periodic extension on the IDFT calculation results of the upper and lower channels respectively, call the complex multiplier accumulator group to perform cross-correlation calculation on the extended upper and lower channel data, and save the results; when the cross-correlation calculation of all the data to be processed is completed, call the CORDIC module to perform modulus calculation and output the final power spectrum.

2. The fast coprime spectrum analysis system according to claim 1, characterized in that: The method for step one is as follows: Set the equivalent Nyquist frequency of the spectrum detection upper limit, set the coprime integer pair M and N and the spectrum resolution, and design the convolution window filter group according to the parameters M and n; perform fixed-point quantization on the filter coefficients, and for the upper channel, concatenate every M data; for the lower channel, concatenate every N data, and generate the rotation factors of M1, M2, N1, and n2 points respectively; Construct M and N complex multiplier accumulators in the upper and lower channels respectively to prepare for global calls at each stage; convert the polyphase filter coefficients H0(z),...,H M-1 (z), the lower channel polyphase filter coefficients G0(z), ..., G N-1 (z) are stored in a predetermined ROM respectively; the IDFT coefficients at point M1 and the IDFT coefficients at point M2 of the upper channel are stored in a predetermined ROM, and the IDFT coefficients at point N1 and the IDFT coefficients at point N2 of the lower channel are stored in a predetermined ROM.

3. The fast coprime spectrum analysis system according to claim 1, characterized in that: The method for step 2 is as follows: Multiphase decomposition data cache: In the upper channel, the coprime undersampled serial high-speed data stream x(Nn) is converted into M parallel low-speed data streams Cache into RAM; Similarly, in the lower channel, the mutually prime undersampled serial high-speed data stream x(Mn) is converted into N parallel low-speed data streams Cache to RAM; Call the M complex multiplier accumulators on the upper channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients H m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 1, m=0,…,M-1; Call the N complex multiplier accumulators of the lower channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients G m (z), the filtered output result after completing the multiplication and accumulation of the two Store in storage area 2, m=0,…,N-1.

4. The fast coprime spectrum analysis system according to claim 1, characterized in that: The method for step three is as follows: On the upper channel, read the filtering result data of storage area 1 It is further decomposed into M2 M1 point sub-data streams, and the M complex multiplier accumulators are correspondingly divided into M2 groups, each group containing M1 complex multiplier accumulators; for each group of complex multiplier accumulators in the M2 group, the IDFT coefficients of the corresponding M1 point are fed into, and after the multiplication and accumulation operation is performed, the corresponding transformation result is stored in the IDFT first stage output storage area 1; In the lower channel, read the filtering result data of storage area 2 It is further decomposed into N2 N1-point sub-data streams, and the N complex multiplier accumulators are correspondingly divided into N2 groups, each group containing N1 complex multiplier accumulators; for each group of complex multiplier accumulators in the N2 group, the corresponding N1-point IDFT coefficients are fed into it, and after the multiplication and accumulation operations are performed, the corresponding transformation results are stored in the IDFT first stage output storage area 2.

5. The fast coprime spectrum analysis system according to claim 1, characterized in that: The method for step 4 is as follows: In the upper channel, the transformed output data stored in the storage area 1 in the first stage of IDFT is read and further decomposed into M1 M2 point sub-data streams. Correspondingly, the M complex multiplier accumulators are divided into M1 groups, each group contains M2 complex multiplier accumulators, and for each group of complex multiplier accumulators in the M1 group, the IDFT coefficients of the corresponding M2 points are fed into the corresponding calculation results after multiplication and accumulation operations. Store in IDFT second stage output storage area 1; Each time the data outputted from the second stage of IDFT to the storage area 1 is read out, the data inside the storage area 1 is cyclically shifted with a step length of N, and a data write operation is performed on its extension area; the extension area of ​​the storage area 1 includes N data registers, which are sequentially connected to the 1st to the Nth data registers of the storage area 1; Whenever new data is written into IDFT phase 2 storage area 1 or an internal shift is performed, the new data in storage area 1 is synchronously written into the connected registers in the extension area; the output of IDFT phase 2 output storage area 1 after expansion is in in((.)) M is the modulo operation on M; In the lower channel, the transformed output data of storage area 2 is read and further decomposed into N1 N2 point sub-data streams. Correspondingly, the N complex multiplier accumulators are divided into N1 groups, each group containing N2 complex multiplier accumulators; for each group of complex multiplier accumulators in the N1 group, the corresponding N2 point IDFT coefficients are fed into the corresponding N2 point IDFT coefficients, and after the multiplication and accumulation operation, the corresponding transformed result is Store in IDFT second stage output storage area 2; After each reading of the data in the IDFT second-stage output storage area 2, the data inside the storage area 2 is cyclically shifted with a step length of M, and a data write operation is performed on its extension area; the extension area of ​​the storage area 2 includes M data registers, wherein the first N registers are sequentially connected to the 1st register to the Nth data register of the storage area 2, and the N+1th register to the Mth register in the extension area of ​​the storage area 2 are sequentially connected to the 1st register to the MNth register in the storage area 2; Whenever new data is written into IDFT phase 2 storage area 2 or an internal shift is performed, the new data in storage area 2 is synchronously written into the connected registers in the extension area; the output of the IDFT phase 2 output storage area 2 after expansion is in in((.)) N It is the modulo operation of N.

6. The fast coprime spectrum analysis system according to claim 1, characterized in that: The method for step five is as follows: At the same time, the M complex multiplier accumulators of the upper channel and the N complex multiplier accumulators of the lower channel are called, and the M+N storage area data after the upper channel extension and the M+N storage area data after the lower channel extension are respectively fed into the above M+N complex multiplier accumulators, and the following cross-correlation operation is performed: The M+N cross-correlation operation results are stored in the cross-correlation storage area with a capacity of MN. When the spectrum number i has traversed all MN integer values, the M+N complex multiplication accumulators are continuously called. A total of Round parallel cross-correlation operation; The data in the cross-correlation storage area with a cross-correlation output capacity of MN is read out in serial form, and the CORDIC module is called to perform a modulo operation on the cross-correlation results of MN complex values; first, the complex coordinates of the input complex data are judged. If the real input sign bit is 1 and the imaginary input sign bit is 0, it corresponds to the second quadrant, or the real input sign bit is 1 and the imaginary input sign bit is 1, it corresponds to the third quadrant. The real input and the imaginary input are inverted at the same time, and the complex coordinates are corrected from the second quadrant to the fourth quadrant or from the third quadrant to the first quadrant; if the complex coordinates of the input data are in the first quadrant or the fourth quadrant, no correction is made; set the CORDIC module The number of iterations is K, and the corrected data passes through K CORDIC units for coordinate transformation in succession; for the nth CORDIC unit, first determine the sign bit of the imaginary input. If the sign bit of the imaginary input is 1, the real output is the real input plus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input minus the imaginary input shifted right by n bits; if the sign bit of the imaginary input is 0, the real output is the real input minus the imaginary input shifted right by n bits, and the imaginary output is the imaginary input plus the imaginary input shifted right by n bits; after 16 transformations, the real output of the last unit multiplied by the scale factor is regarded as the modulus result, and MN power spectrum values ​​|S are output one by one. xx (ω i )|, i∈[0,MN-1].

Images