Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a real-time processing method for point-by-point scanning of satellite-borne VHF transient signals based on a Goertzel algorithm, the method is based on the Goertzel algorithm, a digital multichannel short-time Fourier transform is combined with a point-by-point sliding window scanning mode, the window width is fixed and calculated, the step length is selected as a sampling point, the received signals are subjected to point-by-point multichannel narrow-band channel spectrum analysis, the VHF transient events passing through an ionized layer are captured, the Goertzel algorithm can calculate while receiving the signals, all data participating in calculation are not required to be prepared, the real-time performance is good, the repeated calculation can be avoided by utilizing the sliding window recursive calculation, and the calculation speed is further accelerated.
In order to realize the purpose of the invention, the invention provides the following technical scheme:
a satellite-borne VHF transient signal point-by-point scanning real-time processing method based on a Goertzel algorithm comprises the following steps:
(1) broadband sampling is carried out on the VHF transient signal to obtain an original time domain waveform of a detection frequency band;
(2) calculating Fourier components of a plurality of frequency points in real time by using a Goertzel algorithm and sliding a window point by point;
(3) and judging the signal according to the judging condition.
The invention is further configured to: in the step (1), two band-pass filters are used for screening out signals of two detection frequency bands.
The invention is further configured to: the two detection frequency ranges are respectively a low frequency range of 25MHz-75MHz and a high frequency range of 110MHz-150 MHz.
The invention is further configured to: in step (1), a sampling frequency fS is selected to be 160MHz, two a/D chips are used to simultaneously convert signals of two frequency bands, a low frequency band is directly sampled to obtain a digital signal of 25MHz to 75MHz, and a high frequency band is undersampled to obtain a digital signal of 10MHz to 50 MHz.
The invention is further configured to: in the step (1), the sampling time length is 256 us.
The invention is further configured to: in the step (2), the method comprises the following steps:
1) setting the number of channels to 20;
2) selecting 101 time domain signal sampling points with the calculation window width N being 101 to obtain a Fourier component frequency point amplitude, and calculating the Fourier component frequency point amplitude when moving one sampling point;
3) selecting a Fourier component frequency;
4) optimizing an iterative formula of a Goertzel algorithm;
5) and calculating Fourier components of a plurality of frequency points in real time by sliding a window point by point.
The invention is further configured to: the optimized Goertzel algorithm has the iteration formula as follows:
where N is the window width and X (N) is the sample value.
The invention is further configured to: in the step (3), the threshold and the duration of the threshold are used as the judgment conditions, and when the frequency point number meeting the judgment conditions reaches M and has the characteristic of high frequency first arrival, the signal is considered to be effective.
Advantageous effects
Compared with the known public technology, the technical scheme provided by the invention has the following beneficial effects:
the method calculates Fourier components of 20 frequency points in real time in two detection frequency bands of 25MHz-75MHz and 110MHz-150MHz, takes a threshold value, a duration time of the threshold value and the like as discrimination conditions, and when the frequency points meeting the discrimination conditions reach M (M is less than or equal to 20 and can be set), and has the characteristic of high frequency first arrival, the signals can be considered to be effective.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
In order to achieve the effects of the present invention, the present invention will be further described with reference to the following examples.
Example 1:
referring to fig. 1-3, the real-time processing method for point-by-point scanning of satellite-borne VHF transient signals based on Goertzel algorithm includes the following steps:
(1) and carrying out broadband sampling on the VHF transient signal to obtain an original time domain waveform of the detection frequency band.
In this embodiment, AD9642(14bits) is selected for sampling.
1) Selecting a sampling frequency band: the target VHF transient signal frequency is mainly distributed from dozens of kHz to hundreds of MHz, and signals larger than 20MHz can pass through the ionosphere due to the influence of ionosphere propagation effect. Under the influence of the ionized layer TEC, the group delay of electromagnetic waves with different frequencies in the ionized layer propagation process has obvious difference, and the delay is larger when the frequency is lower.
In order to reduce the influence of ground radio stations or other natural signals as much as possible, the detection frequency ranges are 25MHz-75MHz in low frequency range and 110MHz-150MHz in high frequency range, and before sampling, signals in the two frequency ranges are screened out by two band-pass filters.
2) Selecting a sampling frequency fS: for the probing Frequency band, the Nyquist Frequency (Nyquist Frequency) should be selected above the highest probing Frequency of the probing Frequency band or between the high Frequency band and the low Frequency band so as to cover the probing Frequency band; meanwhile, the nyquist frequency fN is to avoid selecting the edge of the detection frequency band as much as possible to reduce clutter interference, so the sampling frequency fS is selected to be 160 MHz. Two A/D chips are used for simultaneously converting signals of two frequency bands, low frequency bands are directly sampled to obtain digital signals of 25MHz-75MHz, and high frequency bands are undersampled to obtain digital signals of 10MH-50 MHz.
3) Selecting a sampling duration of an event: according to the ionized layer group delay characteristic, the TEC is 1 multiplied by 1017m-2The time delay difference of the electromagnetic waves of 30MHz and 31MHz is 1.62 multiplied by 103The time delay difference of the electromagnetic waves of ns, 150MHz and 151MHz is 7.9ns, and the total time delay difference of the electromagnetic waves of 30MHz and 150MHz is 1.43 multiplied by 104ns. ExaminationConsidering that the ionospheric concentration changes greatly, the sampling time of one event in the invention is 256us, namely 40960 sampling points.
(2) And calculating Fourier components of a plurality of frequency points in real time by using a Goertzel algorithm and sliding a window point by point.
In this embodiment, the method comprises the following steps:
1) setting the number of channels: in the past, the number of channels set by the satellite load is 5 channels (such as XX satellite number three in China) and 8 channels (such as FORTE satellite in the United states), and in consideration of the reliability of rapid identification of target signals, the number of the channels in the invention is 20, namely, Fourier components of 20 frequency points need to be calculated in real time when the acquired signals are processed.
2) Window width N: the window width N is that N Fourier transform frequency points are distributed in the acquisition frequency range (0-pi), and N times of iterative operation is carried out so as to obtain the Fourier transform frequency spectrum in the window width range. Fourier transform frequency point fk=kfsN, frequency point spacing fsN, wherein fSIs the sampling frequency of the system, N is the window width, fSand/N is the frequency resolution. The larger the N is selected, the higher the spectral resolution is, and the larger the required calculation amount is, so that N is related to calculation accumulated errors and system resources. For real-time signal property judgment, the isolation between selected frequency points and proper frequency point bandwidth are ensured, and frequency interference and detection reliability between the selected frequency points are avoided; but also to ensure computational accuracy and efficient use of resources. Through weighing, the calculation window width N is 101, namely 101 time domain signal sampling points are selected to obtain a Fourier component frequency point amplitude, and a Fourier component frequency point amplitude is calculated when every sampling point is moved.
3) Selecting a fourier component frequency: sampling frequency fS160MHz, so the Fourier component frequency resolution is 1.58416 MHz; the embodiment selects the frequency interval of the Fourier components to be 3.16832MHz, and 2 groups of Fourier components can be artificially arranged. In this embodiment, fourier component frequencies are selected as in table 1.
TABLE 1 Fourier component frequency Table
4) The iterative formula of the Goertzel algorithm is optimized:
iterating the formula according to a Goertzel recursion algorithm:
Vk(n)=2cos(2πk/N)Vk(n-1)-Vk(n-2) + X (n), when inputting a sampling value x (n), iteratively calculating multiplication and addition of real numbers which need to be carried out once, and calculating Vk(N), N calculations are required.
When the sampling rate is low, the iterative computation can be performed in the sampling interval, but when high-speed sampling is performed, the sampling interval is insufficient to complete the iterative computation of the size, so that the optimization of the iterative formula of the Goertzel algorithm is needed.
By derivation, the above iterative formula can be converted into:
where N is the window width and X (N) is the sample value, so that for V, where the window width is calculated to be N
k(N) calculation, can be by V
k(N-1)、V
k(N-2) and X (N) (0. ltoreq. N. ltoreq.n +1), directly calculating V
k(N +1), substituting N101 into the above-described converted equation gives:
in which w is the value of k for each given value of k
k、
Is a constant.
5) Calculating Fourier components of a plurality of frequency points in real time by sliding a window point by point:
in this example, the Zynq-7000 series chip XC7z045ffg 900-2 was used. To implement pipelined operations, 4 multipliers and 5 add (subtract) adders are required, where 2 multipliers and 2 adders are implemented with 2 multiply-add adders of DSP 48.
The computation time of the multiplier adopting a lookup table mode by the FPGA is 4 times of the computation speed (4clk), the computation time of the adding (subtracting) device adopting a fabric mode is 2clk, and the computation time of 2 multiplier-adder devices of the DSP48s is respectively multiplier 1clk, adder 3clk, multiplier 1clk andadder 2 clk.
When the computation speed fC of the FPGA is 2fS 320MHz, the pipeline method may be used to complete the formula
Iterative calculations (the calculation results are delayed by only a dozen clocks relative to the input data). The parallel processing advantage of the FPGA is fully utilized, the frequency domain analysis of multiple frequency points is completed in one sampling clock period, namely, each sampling clock simultaneously generates a group of 20-channel narrow-band multi-frequency-point frequency domain results, so that the real-time capture of transient signals is realized, wherein the point-by-point sliding window scanning detection results are shown in figure 1
(3) And judging the signal according to the judging condition.
The threshold value and the duration of the threshold value are used as discrimination conditions, when the frequency point number meeting the discrimination conditions reaches M (M is less than or equal to 20, the frequency point number can be set), and the signal is considered to be effective due to the characteristic that the frequency point number meets the discrimination conditions and has high frequency first arrival.
In this embodiment, as the data analysis received by the XX satellite No. three electromagnetic radiation detector, some specific areas often receive a single-frequency signal (such as radar, etc.), which has a long duration and is greatly different from the target signal. In order to filter such signals, energy peak value maintaining time judgment is added in an electromagnetic pulse receiver integrated processing unit, if the narrow-band signal continuously exceeds a threshold value and exceeds tPThe signal is not picked up, thus most of the non-target signals are filtered out. Meanwhile, the target source signal is a pulse signal, the frequency spectrum of the pulse signal can cover different detection frequency bands after ionosphere scattering, and the target source signal is judged through 20 narrow band channels and is within a certain time window tWOnly when M (less than or equal to 20, can be set) way meets the condition, the mining is carried outAnd the signal triggering judgment schematic diagram is shown in fig. 2.
The signal triggering process is as follows: writing a digitized threshold value into the latch, comparing the digitized acquisition signal with a set threshold value, and generating an over-threshold signal when the acquisition signal is greater than the set threshold value; when the signal passes through the peak value maintaining time judgment, the signal enters a 20M-selecting module for judgment, and the judgment is carried out in a time window tWAnd storing the event waveform when the over-threshold signal is larger than or equal to M within the time, and latching the event triggering time, wherein the signal triggering judgment process is shown in FIG. 3.
Example 2:
the method of the present invention is compared with the prior art, and the relevant data is shown in table 2.
TABLE 2 comparative data sheet
As can be seen from Table 2, the invention reduces the complexity of hardware circuits, enhances the anti-interference capability of the circuits, and has the advantages of real-time performance, high efficiency, adjustable parameters and the like; in addition, the acquisition of the original time domain signal waveform can be completed, and information is provided for ionospheric propagation group delay and waveform characteristic analysis.