A fast coprime spectrum analysis system based on FPGA

Abstract

The invention discloses a fast coprime spectrum analysis system based on FPGA. The system comprises the following processing procedures: initialization: setting spectrum detection parameters, designing a convolution window filter group according to the set parameters, performing fixed-point quantization and splicing on filter coefficients and rotation factors required in calculation, and storing them in a predetermined ROM; a polyphase filtering stage; an IDFT stage 1: performing a first sequential transformation on a polyphase filtering result, reading out the rotation factors at the same time, calling a complex multiplication accumulator group, completing IDFT stage 1 calculation, and storing the result; an IDFT stage 2: performing a second sequential transformation on the IDFT stage 1 calculation result, reading out the rotation factors at the same time, calling a complex multiplication accumulator group, completing IDFT stage 2 calculation, performing a third sequential transformation on the calculation result, and storing it in a storage area; and a cross-correlation stage, outputting a final power spectrum.

Description

Technical Field

[0001] The present 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. Currently, 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 have broken away from the dependence on high-speed ADCs and significantly reduced 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, resulting in high hardware costs. The multiple coset sampling method requires that the driving clock of the analog-to-digital converter of each branch meet the sub-nanosecond delay relationship. Currently, there are only simulation reports on the multiple coset method, and no reports on engineering implementation.

[0004] Therefore, realizing undersampled broadband spectrum analysis with low hardware cost is an urgent problem to be solved in the field of radio. This invention intends to realize undersampled broadband spectrum analysis function on a general digital signal processing platform by using the coprime undersampled 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 the limited computing resources and storage space in FPGAs, and the increasing hardware resource requirements of coprime spectrum analysis algorithms as the engineering requirements for spectral resolution increase, it is necessary to strike a balance between spectral resolution, hardware resource consumption, and detection speed, and maximize performance in all aspects.

[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 present invention aims to design a parallel processing system based on an FPGA to implement a fast coprime spectrum analysis method. Based on the coprime spectrum, the system sequentially performs polyphase filtering, IDFT, and cross-correlation, creating a multi-stage, high-efficiency MAC parallel digital processing system. Furthermore, the system employs a modular design with an internal pipeline structure for basic operations such as complex multiplication and complex modulo. This not only achieves rapid spectrum analysis response during frequent reuse, but also conserves hardware resources. 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 bank 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 to the predetermined ROM;

[0018] Step 2, polyphase filtering stage: write the sampled data into the polyphase 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 the polyphase decomposition data and filter coefficients, call the complex multiplier accumulator group, complete the polyphase filtering calculation, and save the results;

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

[0020] Step 4, IDFT Phase 2: Perform a second sequential transformation on the IDFT Phase 1 calculation results, read out the rotation factors at the same time, call the complex multiplier accumulator group, complete the IDFT Phase 2 calculation, perform a third sequential 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] Furthermore, 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 bank based on the parameters M and N; perform fixed-point quantization on the filter coefficients, and for the upper channel, splice every M data; for the lower channel, splice every N data to 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 in each stage; 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 of the M1 point and the IDFT coefficients of the M2 point of the upper channel are stored in the predetermined ROM, and the IDFT coefficients of the N1 point and the IDFT coefficients of the N2 point of the lower channel are stored in the predetermined ROM.

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

[0026] Polyphase decomposition data buffer: 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 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 multiplication and accumulation units 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 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 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 divided into M2 groups, each group containing M1 complex multiplier accumulators; for each complex multiplier accumulator in the M2 group, 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 stage 1 output storage area;

[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 it, and after the multiplication and accumulation operations, the corresponding transformation results are stored in the IDFT stage 1 output storage area 2.

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

[0033] In the upper channel, the transformed output data stored in the storage area 1 of the IDFT stage is read and further decomposed into M1 M2 point sub-data streams. 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 the multiplication and accumulation operation. Store in IDFT second stage output storage area 1;

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

[0035] Whenever new data is written into the IDFT second stage 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 the IDFT second stage 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 and N2 point sub-data streams. The N complex multiplier accumulators are divided into N1 groups, each group contains 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 stage 2 output storage area 2;

[0037] Each time the data in the IDFT second-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 ​​the storage area 2 includes M data registers, of which 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 the IDFT second-stage storage area 2 or an internal shift is performed, the new data in the storage area 2 is synchronously written into the connected registers in the extension area; the output of the IDFT second-stage 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-way storage area data after the upper channel extension and the M+N-way storage area data after the lower channel extension are fed into the above M+N complex multiplier accumulators respectively, 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 and accumulator 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 if 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. For the nth CORDIC unit, the sign bit of the imaginary input is first determined. 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 is multiplied by the scale factor and is regarded as the modulus value 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, greatly reducing 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 coprime spectrum analysis cannot apply FFT based on butterfly operation to the upper and lower channels at the same time; at the same time, combined with the FPGA hardware and software design, the calculation amount of IDFT is reduced without the need for additional computing resources, and the calculation speed of IDFT is improved. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0048] Figure 2 Schematic diagram of the hardware structure design of the fast coprime 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] Figure 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 specific experiments 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 (including the polyphase filtering module, IDFT module and 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 of the M1 point and the IDFT coefficients of the M2 point of the upper channel are stored in the predetermined ROM, and the IDFT coefficients of the N1 point and the IDFT coefficients of the N2 point of the lower channel are stored in the predetermined ROM.

[0061] Step 1: Polyphase filtering stage

[0062] (1) Multi-phase 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 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 and accumulation devices 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 the multiplication and accumulation of the two Store in storage area 1, m = 0, ..., M-1;

[0064] Similarly, call the N multiplication and accumulation devices 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 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 divided into M2 groups (each group contains M1 multiplier accumulators). For each multiplier accumulator, the IDFT coefficient of the corresponding M1 point is fed into it, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT stage 1 output storage area;

[0068] Similarly, in the lower channel, read the filtered result data of storage area 2 It is further decomposed into N2 N1-point sub-data streams, and the N multiplier accumulators are divided into N2 groups (each group contains N1 multiplier accumulators). Each multiplier accumulator is fed with the corresponding N1-point IDFT coefficients, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT stage 1 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 and M2 point sub-data streams. The M multiplication accumulators are divided into M1 groups (each group contains M2 multiplication accumulators). For each multiplication accumulator, the IDFT coefficients of the corresponding M2 points are fed into the multiplication accumulator. After the multiplication and accumulation operation, 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 second-stage output storage area 1 is read out, the data inside the storage area 1 is cyclically shifted with a step size 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 second-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 second-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 and 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 converted into Stored in IDFT second stage output storage area 2. Figure 3 As shown in (b), each time the data in the IDFT second-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 the data is written to its extension area at the same time. 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 second-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 second-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 multiplication accumulators of the upper channel and the N multiplication accumulators of the lower channel, feed the M+N-way storage area data after the upper channel extension and the M+N-way storage area data after the lower channel extension into the above M+N multiplication accumulators, and perform the following cross-correlation operation:

[0074]

[0075] (2) Store the M+N cross-correlation operation results into the cross-correlation storage area with a capacity of MN. When the spectrum number i has traversed all MN integer values, continue to call these M+N multiplication accumulators, and 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 the MN cross-correlation results (complex values). First, judge the complex coordinates of the input complex data. If the real input sign bit is 1 and the imaginary input sign bit is 0 (corresponding to the second quadrant), or if the real input sign bit is 1 and the imaginary input sign bit is 1 (corresponding to the third quadrant), invert the real and imaginary inputs at the same time, and correct the complex coordinates 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 number of iterations, in order to ensure the calculation accuracy of the modulo operation, 16 iterations are adopted 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 scale factor is the approximate modulus value result. The CORDIC module receives the serial complex cross-correlation results and outputs MN power spectrum values ​​one by one after calculating the delay. xx (ω i )|,i∈[0,MN-1]

[0077] Then, the principles and processes of the coprime spectrum analysis algorithm and 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 bank 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) The discrete signal x(n) to be analyzed is downsampled M times to obtain the signal x(Mn), and downsampled N times to obtain the signal x(Nn); the filter bank H is used (0) (jω) Perform polyphase filtering on the signal x(Nn) to obtain the branches after filtering Use filter bank G (0) (jω) Perform polyphase 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 the channelization right Perform N-point IDFT to obtain the branches of the lower channel that have completed the 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. Decompose N and let 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 obtained 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] The output index k is determined 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 Euler function value for N2, then subtract 1 from the Euler function value and find the remainder of the Euler function with respect to N1. Same thing.

[0096] When N=12, let N1=4, N2=3, and the transform index of the third order in Good-Thomas FFT is shown in Table 1:

[0097] Original serial number Input index n Middle 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 complex multiplication and accumulation unit designed by the Institute

[0101] The designed complex multiplier accumulator structure is as follows Figure 4 The complex multiplier has 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 a total of 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 register in the accumulator is initialized to 0. Then, in each clock cycle, the accumulator adds the input to the data in the register and stores the result 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 to zero and the register is cleared to zero 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 of 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 signal's complex coordinates. If the signal's complex coordinates are in the second or third quadrant, the real and imaginary parts of the signal are inverted to bring it within the first or fourth quadrants. If the signal's complex coordinates are in the first or fourth quadrant, no correction is made. The data is then fed into the CORDIC unit for pipeline processing.

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

[0116] When b n <0,

[0117]

[0118] When b n >0,

[0119]

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

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

[0122] Quantize the filter coefficients, converting them from floating-point numbers to 16-bit fixed-point signed numbers. For the upper channel, concatenate every M data points; for the lower channel, concatenate every N data points. Save the coefficients as a coe file.

[0123] Generate the rotation factors for points M1, M2, N1, and N2, respectively, and convert them from floating-point numbers to 16-bit fixed-point signed numbers. For the upper channel, concatenate every M data points. For the lower channel, concatenate every N data points. Save the coefficients as 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 multiphase decomposition phase: sampled data is stored, one piece of data per clock cycle, and after delay, it is stored in the multiphase decomposition data cache. The upper channel stores data for 2M clock cycles, while the lower channel stores data for 2N clock cycles. Read enable is controlled by an independent counter. For the multiphase decomposition data cache, within a cycle, data is written to the register with the lowest address bit, while the value of the low-address register is simultaneously transferred to the high-address register. During writes, the column address takes precedence over the row address. After all operations are completed, 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 and accumulator group for calculation after one clock cycle; after a calculation delay of several clock cycles, the multiplication and 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, a read operation is performed on the filter result storage area, and 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 respectively, and the address pointer is switched at the same time; a read operation is performed on the M1-IDFT coefficient storage module and the N1-IDFT coefficient storage module, 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 areas 1 and 2 respectively, and M2 times and N2 times of the data are copied and spliced ​​to form parallel data respectively, and the address pointer is switched when reading the 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 IFFT first stage output storage area is 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, and 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 includes the cross-correlation calculation stage: the IFFT second stage output storage area 1 is read, and each time the data with the expanded length of M+N is read out, the address is incremented when reading the data, and the data in the register group is cyclically shifted with a step length of N; the IFFT second stage output storage area 2 is read, and each time the data with the expanded length of M+N is read out and its conjugate value is taken, the address is incremented when reading the data, and the data in the register group is cyclically shifted with a step length of M; the read command is delayed by 1 clock After the cycle, the output data of the expanded modules 1 and 2 are fed into 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 the accumulated cross-correlation calculation results of the previous snapshot, with M+N data read out each time. If the current snapshot is the first snapshot, this step is skipped. The multiplication result and the accumulated cross-correlation calculation results of the previous snapshot are fed into the adder for addition. The addition result is output and written to the cross-correlation output storage area. This process completes a total of NM parallel complex multiplications and additions, 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 Verification 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 bank coefficients and twiddle factors are generated. The filter coefficients are quantized to 16-bit fixed-point signed numbers, and the twiddle factors are quantized to 16-bit fixed-point signed numbers. The generated data is saved in ROM as a CoE file.

[0135] The spectrum detection range was set to 0 to 6.006 GHz, the downsampling multiples were set to N = 44 and N = 39, and the spectrum resolution Δf was 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, 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 required 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 true frequency of the signal, and it has a faster calculation speed.

[0144] like Figure 9As shown, the embodiment of the present invention implants the proposed FPGA-based fast coprime spectrum analysis system into the FPGA device, thereby completing the generation of high-resolution panoramic spectrum of the signal to be measured.

Claims

1. A fast coprime spectrum analysis system based on FPGA, including the following processing steps: Step 1, initialization: Set the parameters of spectrum detection, design the convolution window filter bank 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 to the predetermined ROM. The method is as follows: Set the equivalent Nyquist frequency of the spectrum detection upper limit and set the coprime integer pair and and spectral resolution, according to the parameters 、 Design a convolution window filter bank; perform fixed-point quantization on the filter coefficients. For the upper channel, each For the lower channel, each The data are spliced ​​and generated respectively 、 、 、 The rotation factor of the point; Construct M and N complex multiplier accumulators in the upper and lower channels respectively to prepare for global calls at each stage; , lower channel polyphase filter coefficients Store them in the predetermined ROM respectively; store the IDFT coefficients of point M1 and point M2 of the upper channel in the predetermined ROM, and store the IDFT coefficients of point N1 and point N2 of the lower channel in the predetermined ROM; Step 2, polyphase filtering stage: write the sampled data into the polyphase decomposition data buffer, and perform the upper and lower channel sampled data and Perform rate conversion and decompose into M-way data stream and N-way data stream 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. The method is as follows: Polyphase decomposition data buffer: In the upper channel, the coprime undersampled serial high-speed data stream Convert to M parallel low-speed data streams ,…, , cached in RAM; similarly, in the lower channel, the coprime undersampled serial high-speed data stream Convert to N parallel low-speed data streams , cached to RAM; Call the M complex multiplication and accumulation units on the upper channel to feed the corresponding low-speed data stream and the corresponding polyphase filter coefficients , the filtered output result after the multiplication and accumulation of the two is 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 , the filtered output result after the multiplication and accumulation of the two is Store in storage area 2, m=0,…,N-1; Step 3, IDFT stage 1: Perform the first sequential transformation on the polyphase filtering result, read out the rotation factors at the same time, call the complex multiply accumulator group, complete the IDFT stage 1 calculation, and save the result. The method is as follows: On the upper channel, read the filtering result data of storage area 1 , which is further decomposed into indivual The sub-data stream of the point corresponds to dividing the M complex multiplier accumulators into Groups, each containing complex multiplier accumulators; for Each complex multiplier accumulator in the group is fed with the IDFT coefficient of the corresponding M1 point, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT stage 1 output storage area; In the lower channel, read the filtering result data of storage area 2 , which is further decomposed into indivual The sub-data stream of the point corresponds to dividing N complex multiplier accumulators into Groups, each containing complex multiplier accumulators; for Each complex multiplier accumulator in the group is fed with the corresponding IDFT coefficient of N1 points, and after the multiplication and accumulation operation, the corresponding transformation result is stored in the IDFT stage 1 output storage area 2; Step 4, IDFT Phase 2: Perform a second sequential transformation on the results of the IDFT Phase 1 calculation, read out the rotation factors at the same time, call the complex multiplier accumulator group, complete the IDFT Phase 2 calculation, perform a third sequential transformation on the calculation results, and save them to the storage area. The method is as follows: In the upper channel, the transformed output data stored in storage area 1 in the IDFT stage is read and further decomposed into indivual The sub-data stream of the point corresponds to dividing the M complex multiplier accumulators into Groups, each containing complex multiplier accumulators, for Each complex multiplication and accumulation device in the group is fed with the IDFT coefficient of the corresponding M2 point, and after the multiplication and accumulation operation, the corresponding calculation result is Store in IDFT second stage output storage area 1; Each time the data output from the second stage of IDFT to storage area 1 is read out, the data inside storage area 1 is cyclically shifted with a step size of N, and a data write operation is performed on its extension area; the extension area of ​​storage area 1 includes N data registers, which are sequentially connected to the first to the Nth data registers of storage area 1; Whenever new data is written into the IDFT second stage 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 the IDFT second stage 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 indivual The sub-data stream of the point corresponds to dividing N complex multiplier accumulators into Groups, each containing complex multiplier accumulators; for Each complex multiplication and accumulation group in the group is fed with the corresponding IDFT coefficients of N2 points, and after the multiplication and accumulation operation, the corresponding transformation result is converted into Store in IDFT stage 2 output storage area 2; Each time the data in the IDFT second-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 ​​the storage area 2 includes M data registers, of which 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 the IDFT second-stage storage area 2 or an internal shift is performed, the new data in the storage area 2 is synchronously written into the connected registers in the extension area; the output of the IDFT second-stage output storage area 2 after expansion is ,in in((.)) N is the modulo operation on N; 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; after 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. The method 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-way storage area data after the upper channel extension and the M+N-way storage area data after the lower channel extension are fed into the above M+N complex multiplier accumulators respectively, 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. After traversing all MN integer values, continue to call this M+N complex multiplication accumulator, 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 if 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 of the block 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 is multiplied by the proportional factor, which is regarded as the modulus value result, and MN power spectrum values ​​are output one by one. .

Images