High-precision frequency measurement method based on sparse Fourier transform

Abstract

The invention discloses a high-precision frequency measurement method based on sparse Fourier transform, which comprises the following steps of: issuing normal and reverse reflection parameters, rearranging the frequency spectrum of a sampled measured signal in an FPGA (Field Programmable Gate Array), and performing window function filtering and down-sampling FFT (Fast Fourier Transform) on the rearranged sequence to obtain frequency spectrum data; and then high-energy frequency points are removed, and Hash mapping and voting are carried out on the remaining frequency points to obtain high-precision frequency information.

Description

Technical Field

[0001] The invention belongs to the technical field of digital oscilloscopes, and more specifically, relates to a high-precision frequency measurement method based on sparse Fourier transform. Background Art

[0002] Discrete Fourier transform (DFT) is one of the most basic theoretical algorithms in signal processing. The most widely used method for calculating DFT is fast Fourier transform (FFT), whose time complexity is O(Nlog2N). Compared with direct calculation of DFT, FFT greatly simplifies the calculation process and can reduce the amount of calculation by 1 to 2 orders of magnitude, which has promoted the revolutionary progress of signal processing technology. However, with the rapid growth of data volume, the calculation speed of FFT is no longer fast enough. In order to improve the time complexity of DFT, researchers have been constantly pursuing the improvement of FFT algorithm and proposed several variants, such as in-situ FFT algorithm and split-base FFT algorithm. These variants are indeed helpful in some hardware implementations, but the complexity has not been significantly improved. So researchers began to focus on the characteristics of the signal itself. The study found that a large number of signals have "sparseness": there are only a small number of large value points in the frequency domain, and the values ​​of most of the other points are close to 0. In real life, most frequencies of many signals contribute almost nothing to the entire signal, and they all have this "sparse" characteristic. Based on the characteristic that the signal data is dense and sparse, some researchers have proposed the sparse Fourier transform algorithm (SFT), which is widely used in image and speech signals, machine learning theory, Boolean function analysis, compressed sensing, broadband spectrum sensing and other fields. Related research shows that if the data is sparse in the frequency domain, the million-point transformation can be effectively calculated. SFT meets the needs of fast and efficient Fourier transform of very large data volumes.

[0003] In the field of electronic measurement, digital oscilloscopes are commonly used equipment, and their main function is to observe and analyze time domain waveforms. With the development and progress of technology, more and more oscilloscopes have been integrated with frequency domain analysis functions. Typical mixed domain oscilloscopes can simultaneously observe time domain waveforms and spectrum distributions, and have local time-frequency joint analysis capabilities. The spectrum analysis function of mixed domain oscilloscopes is realized through software FFT operations. After a small amount of data is collected by the acquisition system, the spectrum distribution information is obtained based on the butterfly fast Fourier transform algorithm. Considering that the software takes a long time to execute FFT operations, the rate of FFT spectrum updates is slow. In addition, software FFT operations will increase software overhead, resulting in a longer cycle of the entire software operation process, which ultimately manifests as slow waveform and spectrum refresh speeds. Software-based FFT operations are mainly limited by three points: computing capacity, computing efficiency, and computing accuracy. First of all, the problem of computing capacity is considered. With the development of intelligent technology, the demand for large data volume analysis is increasing. The software end is difficult to load when analyzing and processing data that is often on the G-scale. A more efficient spectrum analysis algorithm is needed to achieve spectrum perception of massive data. Secondly, computing efficiency is also a very important factor. The software uses a process-based processing method to perform spectrum analysis, which occupies a large amount of system computing resources and affects the overall operating efficiency of the system. In addition, computing accuracy is also a very important factor; computing accuracy ultimately affects the reliability and credibility of spectrum perception; limited by the software's ability to analyze the spectrum, its analysis accuracy and accuracy are difficult to significantly improve. Therefore, there are great challenges in the software implementation of FFT operations on large amounts of data, so the trend is to implement FFT operations on large amounts of data based on hardware.

[0004] Implementing spectrum analysis based on hardware can reduce the working pressure of the system and improve the operating efficiency of the system. At present, the main spectrum analysis algorithm based on hardware is the fast Fourier transform algorithm. In today's oscilloscope design, due to the development of deep storage technology, it has a large storage depth and its data throughput is very large. With the increase of data volume, the computational complexity increases linearly, and the spectrum perception of massive data is difficult to solve based on existing solutions. For FFT calculation of large data volumes, segmented FFT or parallel FFT is currently commonly used. The former divides the collected data into multiple smaller segments, and then performs FFT calculation on each segment; the latter uses parallel processing technology to process multiple FFT calculation tasks at the same time to improve the FFT calculation speed. However, these two methods are essentially still processing all data of large data volumes, and the ability to improve the calculation speed is limited, resulting in limited spectrum analysis capabilities of the oscilloscope and affecting the frequency resolution of spectrum analysis. Summary of the invention

[0005] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a high-precision frequency measurement method based on sparse Fourier transform, which utilizes the sparsity of signals in the frequency domain to solve the problem of low resolution of large-scale spectrum analysis.

[0006] To achieve the above-mentioned object of the invention, the present invention provides a high-precision frequency measurement method based on sparse Fourier transform, characterized in that it comprises the following steps:

[0007] (1) Power on the digital oscilloscope and connect it to the signal to be measured;

[0008] (2) The host computer generates a set of forward mapping parameter sequences σ with a length of L, and the corresponding reverse mapping parameter sequence σ -1 ; Then send the forward and reverse mapping parameter sequences to the sparse Fourier module;

[0009] (3) Set the maximum number of iterations of the sparse Fourier transform to L, and initialize the current number of iterations to l = 1, l = 1, 2, ..., L;

[0010] (4) Performing sparse Fourier transform on the measured signal through a cyclic iteration method;

[0011] (4.1) At the first iteration, the acquisition system of the digital oscilloscope samples the measured signal and obtains a time domain sequence of length N, which is recorded as x l (n), n=0,1,2,…,N-1;

[0012] (4.2), read the lth positive mapping parameter σ in the positive mapping parameter sequence σ l , then for the time domain series x l (n) do σ l Uniform extraction of the multiples to obtain the sequence after time-order rearrangement

[0013]

[0014] Among them, δ is the starting address of the extraction;

[0015] (4.3), sequence Perform window function filtering;

[0016] (4.3.1), define the flat window function y(n);

[0017]

[0018] in, is a convolution operation, p(n) is a rectangular window function, and satisfies: μ is the width of the rectangular window; g(n) is a Gaussian function that satisfies: b is the mean of the Gaussian function, a is the variance of the Gaussian function, and x is the Gaussian function variable;

[0019] (4.3.2), based on fixed-point multiplication Multiply it point by point with the window function y(n) to get the sequence

[0020] (4.4), for the sequence Perform decimation FFT processing;

[0021] (4.4.1) and the sequence with sampling rate B Perform downsampling and extract to get the sequence

[0022]

[0023] (4.4.2) for the sequence Perform FFT transformation to obtain the spectrum X l (k);

[0024]

[0025] Wherein, k represents the frequency index value of the output spectrum;

[0026] (4.5), frequency point hash mapping;

[0027] (4.5.1), based on spectrum X l (k), filter out the energy exceeding the threshold A Thres The frequency points, and then record the frequency index number k of these frequency points l and amplitude value A(k l );

[0028] (4.5.2), read the reverse mapping parameter sequence σ -1 The lth positive mapping parameter in Then perform hash mapping on each recorded frequency point to obtain a frequency point with higher frequency accuracy.

[0029]

[0030] Where N is the number of sparse Fourier transform points;

[0031] (4.5.3) Record frequency The corresponding amplitude Then make

[0032] (4.6) Determine whether the current number of iterations l reaches the maximum number of iterations L. If l < L, set l = l + 1, and then return to step (4.1) for the next round of iterations; otherwise, go to step (5);

[0033] (5) Frequency voting;

[0034] (5.1) All the frequency points obtained after L rounds of iterations Carry out a merger;

[0035] If a frequency If it appears more than once, all the amplitudes Accumulate and then keep a frequency index number The corresponding amplitude is recorded as Obtain a set of frequency index vectors with non-repeated frequency index numbers and corresponding cumulative amplitudes;

[0036] (5.2) Traverse each frequency index number in the frequency index vector And the corresponding cumulative amplitude Record amplitude Greater than L*A thres The frequency index number is denoted by k v , the corresponding amplitude is recorded as A(k v );

[0037] (5.3), the frequency index k v Convert to real frequency value f;

[0038]

[0039] Among them, F s is the sampling rate of the measured signal;

[0040] (6) Spectrum visualization;

[0041] The spectrum distribution (k v ,A(k v )) is displayed on the screen in real time.

[0042] The object of the invention of the present invention is achieved in this way:

[0043] The high-precision frequency measurement method based on sparse Fourier transform of the present invention first sends forward and reverse mapping parameters, then rearranges the spectrum of the sampled measured signal in FPGA, and obtains spectrum data through window function filtering and down-sampling FFT transformation of the rearranged sequence; then removes high-energy frequency points, and performs hash mapping and voting on the remaining frequency points to obtain high-precision frequency information.

[0044] At the same time, the high-precision frequency measurement method based on sparse Fourier transform of the present invention also has the following beneficial effects:

[0045] (1) The present invention greatly reduces the computational complexity of spectrum analysis by rearranging the signal spectrum and detecting the frequency point hash mapping, while improving the frequency resolution of large-scale spectrum analysis;

[0046] (2) The present invention reduces the amount of calculation for spectrum analysis by performing FFT on the extracted signal;

[0047] (3) The present invention improves the reliability of spectrum identification through multiple rounds of hash mapping and voting. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a principle block diagram of the high-precision frequency measurement based on sparse Fourier transform of the present invention;

[0049] Figure 2 It is a flow chart of the high-precision frequency measurement method based on sparse Fourier transform of the present invention;

[0050] Figure 3 It is a schematic diagram of the time-division basket when downsampling is performed. DETAILED DESCRIPTION

[0051] The specific implementation of the present invention is described below in conjunction with the accompanying drawings so that those skilled in the art can better understand the present invention. It should be noted that in the following description, when the detailed description of known functions and designs may dilute the main content of the present invention, these descriptions will be omitted here.

[0052] Example

[0053] In this embodiment, if Figure 1 As shown, the high-precision frequency measurement method based on sparse Fourier transform of the present invention mainly includes three steps of binning, reconstruction and iteration. By performing multiple rounds of binning and reconstruction iterations on the measured signal x(n), high-precision frequency data x(k) is obtained.

[0054] Next we combine Figure 2 , the specific operation process is described in detail as follows:

[0055] (1) Power on the digital oscilloscope and connect it to the signal to be measured;

[0056] (2) The host computer generates a set of forward mapping parameter sequences σ with a length of L, and the corresponding reverse mapping parameter sequence σ -1 ; Then send the forward and reverse mapping parameter sequences to the sparse Fourier module;

[0057] In this embodiment, the positive mapping parameter sequence σ is a prime number starting from 3, including 3, 5, 7, 11, 13, 17, 19, 23, 29, ..., and the length thereof depends on the number of iterations parameter L of the sparse Fourier.

[0058] The inverse mapping parameter sequence σ -1 is achieved by mapping each parameter σ in the parameter sequence σ l Perform modular inverse operation to obtain the inverse mapping parameter σl -1 ;

[0059]

[0060] Wherein, N is the number of sparse Fourier transform points, and N=65536;

[0061] In this embodiment, assuming that the total number of iterations of the sparse Fourier module is set to L = 8, the forward mapping parameter sequence σ = {3, 5, 7, 11, 13, 17, 19, 23}, and the reverse mapping parameter sequence σ -1 ={43691,52429,28087,35747,20165,61681,51739,14247};

[0062] (3) Set the maximum number of iterations of the sparse Fourier transform to L, and initialize the current number of iterations to l = 1, l = 1, 2, ..., L;

[0063] (4) Performing sparse Fourier transform on the measured signal through a cyclic iteration method;

[0064] (4.1) At the first iteration, the acquisition system of the digital oscilloscope samples the measured signal and obtains a time domain sequence of length 65536, which is recorded as x l (n), n = 0, 1, 2, ..., 65535;

[0065] (4.2), read the lth positive mapping parameter σ in the positive mapping parameter sequence σ l , then for the time domain series x l (n) do σ l Uniform extraction of the multiples to obtain the sequence after time-order rearrangement

[0066]

[0067] Among them, δ is the starting address of the extraction, and δ=0;

[0068] In this embodiment, when the first round of iteration is performed, the first positive mapping parameter σ1=3 is read, and the corresponding negative mapping parameter is

[0069] (4.3), sequence Perform window function filtering;

[0070] (4.3.1), define the flat window function y(n);

[0071]

[0072] in, is a convolution operation, p(n) is a rectangular window function, and satisfies: μ is the width of the rectangular window, μ=128; g(n) is a Gaussian function that satisfies: b is the mean of the Gaussian function, a is the variance of the Gaussian function, b=0,

[0073] (4.3.2), based on fixed-point multiplication Multiply it point by point with the window function y(n) to get the sequence

[0074] (4.4), for the sequence Perform decimation FFT processing;

[0075] (4.4.1) and the sequence with sampling rate B Perform downsampling and extract to get the sequence

[0076]

[0077] In this embodiment, if Figure 3 As shown in the figure, after the sequence is downsampled and decimated, the N high-frequency frequency points are mapped to B low-frequency baskets. Subsequently, only the FFT operation of the decimated data is required, and the frequency points in the B baskets are detected, which greatly reduces the complexity of the operation.

[0078] (4.4.2) for the sequence Perform FFT transformation to obtain the spectrum X l (k);

[0079]

[0080] Wherein, k represents the frequency index value of the output spectrum;

[0081] (4.5), frequency point hash mapping;

[0082] (4.5.1), based on spectrum X l (k), filter out the energy exceeding the threshold A Thres The frequency points, and then record the frequency index number k of these frequency points l and amplitude value A(k l );

[0083] (4.5.2), read the reverse mapping parameter sequence σ -1 The lth positive mapping parameter in Then perform hash mapping on each recorded frequency point to obtain a frequency point with higher frequency accuracy.

[0084]

[0085] Where N is the number of sparse Fourier transform points;

[0086] (4.5.3) Record frequency The corresponding amplitude Then make

[0087] (4.6) Determine whether the current number of iterations l reaches the maximum number of iterations L. If l < L, set l = l + 1, and then return to step (4.1) for the next round of iterations; otherwise, go to step (5);

[0088] (5) Frequency voting;

[0089] (5.1) All the frequency points obtained after L rounds of iterations Carry out a merger;

[0090] If a frequency If it appears more than once, all the amplitudes Accumulate and then keep a frequency index number The corresponding amplitude is recorded as Obtain a set of frequency index vectors with non-repeated frequency index numbers and corresponding cumulative amplitudes;

[0091] (5.2) Traverse each frequency index number in the frequency index vector And the corresponding cumulative amplitude Record amplitude Greater than L*A thres The frequency index number is denoted by k v , the corresponding amplitude is recorded as A(k v );

[0092] (5.3), the frequency index k v Convert to real frequency value f;

[0093]

[0094] Among them, F s is the sampling rate of the measured signal, and this embodiment corresponds to F s =20×10 9 ;

[0095] (6) Spectrum visualization;

[0096] The spectrum distribution (k v ,A(k v )) is displayed on the screen in real time.

[0097] Although the above describes the illustrative specific embodiments of the present invention to facilitate those skilled in the art to understand the present invention, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations using the concept of the present invention are protected.

Claims

1. A high-precision frequency measurement method based on sparse Fourier transform, characterized in that: The following steps are involved: (1) Power on the digital oscilloscope and connect it to the signal to be measured; (2) The host computer generates a set of forward mapping parameter sequences σ with a length of L, and the corresponding reverse mapping parameter sequence σ -1 ; Then send the forward and reverse mapping parameter sequences to the sparse Fourier module; (3) Set the maximum number of iterations of the sparse Fourier transform to L, and initialize the current number of iterations to l = 1, l = 1, 2, ..., L; (4) Performing sparse Fourier transform on the measured signal through a cyclic iteration method; (4.1) At the first iteration, the acquisition system of the digital oscilloscope samples the measured signal and obtains a time domain sequence of length N, which is recorded as x l (n), n=0,1,2,…,N-1; (4.2), read the lth positive mapping parameter σ in the positive mapping parameter sequence σ l , then for the time domain series x l (n) do σ l Uniform extraction of the multiples to obtain the sequence after time-order rearrangement Among them, δ is the starting address of the extraction; (4.3), sequence Perform window function filtering; (4.3.1), define the flat window function y(n); in, is a convolution operation, p(n) is a rectangular window function, and satisfies: μ is the width of the rectangular window; g(n) is a Gaussian function that satisfies: b is the mean of the Gaussian function, a is the variance of the Gaussian function, and x is the Gaussian function variable; (4.3.2), based on fixed-point multiplication Multiply it point by point with the window function y(n) to get the sequence (4.4), for the sequence Perform decimation FFT processing; (4.4.1) and the sequence with sampling rate B Perform downsampling and extract to obtain the sequence (4.4.2), for the sequence Perform FFT transformation to obtain the spectrum X l (k); Wherein, k represents the frequency index value of the output spectrum; (4.5), frequency point hash mapping; (4.5.1), based on spectrum X l (k), filter out the energy exceeding the threshold A Thres The frequency points, and then record the frequency index number k of these frequency points l and amplitude value A(k l ); (4.5.2), read the reverse mapping parameter sequence σ -1 The lth positive mapping parameter in Then perform hash mapping on each recorded frequency point to obtain a frequency point with higher frequency accuracy. Where N is the number of sparse Fourier transform points; (4.5.3) Record frequency The corresponding amplitude Then make (4.6) Determine whether the current number of iterations l reaches the maximum number of iterations L. If l < L, set l = l + 1, and then return to step (4.1) for the next round of iterations; otherwise, go to step (5); (5) Frequency voting; (5.1) All the frequency points obtained after L rounds of iterations Carry out a merger; If a frequency If it appears more than once, all the amplitudes Accumulate and then keep a frequency index number The corresponding amplitude is recorded as Obtain a set of frequency index vectors with non-repeated frequency index numbers and corresponding cumulative amplitudes; (5.2) Traverse each frequency index number in the frequency index vector And the corresponding cumulative amplitude Record amplitude Greater than L*A thres The frequency index number is denoted by k v , the corresponding amplitude is recorded as A(k v ); (5.3), the frequency index k v Convert to real frequency value f; Among them, F s is the sampling rate of the measured signal; (6) Spectrum visualization; The spectrum distribution (k v ,A(k v )) is displayed on the screen in real time.

2. The high-precision frequency measurement method based on sparse Fourier transform according to claim 1, characterized in that: The positive mapping parameter sequence σ is a prime number starting from 3, including 3, 5, 7, 11, 13, 17, 19, 23, 29, ..., and the sequence length depends on the number of iterations L of the sparse Fourier transform.

3. The high-precision frequency measurement method based on sparse Fourier transform according to claim 1, characterized in that: The inverse mapping parameter sequence σ -1 is achieved by mapping each parameter σ in the parameter sequence σ l Perform modular inverse operation to obtain the inverse mapping parameters Among them, l represents the parameter number and N is the number of sparse Fourier transform points.

Images