CN115510787B - Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization - Google Patents
Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization Download PDFInfo
- Publication number
- CN115510787B CN115510787B CN202211164323.5A CN202211164323A CN115510787B CN 115510787 B CN115510787 B CN 115510787B CN 202211164323 A CN202211164323 A CN 202211164323A CN 115510787 B CN115510787 B CN 115510787B
- Authority
- CN
- China
- Prior art keywords
- window function
- function
- constraint
- optimization
- amplitude
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000005457 optimization Methods 0.000 title claims abstract description 56
- 238000000034 method Methods 0.000 title claims abstract description 33
- 238000013461 design Methods 0.000 title claims abstract description 24
- 238000001228 spectrum Methods 0.000 claims abstract description 26
- 238000010606 normalization Methods 0.000 claims abstract description 7
- 238000005070 sampling Methods 0.000 claims description 32
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 230000004044 response Effects 0.000 claims description 7
- 230000003595 spectral effect Effects 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000009467 reduction Effects 0.000 claims description 2
- 238000005562 fading Methods 0.000 claims 3
- 230000009466 transformation Effects 0.000 claims 1
- 230000001427 coherent effect Effects 0.000 description 24
- 230000006872 improvement Effects 0.000 description 9
- 238000010183 spectrum analysis Methods 0.000 description 8
- 230000000694 effects Effects 0.000 description 3
- 238000012795 verification Methods 0.000 description 2
- 230000002349 favourable effect Effects 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/30—Circuit design
- G06F30/32—Circuit design at the digital level
- G06F30/337—Design optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Landscapes
- Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Monitoring And Testing Of Transmission In General (AREA)
Abstract
The invention relates to a peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization, wherein the window function designed by the method is used for estimating the amplitude spectrum of a signal; the method comprises the following steps: according to the number of window function coefficients to be designed, determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak sidelobe level constraint condition, a rapid attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficients; and obtaining window function coefficients through optimization solution according to the determined optimization target and each constraint condition, thereby completing the design of the window function. According to the window function design method, the minimized maximum amplitude estimation error is used as an optimization target, and the peak sidelobe level and the attenuation rate are used as constraint conditions, so that a rapid attenuation window function with the minimum amplitude estimation error can be obtained on the premise of preferentially meeting the constraint of the desired sidelobe characteristic.
Description
Technical Field
The invention relates to the field of design of window functions required by signal spectrum analysis, in particular to a window function with good side lobe characteristics (high attenuation rate and low peak side lobe level) and low amplitude estimation error, which is used for amplitude spectrum estimation of an amplitude signal with a high dynamic range, and particularly relates to a peak side lobe constraint rapid attenuation window function design method based on multi-objective optimization.
Background
In real life, coherent sampling is difficult (i.e. the sampling time must be exactly an integer multiple of the signal period). Even though the sampling time is approximately an integer multiple of the signal period (only a few sample points apart), it is also an incoherent sampling, where spectral leakage and a fence effect still occur when performing spectral analysis with FFT processing, resulting in severe distortion of the spectral estimate, especially for signals with high dynamic range amplitudes, i.e. signals with amplitude ranges spanning 20 orders of magnitude or even more. Spectral leakage will cause frequency components with small amplitudes to be masked by leakage of frequency components with large amplitudes, resulting in failure to recover them. For better spectral analysis, more accurate estimation of the signal amplitude, a suitable window function needs to be chosen to reduce the effects of leakage.
Each window function has its own characteristics, and different window function types are suitable for different application scenarios. By selecting a window function with good side lobe characteristics, the interference of adjacent frequency components and the influence of long-distance leakage can be reduced. The dynamic range of the analyzed signal is also a factor to consider when making window function selections. For voltage and current signals, typically the amplitude ranges from only a few volts to a few hundred volts, i.e. spans around 40 dB; for GEO600 tasks, its dynamic range exceeds 120dB; whereas for inter-satellite distance signals in deep space exploration tasks the dynamic range typically spans more than 20 orders of magnitude, i.e. more than 400dB. In spectral analysis, the high dynamic range makes frequency components with large amplitudes highly affected by frequency components with small amplitudes. Thus, for signals with high dynamic range at incoherent sampling, a high decay rate window function is very important for achieving a more accurate spectral analysis.
The sidelobe characteristics of the current window function are insufficient to meet the requirements of amplitude spectrum estimation at that time, and therefore, it is necessary to design a window function having a fast attenuation rate, while having a low peak sidelobe level and a low amplitude estimation error. In summary, currently, there is no suitable window function for the problem of amplitude spectrum estimation for signals with high dynamic range amplitudes.
Disclosure of Invention
The invention aims to solve the problem of insufficient window function capability in the amplitude spectrum estimation of an amplitude signal with a high dynamic range, and provides a peak sidelobe constraint rapid attenuation window function design method based on constraint multi-objective optimization, which shows the properties of the obtained window function and the application of the window function in the amplitude spectrum estimation of the amplitude signal with the high dynamic range.
In order to achieve the above purpose, the present invention is realized by the following technical scheme.
The invention provides a peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization, wherein the window function designed by the method is used for estimating the amplitude spectrum of a signal; the method comprises the following steps:
according to the number of window function coefficients to be designed, determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak sidelobe level constraint condition, a rapid attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficients;
and obtaining window function coefficients through optimization solution according to the determined optimization target and each constraint condition, thereby completing the design of the window function.
As one of the improvements of the above technical solutions, the window function designed by the method includes a combined cosine window; the discrete form of the combined cosine window function g (n) is:
wherein M+1 is the number of window function coefficients, N is the number of samples of the signal sample, a m For the m+1th coefficient, m=0, 1,2, …, M, n represents the nth sample.
As one of the improvements of the above technical solution, the number of the fast decay constraints requires: m or less and 3 or more; the fast decay constraint conditions all satisfy the following formula:
under the fast decay constraint, the window function is used for
Attenuation, h is normalized frequency.
As one of the improvements of the above technical solution, the fast attenuation constraint condition is derived by discrete fourier transform of a window function, and specifically includes:
the window function g (n) is converted into by the euler equation:
wherein the variables are
Converting the window function g (n) after conversion into the following discrete Fourier transform
The expression is:
wherein the function is
Function of
Further simplifying and obtaining:
will be 1/(h) 2 -m 2 ) And (5) expanding power series to obtain:
the attenuation rate of the side lobe is determined according to the first non-zero term of the r-series.
As one of the improvements of the above technical solution, the expression for minimizing the maximum amplitude estimation error is:
where k represents the kth frequency component, the variable
h k For the position of the maximum spectral line of the kth frequency component,/->
[·]Representing an integer part function; sinc (·) represents a sine function, abs (·) represents an absolute value function, max (·) represents a maximum value function, and min (·) represents a minimum value function; a, a 0 Representing the first coefficient of the window function.
As one of the improvements of the above technical solution, the optimization objective of minimizing the maximum amplitude estimation error is obtained according to the discrete fourier transform of the signal to be estimated, and specifically includes the following steps:
setting a signal x to be estimated 1 (n) consists of a plurality of frequency components, expressed as:
wherein ,x1k (n) is the kth frequency component of the signal to be estimated, expressed as:
wherein ,
represents the amplitude of the kth frequency component, kf represents the frequency of the kth frequency component, f is the fundamental frequency, i.e. the frequency of the first frequency component, Δt represents the sampling interval, α 1k Representing the phase of the kth frequency component;
calculating a signal x 1 (n) windowed discrete Fourier transform
The calculation formula is as follows:
wherein the function is
h is the normalized frequency; />
Is the spectrum of the kth harmonic component, is the sum of the two spectra, and is expressed as:
wherein ,
discrete fourier transform of window function g (n) converted by euler formula;
focusing only on the first spectrum, namely:
wherein ,
for phasor, add>
Δf is the frequency resolution, +.>
wherein fs Is the sampling frequency; />
When (when)
When the number is a non-integer, where h=h k Maximum spectral line->
The expression is:
wherein ,
subtracting the true amplitude to obtain an amplitude estimation error as follows:
wherein ,|δk |≤0.5
Further reduction to using sinc (·) function
The optimization objective to minimize the maximum amplitude estimation error is:
as one of the improvements of the above technical solution, the expression of the normalization constraint condition is:
as one of the improvements of the above-described technical means, the peak sidelobe level constraint is observed by an amplitude response;
the amplitude response is defined as:
the peak sidelobe level constraint is defined as:
where h > M+1, sidelobe is a predefined peak sidelobe level constraint value.
As one of the improvements of the above technical solution, the constraint conditions of the upper and lower limits of the window function coefficients are: for each coefficient, the requirement is:
0≤a m ≤1,0≤m≤M
wherein M+1 is a coefficient number, a m Representing the m+1th coefficient.
As one of the improvement of the technical scheme, the method uses the fmincon function to carry out optimization solution to obtain window function coefficients, and specifically comprises the following steps:
s1, setting the number of window function coefficients, an initial value of the window function coefficients and an initial value of a desirobe;
s2, setting parameters required by an fmincon function, wherein an Algorithm Algorithm selects a sqp Algorithm;
s3, inputting the set window function coefficient initial value, the ideal value and parameters required by the fmincon function, the optimization objective function and each constraint condition into the fmincon function, and carrying out optimization solution;
s4, judging a return value exitflag of the fmincon function: if the exitflag= -2, the feasible solution cannot be found, at the moment, the value sidelobe of the peak sidelobe level constraint is automatically adjusted, the sidelobe+1 is assigned to the sidelobe, the feasible region of the solution is gradually enlarged, the step S3 is carried out, and the optimization solution is carried out again; if exitflag is not equal to-2, step S5 is carried out;
s5, further judging the exitflag: if exitflag is not equal to 1 and exitflag is not equal to 2, the iteration times are beyond the setting or the function is not converged, and the parameters need to be reset at the moment, namely, the step S2 is skipped; if exitflag=1 or exitflag=2, it means that a feasible solution is found, and the window function coefficients obtained by optimization are output.
Compared with the prior art, the invention has the advantages that:
1. the invention is a more general window function design method, can design 5 coefficients and above window functions with rapid attenuation characteristics, and meets the requirements of different application scenes on a plurality of window functions with different coefficients;
2. the window function designed by the invention has wide applicability and is also suitable for spectrum analysis of other signals with high dynamic range amplitude;
3. according to the window function design method, the minimized maximum amplitude estimation error is used as an optimization target, and the peak sidelobe level and the attenuation rate are used as constraint conditions, so that a rapid attenuation window function with the minimum amplitude estimation error can be obtained on the premise of preferentially meeting the constraint of the desired sidelobe characteristic. The method is more suitable for spectrum analysis of signals sampled by approximate integer times of the period;
4. by estimating the amplitude spectrum of the inter-satellite distance signal, it can be seen that for signals having amplitudes spanning 20 orders of magnitude, a good amplitude estimation result can be obtained with a maximum relative estimation error within 0.44%.
Drawings
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2 is an amplitude response of a signal estimated using the method of the present invention;
FIG. 3 is a high dynamic range signal amplitude spectrum estimate of whether or not coherent sampling is performed, wherein FIG. 3 (a) is a graph comparing amplitude estimates of coherent and incoherent sampling without windowing, and FIG. 3 (b) is a graph comparing amplitude estimates of 6 fast decay windows of coherent and incoherent sampling with updating;
FIG. 4 is a graph of the difference in amplitude estimation of a high dynamic range signal for coherent sampling, wherein FIG. 4 (a) is a graph of the difference in amplitude estimation for coherent and incoherent sampling, and FIG. 4 (b) is a graph of the difference in amplitude estimation for 6 fast decay windows for coherent and incoherent sampling;
fig. 5 is a graph of the percentage of difference between the amplitude estimates of the high dynamic range signal of whether the coherent samples are relatively different, wherein fig. 5 (a) is a graph of the percentage of difference between the amplitude estimates of the coherent samples and the non-coherent samples, and fig. 5 (b) is a graph of the percentage of difference between the amplitude estimates of the 6 fast decay windows of the coherent samples and the non-coherent samples.
Detailed Description
The application provides a peak sidelobe constraint rapid attenuation window function design method based on epsilon-constraint multi-objective optimization, which solves the problem of amplitude spectrum estimation of signals with high dynamic range amplitude by designing a window function with rapid attenuation rate, low peak sidelobe level and low amplitude estimation error, and realizes the estimation of the amplitudes of signals with amplitude ranges exceeding more than 20 orders of magnitude. In addition, verification is performed by the simulated inter-satellite distance signal to evaluate the validity of the obtained window function.
The design of window functions involves a multi-objective optimization problem that solves in a number of ways. Wherein the epsilon-constrained multi-objective optimization model converts the original multi-objective optimization problem into a single-objective optimization problem by converting one of the objectives into the form of constraint conditions, and then the problem can be solved by a single-objective optimization method.
The window function of interest in the present invention is a combined cosine window in discrete form
Wherein the coefficient term number is M+1, and N is the sample number of the sample.
It is first necessary to determine the optimization objectives and constraints of the window function. Since a window function having good side lobe characteristics is focused, by setting both the peak side lobe level and the attenuation rate as constraints, it is possible to more conveniently obtain the desired side lobe characteristics and minimize the amplitude estimation error as much as possible in the feasible region as an optimization target.
Optimization target: minimizing maximum amplitude estimation error
Constraint:
(1) fast decay constraint
(2) Peak sidelobe level constraint
(3) Normalized constraint
The following gives a formulation modeling of optimization objectives and constraints.
The fast decay constraint can be derived by the discrete fourier transform DFT of the window function, and equation (1) can be converted into by the euler equation
Consider the following DFT transform
DFT transforming window function g (n)
wherein 
Further simplify the approximation to obtain
Will be 1/(h) 2 -m 2 ) Expanding the power series, substituting the power series to obtain
The decay rate of the side lobe depends on the first non-zero term of the r-series, so when
When the first non-zero term of the series is obtained when r=0, then the window function decays at 1/h. When (when)
When the first non-zero term is obtained at least when r=1, then the window function is at least 1/h 3 Attenuation. At the same time, when
When the first non-zero term is obtained at least when r=2, then the window function is at least 1/h 5 Attenuation. At the same time, when
When the first non-zero term is obtained at least when r=3, then the window function is at least 1/h 7 Attenuation. At the same time, when
When the first non-zero term is obtained at least when r=4, the window function is at least 1/h 9 Attenuation. At the same time, when
When the first non-zero term is obtained at least when r=5, the window function is at least 1/h 11 Attenuation. And so on, it can continue to deduce
At least 1/h of window function 13 Attenuation, at which time a very high attenuation rate has been reached. The constraint of a better decay rate can be continued to be derived if desired.
The peak sidelobe level of the window function can be observed by an amplitude response, which is defined as
The peak sidelobe level constraint can thus be defined as
Where the sidelobe is a predefined hypothesis, different sidelobes result in different feasible regions.
The normalized constraint is defined as
The amplitude estimation error is obtained according to DFT after signal windowing, and the signal x to be estimated is assumed 1 (n) is composed of a plurality of frequency components
wherein ,
represents the amplitude of the kth frequency component, k f represents the frequency of the kth frequency component, f is the fundamental frequency, i.e. the frequency of the first frequency component, Δt represents the sampling interval, N represents the nth sample, N represents the total number of samples sampled, α 1k Representing the phase of the kth frequency component;
wherein order
Is the conjugate of (I)>
Δf is the frequency resolution, +.>
wherein fs Is the sampling frequency. Can obtain
Calculating a signal x 1 (n) windowed DFT is expressed as
Where h is the normalized frequency;
is the signal x 1 (n) form of windowed DFT;
wherein ,
it can be seen that the spectrum of the kth harmonic component is the sum of the two spectra, we focus on only the first spectrum, i.e
When (when)
If the spectrum is a non-integer, leakage occurs, and if the influence of interference superposition is ignored, the maximum spectrum line is represented by h=h k Let->
wherein />
Is an integer part of>
Can be obtained according to the formula (5)
Then for the kth frequency component the estimated amplitude is
Whereas the true amplitude is +.>
The expression that the amplitude estimation error can be obtained by subtracting the two is that
And then deriving the amplitude estimation error for the kth frequency component as
The use of sinc (delta) functions can be further simplified into
The optimization objective of minimizing the maximum amplitude estimation error can thus be defined as
The technical scheme provided by the invention is further described below with reference to the accompanying drawings and the embodiments.
The window function design flow according to the embodiment of the present invention is shown in fig. 1, and the flow according to the present invention is described below with reference to fig. 1.
1. First, the number of window function coefficients to be optimized is determined and is denoted as the (M+1) term. In order to obtain a fast decaying window function, more than 5 coefficients are usually selected to apply more fast decaying constraints, but the upper limit is not limited, but is usually less than 10, so that the fast decaying window function is more favorable for obtaining fewer signal periods and better spectrum analysis can be performed.
2. According to the number of coefficients determined in step 1, the initial value of the coefficient of the optimization time window function is set to be a random number, that is, a0=rand (m+1, 1).
3. The initial value sidelobe of the peak sidelobe level constraint is predefined, and the initial value is negative, usually-200 to-80 dB, but is not limited to the interval, and the initial value is determined according to the number (M+1) of the coefficients selected in the step 1, and can be roughly set according to the relation of-20 x (M+1) or can be set by self.
4. Based on the number of coefficients determined in step 1, an expression for optimizing the objective-minimizing the maximum amplitude estimation error, equation (27), is determined.
5. And (3) determining a normalization constraint condition, namely a formula (16), according to the number of coefficients determined in the step (1). This condition is a constraint necessary to design the window function.
6. And (3) determining the required rapid attenuation constraint conditions according to the number of the coefficients determined in the step (1). The number of fast decay constraints is required to be equal to or less than M, whereas in order to obtain a window function with fast decay characteristics, the required fast decay constraints generally need to be equal to or greater than 3, i.e. equations (8) - (10), i.e. at least guarantee a window function of 1/h 7 Attenuation.
7. Determining upper and lower bound constraint of coefficients, wherein each coefficient is 0.ltoreq.a m ≤1,0≤m≤M。
8. The nonlinear optimization function fmincon in matlab is utilized to carry out optimization solution of the window function, so that parameters of the window function need to be set. The Algorithm selects sqp, namely a sequence quadratic programming Algorithm, and the rest parameters can be adjusted by themselves or default values can be adopted.
9. After the parameter setting is completed, the optimization solution is started.
10. If the exit condition exitflag= -2 is optimized, the feasible solution cannot be found, at this time, the peak sidelobe level constraint sidelobe value is automatically adjusted, the sidelobe=sidelobe+1 is caused to gradually expand the feasible region of the solution, and then the step 9 is returned to for the optimization solution again. If exitflag is not equal to-2, proceed to step 11.
11. If the exit condition exitflag is not equal to 1 and exitflag is not equal to 2, the iteration number exceeds the setting or the function is not converged, and the parameters need to be reset at this time, namely, the step 8 is skipped.
12. If the exit condition exitflag=1 or exitflag=2 is optimized, the feasible solution is found, and the window function coefficient obtained by optimization is output.
13. The characteristics of the window function obtained through calculation comprise peak sidelobe level and maximum amplitude estimation error, and the attenuation rate of the peak sidelobe level and the maximum amplitude estimation error is fitted to judge whether the window function obtained through optimization meets the rapid attenuation constraint condition.
Taking 6 coefficients as an example, according to the above procedure, a new window function with rapid decay characteristic is obtained by optimization.
For the problem of amplitude spectrum estimation of a signal with sampling time approximately an integer multiple of the signal period (approximately 8 periods) and high dynamic range amplitude, 6 coefficients are selected, an initial value side lobe= -120dB of peak sidelobe level constraint is predefined, and an optimization target-minimum maximum amplitude estimation error is determined to be
The normalization constraint is that
To obtain a fast decaying window function, 4 decay conditions are chosen, i.e
Determining upper and lower bound limits of coefficients, for each coefficient
0≤a m ≤1,0≤m≤5
Setting fmincon function parameters, wherein an Algorithm Algorithm selects sqp, namely a sequence quadratic programming Algorithm; the tolerance TolCon for constraint violations is set to 1e-15; the first order optimal termination tolerance TolFun is set to 1e-10; the termination tolerance of x is set to 1e-14; the maximum number of allowed function calculations MaxFunEvals is set to 200000; the maximum number of iterations allowed MaxIter is set to 1500. The remaining parameters are default values.
Obtaining a 6-term coefficient window function with rapid attenuation characteristic through optimizing and solving, wherein the coefficient is
The new window function has peak sidelobe level of-95 dB, maximum amplitude estimation error of about 4.78%, and rapid attenuation characteristic of 1/h 9 Attenuation, the amplitude response of which is shown in figure 2.
Verification using a simulated signal with a high dynamic range amplitude does not fall within the scope of the technical solution of the present application, but it is necessary to explain the validity of a new window obtained by the window function design method of the present application because it needs to be verified.
For signals with high dynamic range amplitude, the space-to-star signal of the task "taiji" is detected using spatial gravitational waves considering only the gravitational field of the sun. The inter-satellite signal at this time is a periodic signal, but the period is non-integer, but the sampling time is just an integer multiple of the signal period by setting the sampling interval to be a decimal during simulation, namely, coherent sampling is performed, so that the amplitude of the inter-satellite signal can be accurately estimated at this time; in reality, the sampling interval is usually an integer, so that the sampling time is not exactly an integer multiple of the signal period, which causes spectrum leakage, and windowing is needed to reduce the influence of the leakage. Therefore, the new 6 fast decay windows obtained by optimizing the method of the invention are used for estimating the amplitude of the signal when leakage exists and comparing the amplitude with the result when coherent sampling exists.
The signal amplitude can be estimated very well without windowing during coherent sampling, and can be used as a true value for comparison; estimating the amplitude of the signal by using a new 6-item rapid attenuation window obtained by optimization during incoherent sampling; the non-windowed amplitude estimation results are also shown for comparison during non-coherent sampling.
The sampling time length is selected to be 8 periods, and the fundamental frequency of the inter-satellite distance signal
f 0 ＝3.168748934534878517772678686201385×10 -8 Hz
Comparing the first 16 multiplied amplitude estimates, i.e. f k ＝k·f 0 K=1, 2,... With an amplitude dynamic range of about 8 x 10 6 m～4×10 -14 m, i.e. spanning 20 orders of magnitude.
The obtained amplitude estimation results and differences are shown in fig. 3 and fig. 4, wherein fig. 3 (a) is an amplitude estimation comparison graph of coherent sample non-windowing and incoherent sample non-windowing, fig. 3 (b) is an amplitude estimation comparison graph of coherent sample non-windowing and incoherent sample non-windowing 6 fast attenuation windows, fig. 4 (a) is an amplitude estimation difference graph of coherent sample non-windowing and incoherent sample non-windowing, and fig. 4 (b) is an amplitude estimation difference graph of coherent sample non-windowing and incoherent sample non-windowing 6 fast attenuation windows. As can be seen from fig. 3 (a) and fig. 4 (a), when the incoherent sampling is not windowed, the amplitude estimation result is very different from the amplitude estimation result of the coherent sampling, and as can be seen from fig. 3 (b) and fig. 4 (b), when the incoherent sampling is added with a new 6-term rapid attenuation window, the amplitude estimation result is very different from the amplitude estimation result of the coherent sampling, so that the influence of leakage is effectively reduced.
The percentage of amplitude estimation error versus true value for both is shown in fig. 5, where fig. 5 (a) is a graph of the percentage of amplitude estimation difference for coherent sample un-windowed and incoherent sample un-windowed, and fig. 5 (b) is a graph of the percentage of amplitude estimation difference for coherent sample un-windowed and incoherent sample updated 6 fast decay windows; it can be seen that the new 6 fast decay windows obtained by using the method of the present application, when non-coherently sampled, have a relative estimation error within 0.44% for the amplitude of the signal with high dynamic range (amplitude range spans about 20 orders of magnitude), well reducing the effect of leakage.
As can be seen from the above detailed description of the present invention, the window function design method of the present invention uses the minimized maximum amplitude estimation error as the optimization target, and uses the peak sidelobe level and the attenuation rate as the constraint conditions, so that the rapid attenuation window function with the minimum amplitude estimation error can be obtained on the premise of preferentially satisfying the desired sidelobe characteristic constraint. The window function designed by the method is used for estimating the amplitude spectrum of the signal, especially for the signal with high dynamic range, and for the estimation of the amplitude which spans over 20 orders of magnitude, the relative error between the amplitude spectrum obtained by estimation and the true value can be realized to be within 0.44 percent.
Finally, it should be noted that the above embodiments are only for illustrating the technical solution of the present invention and are not limiting. Although the present invention has been described in detail with reference to the embodiments, it should be understood by those skilled in the art that modifications and equivalents may be made thereto without departing from the spirit and scope of the present invention, which is intended to be covered by the appended claims.
Claims (3)
1. A peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization is provided, and the window function designed by the method is used for estimating the amplitude spectrum of a signal; the method comprises the following steps:
according to the number of window function coefficients to be designed, determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak sidelobe level constraint condition, a rapid attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficients;
obtaining window function coefficients through optimization solution according to the determined optimization targets and each constraint condition, thereby completing the design of window functions;
the window function designed by the method comprises a combined cosine window; the discrete form of the combined cosine window function g (n) is:
wherein M+1 is the number of window function coefficients, N is the number of samples of the signal sample, a m For the m+1th coefficient, m=0, 1,2, …, M, n represents the nth sample;
the number of the fast decay constraint conditions requires: m or less and 3 or more; the fast decay constraint conditions all satisfy the following formula:
under the fast decay constraint, the window function is used for
Attenuation, h is normalized frequency;
the expression for minimizing the maximum amplitude estimation error is:
wherein k represents the kth, variable
h k Kf represents the frequency of the kth frequency component, Δf is the frequency resolution, for the position of the maximum spectral line of the kth frequency component>
wherein fs Is the sampling frequency, +.>
[·]Representing an integer part function; sinc (·) represents a sine function, abs (·) represents an absolute value function, max (·) represents a maximum value function, and min (·) represents a minimum value function; a, a 0 A first coefficient representing a window function;
the expression of the normalization constraint condition is:
the peak sidelobe level constraint is observed through an amplitude response;
the amplitude response is defined as:
the peak sidelobe level constraint is defined as:
wherein, h is more than M+1, and sidelobe is the value of the peak sidelobe level constraint defined in advance;
the constraint conditions of the upper limit and the lower limit of the window function coefficients are as follows: for each coefficient, the requirement is:
0≤a m ≤1,0≤m≤M
wherein M+1 is a coefficient number, a m Represents the m+1th coefficient;
the method uses fmincon function to carry out optimization solution to obtain window function coefficients, and specifically comprises the following steps:
s1, setting the number of window function coefficients, an initial value of the window function coefficients and an initial value of a desirobe;
s2, setting parameters required by an fmincon function, wherein an Algorithm Algorithm selects a sqp Algorithm;
s3, inputting the set window function coefficient initial value, the ideal value and parameters required by the fmincon function, the optimization objective function and each constraint condition into the fmincon function, and carrying out optimization solution;
s4, judging a return value exitflag of the fmincon function: if the exitflag= -2, the feasible solution cannot be found, at the moment, the value sidelobe of the peak sidelobe level constraint is automatically adjusted, the sidelobe+1 is assigned to the sidelobe, the feasible region of the solution is gradually enlarged, the step S3 is carried out, and the optimization solution is carried out again; if exitflag is not equal to-2, step S5 is carried out;
s5, further judging the exitflag: if exitflag is not equal to 1 and exitflag is not equal to 2, the iteration times are beyond the setting or the function is not converged, and the parameters need to be reset at the moment, namely, the step S2 is skipped; if exitflag=1 or exitflag=2, it means that a feasible solution is found, and the window function coefficients obtained by optimization are output.
2. The peak sidelobe constrained fast fading window function design method based on multi-objective optimization as set forth in claim 1, wherein the fast fading constraint is derived by discrete fourier transform of a window function, and specifically comprises:
the window function g (n) is converted into by the euler equation:
wherein the variables are
Will be converted intoIs transformed by discrete fourier transformation into
The expression is:
wherein the function is
Function of
Further simplifying and obtaining:
will be 1/(h) 2 -m 2 ) And (5) expanding power series to obtain:
the attenuation rate of the side lobe is determined according to the first non-zero term of the r-series.
3. The peak sidelobe constrained fast fading window function design method based on multi-objective optimization according to claim 2, wherein the minimized maximum amplitude estimation error optimization objective is obtained according to a windowed discrete fourier transform of a signal to be estimated, and specifically comprises the following steps:
setting a signal x to be estimated 1 (n) consists of a plurality of frequency components, expressed as:
wherein ,x1k (n) is the kth frequency component of the signal to be estimated, expressed as:
wherein ,
representing the amplitude of the kth frequency component, f being the fundamental frequency, i.e. the frequency of the first frequency component, Δt representing the sampling interval, α 1k Representing the phase of the kth frequency component;
calculating a signal x 1 (n) windowed discrete Fourier transform
The calculation formula is as follows:
wherein the function is
h is the normalized frequency; />
Is the spectrum of the kth harmonic component, is the sum of the two spectra, and is expressed as:
wherein ,
is meridian passageDiscrete fourier transform of window function g (n) after euler formula conversion;
focusing only on the first spectrum, namely:
wherein ,
for phasor, add>
When (when)
When the number is a non-integer, where h=h k Maximum spectral line->
The expression is:
wherein ,
subtracting the true amplitude to obtain an amplitude estimation error as follows:
wherein ,|δk |≤0.5
Further reduction to using sinc (·) function
The optimization objective to minimize the maximum amplitude estimation error is:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211164323.5A CN115510787B (en) | 2022-09-23 | 2022-09-23 | Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211164323.5A CN115510787B (en) | 2022-09-23 | 2022-09-23 | Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115510787A CN115510787A (en) | 2022-12-23 |
| CN115510787B true CN115510787B (en) | 2023-04-21 |
Family
ID=84505887
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211164323.5A Active CN115510787B (en) | 2022-09-23 | 2022-09-23 | Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115510787B (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109871575A (en) * | 2018-12-29 | 2019-06-11 | 陕西海泰电子有限责任公司 | A kind of design method of the electromagnetic interference receiver window function based on time domain FFT |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101692747A (en) * | 2009-09-29 | 2010-04-07 | 天津理工大学 | Design of coexisting filter on the basis of DCS1800 system and TD-SCDMA system interference resistance |
| CN109598093B (en) * | 2018-12-29 | 2020-12-04 | 北京化工大学 | Fitting window function-based seismic vector wave field numerical simulation method and system |
| US10915677B1 (en) * | 2020-07-22 | 2021-02-09 | North China Electric Power University | General design method for phasor estimation in different applications |
| CN112634800A (en) * | 2020-12-22 | 2021-04-09 | 北方液晶工程研究开发中心 | Method and system for rapidly and automatically testing refresh frequency of light-emitting diode display screen |
| CN114116370B (en) * | 2021-08-31 | 2024-04-26 | 西南电子技术研究所（中国电子科技集团公司第十研究所） | Complex electronic system operation health state monitoring point optimization method |
| CN114966201B (en) * | 2022-05-05 | 2023-07-04 | 威胜能源技术股份有限公司 | Harmonic voltage and current measurement and characteristic analysis method in typical scene |
-
2022
- 2022-09-23 CN CN202211164323.5A patent/CN115510787B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109871575A (en) * | 2018-12-29 | 2019-06-11 | 陕西海泰电子有限责任公司 | A kind of design method of the electromagnetic interference receiver window function based on time domain FFT |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115510787A (en) | 2022-12-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106597408B (en) | High-order PPS signal parameter estimation method based on time-frequency analysis and instantaneous frequency curve fitting | |
| CN108037361B (en) | High-precision harmonic parameter estimation method based on sliding window DFT | |
| CN108470089B (en) | Complex signal time delay estimation method based on least square sample fitting | |
| Vanbeylen et al. | Blind maximum-likelihood identification of Wiener systems | |
| JP2001319178A (en) | Excitation signal used for extracting nonlinear black box behavior model and radial basis function method | |
| CN111222088B (en) | Improved method for estimating weighted power harmonic amplitude of flat-top self-convolution window | |
| O'Shea | A high-resolution spectral analysis algorithm for power-system disturbance monitoring | |
| Tamim et al. | Techniques for optimization in time delay estimation from cross correlation function | |
| CN114785379B (en) | Method and system for estimating parameters of underwater sound JANUS signals | |
| CN105071858A (en) | Dispersion estimation method in fiber communication system | |
| CN111159888B (en) | Covariance matrix sparse iteration time delay estimation method based on cross-correlation function | |
| US6892155B2 (en) | Method for the rapid estimation of figures of merit for multiple devices based on nonlinear modeling | |
| Shevgunov et al. | Averaged absolute spectral correlation density estimator | |
| CN114842867A (en) | DFT-based audio sinusoidal signal frequency estimation method and system | |
| CN115510787B (en) | Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization | |
| CN112910533B (en) | Broadband signal array system with parallel structure | |
| US10393865B2 (en) | Phase retrieval algorithm for generation of constant time envelope with prescribed fourier transform magnitude signal | |
| Randall et al. | Updating modal models from response measurements | |
| CN111917676B (en) | Linear frequency modulation interference cancellation method | |
| Xianmin | A new method with high confidence for validation of computer simulation models of flight systems | |
| CN112883787B (en) | Short sample low-frequency sinusoidal signal parameter estimation method based on spectrum matching | |
| Tamim et al. | Hilbert transform of FFT pruned cross correlation function for optimization in time delay estimation | |
| Wang et al. | Accurate and computationally efficient interpolation-based method for two-dimensional harmonic retrieval | |
| CN112557751A (en) | Harmonic parameter estimation method based on DFT iteration method | |
| CN103441975A (en) | Two-phase coding signal parameter estimation method based on power spectrum |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |