CN115510787A - Peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization - Google Patents

Peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization Download PDF

Info

Publication number
CN115510787A
CN115510787A CN202211164323.5A CN202211164323A CN115510787A CN 115510787 A CN115510787 A CN 115510787A CN 202211164323 A CN202211164323 A CN 202211164323A CN 115510787 A CN115510787 A CN 115510787A
Authority
CN
China
Prior art keywords
window function
constraint
function
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.)
Granted
Application number
CN202211164323.5A
Other languages
Chinese (zh)
Other versions
CN115510787B (en
Inventor
韩晓晴
唐文林
彭晓东
杨震
马晓珊
强丽娥
张玉珠
高辰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Space Science Center of CAS
Original Assignee
National Space Science Center of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by National Space Science Center of CAS filed Critical National Space Science Center of CAS
Priority to CN202211164323.5A priority Critical patent/CN115510787B/en
Publication of CN115510787A publication Critical patent/CN115510787A/en
Application granted granted Critical
Publication of CN115510787B publication Critical patent/CN115510787B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

Abstract

The invention relates to a peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization, wherein a window function designed by the method is used for estimating an amplitude spectrum of a signal; the method comprises the following steps: determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak side lobe level constraint condition, a fast attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficient according to the number of the window function coefficients to be designed; and obtaining a window function coefficient 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 minimum maximum amplitude estimation error is used as an optimization target, the peak side lobe level and the attenuation rate are used as constraint conditions, and the fast attenuation window function with the minimum amplitude estimation error can be obtained on the premise that the desired side lobe characteristic constraint is met preferentially.

Description

Peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization
Technical Field
The invention relates to the field of design of window functions required by signal spectrum analysis, in particular to a method for designing a window function with good side lobe characteristics (high attenuation rate and low peak side lobe level) and low amplitude estimation error for amplitude spectrum estimation of amplitude signals 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
When a signal is sampled in real life, coherent sampling is difficult to achieve (namely, the sampling time must be exactly an integral multiple of the signal period). Even if the sampling time is approximately an integer multiple of the signal period (only a few samples apart), it is also incoherent, and in this case, when the FFT process is used for spectrum analysis, the spectrum leakage and the fence effect still occur, which causes serious distortion of the spectrum estimation, especially for signals with high dynamic range amplitude, i.e. signals with amplitude range spanning 20 orders or more. The spectral leakage will cause frequency components with small amplitudes to be masked by the leakage of frequency components with large amplitudes, resulting in their being unrecoverable. For better spectral analysis and more accurate estimation of the amplitude of the signal, it is necessary to select a suitable window function to reduce the effect 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, interference of adjacent frequency components and the influence of remote leakage can be reduced. The dynamic range of the analyzed signal is also a factor to be considered when making the window function selection. For voltage and current signals, the amplitude typically ranges only from a few volts to hundreds of volts, i.e., spanning around 40 dB; for the GEO600 task, its dynamic range exceeds 120dB; for the inter-satellite distance signal in deep space exploration tasks, the dynamic range usually spans more than 20 magnitude, i.e. more than 400dB. In the spectral analysis, the high dynamic range causes frequency components having large amplitudes to greatly affect frequency components having small amplitudes. Therefore, for signals with high dynamic range under incoherent sampling, a window function with high attenuation rate is very important for realizing more accurate spectrum analysis.
The side lobe characteristic of the current window function is not enough to meet the requirement of amplitude spectrum estimation at the moment, so that the window function with a fast attenuation rate, a lower peak side lobe level and a lower amplitude estimation error is required to be designed. In summary, no suitable window function is currently available for the problem of amplitude spectrum estimation of signals with high dynamic range amplitudes.
Disclosure of Invention
The invention aims to overcome the problem of insufficient capacity of a window function in amplitude spectrum estimation of an amplitude signal with a high dynamic range at present, provides a peak sidelobe constraint fast attenuation window function design method based on constraint multi-objective optimization, and shows the properties of the obtained window function and the application of the window function in amplitude spectrum estimation of the amplitude signal with the high dynamic range.
In order to achieve the above purpose, the invention is realized by the following technical scheme.
The invention provides a peak sidelobe constraint fast 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:
determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak side lobe level constraint condition, a fast attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficient according to the number of the window function coefficients to be designed;
and obtaining a window function coefficient through optimization solution according to the determined optimization target and each constraint condition, thereby completing the design of the window function.
As one improvement of the above technical solution, the window function designed by the method includes a combination cosine window; the discrete form of the combined cosine window function g (n) is:
Figure BDA0003861524950000021
wherein M +1 is the number of window function coefficients, N is the number of samples of signal samples, a m M =0,1,2, \ 8230for the M +1 coefficient, M, n denotes the nth sample.
As one improvement of the above technical solution, the number of the fast fading constraints requires: less than or equal to M and more than or equal to 3; the fast decay constraints all satisfy the following equation:
Figure BDA0003861524950000022
under fast fading constraints, a window function of
Figure BDA0003861524950000023
Attenuation, h is the normalized frequency.
As an improvement of the above technical solution, the fast attenuation constraint is derived by discrete fourier transform of a window function, and specifically includes:
the window function g (n) is converted by the euler formula into:
Figure BDA0003861524950000024
wherein, variable
Figure BDA0003861524950000025
Transforming the transformed window function g (n) by discrete Fourier transform
Figure BDA0003861524950000031
The expression is as follows:
Figure BDA0003861524950000032
wherein the function
Figure BDA0003861524950000033
Function(s)
Figure BDA0003861524950000034
Further simplification results in:
Figure BDA0003861524950000035
will be 1/(h) 2 -m 2 ) The power series expansion yields:
Figure BDA0003861524950000036
and determining the decay rate of the side lobe according to the first non-zero term of the r series.
As an improvement of the above technical solution, the expression for minimizing the maximum amplitude estimation error is:
Figure BDA0003861524950000037
where k represents the kth frequency component, variable
Figure BDA0003861524950000038
h k The position of the k-th frequency component maximum spectral line,
Figure BDA0003861524950000039
[·]expressing an integer part taking function; sinc (-) represents a sine function, abs (-) represents an absolute value solving function, max (-) represents a maximum solving function, and min (-) represents a minimum solving function; a is a 0 Representing the first coefficient of the window function.
As an improvement of the above technical solution, the objective of optimizing the minimum maximum amplitude estimation error is obtained according to a windowed discrete fourier transform of a signal to be estimated, and specifically includes the following steps:
setting a signal x to be estimated 1 (n) is composed of a plurality of frequency components, and the expression is:
Figure BDA00038615249500000310
wherein ,x1k (n) is the kth frequency component of the signal to be estimated, and the expression is:
Figure BDA0003861524950000041
wherein ,
Figure BDA0003861524950000042
representing the amplitude of the kth frequency component, kf representing the frequency of the kth frequency component, f being the fundamental frequency, i.e. the frequency of the first frequency component, at representing the sampling interval, α 1k Represents the phase of the k-th frequency component;
calculating the signal x 1 (n) windowed discrete Fourier transform
Figure BDA0003861524950000043
The calculation formula is:
Figure BDA0003861524950000044
wherein the function
Figure BDA0003861524950000045
h is the normalized frequency;
Figure BDA0003861524950000046
is the spectrum of the kth harmonic component, is the sum of two spectra, and the expression is:
Figure BDA0003861524950000047
wherein ,
Figure BDA0003861524950000048
discrete Fourier transform of window function g (n) converted by Euler formula;
only the first spectrum is of interest, namely:
Figure BDA0003861524950000049
wherein ,
Figure BDA00038615249500000410
is a phase quantity of the raw materials,
Figure BDA00038615249500000411
Δ f is the frequency resolution of the frequency,
Figure BDA00038615249500000412
wherein fs Is the sampling frequency;
when in use
Figure BDA00038615249500000413
When it is a non-integer, h = h k At maximum spectral line
Figure BDA00038615249500000414
The expression is as follows:
Figure BDA00038615249500000415
wherein ,
Figure BDA00038615249500000416
subtracting the true amplitude to obtain an amplitude estimation error as:
Figure BDA00038615249500000417
wherein ,|δk |≤0.5
Further reduction to using sinc (·) function
Figure BDA0003861524950000051
The optimization objective to minimize the maximum amplitude estimation error is:
Figure BDA0003861524950000052
as an improvement of the above technical solution, the expression of the normalized constraint condition is:
Figure BDA0003861524950000053
as an improvement of the above technical solution, the peak sidelobe level constraint condition is observed by an amplitude response;
the amplitude response is defined as:
Figure BDA0003861524950000054
the peak sidelobe level constraint is defined as:
Figure BDA0003861524950000055
where h > M +1, sidelobe is the value of a predefined peak sidelobe level constraint.
As an improvement of the above technical solution, the upper and lower bound constraints of the window function coefficient are: for each coefficient, it is required:
0≤a m ≤1,0≤m≤M
wherein M +1 is the number of coefficients, a m Represents the m +1 th coefficient.
As one of the improvements of the above technical solution, the method uses fmincon function to perform optimization solution to obtain window function coefficients, and specifically includes the following steps:
s1, setting the number of window function coefficients, initial values of the window function coefficients and initial values of sideliobe;
s2, setting parameters required by an fmincon function, wherein the Algorithm Algorithm selects an sqp Algorithm;
s3, inputting the initial value of the set window function coefficient, the sidelibe value, the parameters required by the fmincon function, the optimization objective function and each constraint condition into the fmincon function, and performing optimization solution;
s4, judging a return value exitflag of the fmincon function: if exitflag = -2, the feasible solution cannot be found, at the moment, a peak side lobe level constraint value sildelobe is automatically adjusted, sildelobe +1 is assigned to the sildelobe, the feasible domain of the solution is gradually expanded, the step S3 is carried out, and the optimization solution is carried out again; if exitflag is not equal to-2, the step S5 is carried out;
s5, further judging exitflag: if the exitflag is not equal to 1 and the exitflag is not equal to 2, the iteration times exceed the setting or the function is not converged, and at the moment, the parameters need to be reset, namely, the step S2 is skipped; and if the exitflag =1 or the exitflag =2, the feasible solution is found, and the optimized window function coefficient is output.
Compared with the prior art, the invention has the advantages that:
1. the invention is a relatively universal window function design method, can design 5 coefficients and more than window functions with rapid attenuation characteristics, and meet the requirements of different application scenes on different coefficient number window functions;
2. the window function designed by the invention has wide applicability, and is also suitable for the spectrum analysis of other signals with high dynamic range amplitude;
3. according to the window function design method, the minimum maximum amplitude estimation error is used as an optimization target, the peak side lobe level and the attenuation rate are used as constraint conditions, and the fast attenuation window function with the minimum amplitude estimation error can be obtained on the premise that the desired side lobe characteristic constraint is met preferentially. The method is more suitable for the spectrum analysis of the signals which are sampled at approximate integral multiple periods;
4. by estimating the amplitude spectrum of the inter-satellite distance signal, it can be seen that for signals with amplitudes spanning 20 orders of magnitude, very good amplitude estimation results can be obtained with maximum relative estimation errors 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 graph of amplitude spectrum estimation of a coherently sampled high dynamic range signal, wherein FIG. 3 (a) is a graph of amplitude estimation comparison of coherent sample un-windowed and non-coherent sample un-windowed, and FIG. 3 (b) is a graph of amplitude estimation comparison of coherent sample un-windowed and non-coherent sample with a new 6-term fast decay window;
FIG. 4 is a graph of amplitude estimation difference for coherent samples with no windowing for coherent samples and no windowing for incoherent samples, and FIG. 4 (b) is a graph of amplitude estimation difference for coherent samples with no windowing for incoherent samples and a new 6-term fast fading window for incoherent samples;
fig. 5 is a graph of relative percentage of difference in amplitude estimates for a high dynamic range signal with or without coherent sampling, where fig. 5 (a) is a graph of percentage of difference in amplitude estimates for coherent sample without windowing and incoherent sampling without windowing, and fig. 5 (b) is a graph of percentage of difference in amplitude estimates for coherent sample without windowing and incoherent sampling with a new 6-term fast fading window.
Detailed Description
The application provides a peak sidelobe constrained fast attenuation window function design method based on epsilon-constrained multi-objective optimization, and solves the amplitude spectrum estimation problem of signals with high dynamic range amplitude by designing a window function with fast attenuation rate, lower peak sidelobe level and lower amplitude estimation error, thereby realizing the estimation of the amplitude of signals with the amplitude range spanning more than 20 magnitude levels. In addition, verification is performed by simulated inter-satellite distance signals to evaluate the validity of the obtained window function.
The design of the window function involves a multi-objective optimization problem, which is solved in a variety of ways. The epsilon-constraint multi-objective optimization model converts an original multi-objective optimization problem into a single-objective optimization problem by converting one of the objectives into a constraint condition form, and then the problem can be solved by using a single-objective optimization method.
The window function of interest for the present invention is a combined cosine window, discrete form of which
Figure BDA0003861524950000071
Wherein, the coefficient term number is M +1, and N is the number of samples of the sample.
First, the optimization objectives and constraints of the window function need to be determined. 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 desired side lobe characteristics and minimize the optimization target of amplitude estimation error as much as possible in a feasible region.
Optimizing the target: minimizing maximum amplitude estimation error
And (3) constraint:
(1) fast attenuation constraint
(2) Peak sidelobe level constraints
(3) Normalized constraint
The following provides a formulation modeling of the optimization objectives and constraints.
The fast attenuation constraint can be derived by Discrete Fourier Transform (DFT) of the window function, and formula (1) can be converted into a fast attenuation constraint by Euler formula
Figure BDA0003861524950000072
Consider the following DFT transform
Figure BDA0003861524950000073
DFT transform on window function g (n)
Figure BDA0003861524950000081
wherein
Figure BDA0003861524950000082
Figure BDA0003861524950000083
Further simplified approximation
Figure BDA0003861524950000084
Will be 1/(h) 2 -m 2 ) Expansion of power series and substitution of the above equation
Figure BDA0003861524950000085
The attenuation rate of the side lobe is determined by the first non-zero term of the r series, so that
Figure BDA0003861524950000086
The first non-zero term of the series is obtained when r =0, then the window function is attenuated by 1/h. When in use
Figure BDA0003861524950000087
When the first non-zero term is obtained at least at r =1, the window function is then at least 1/h 3 And (4) attenuation. At the same time, when
Figure BDA0003861524950000088
When the first non-zero term is obtained at least at r =2, the window function is then at least 1/h 5 And (4) attenuation. At the same time when
Figure BDA0003861524950000089
When the first non-zero term is obtained at least at r =3, the window function is then at least 1/h 7 And (4) attenuation. At the same time, when
Figure BDA00038615249500000810
When the first non-zero term is obtained at least at r =4, the window function is at least 1/h 9 And (4) attenuation. At the same time when
Figure BDA0003861524950000091
Then, a first non-zero term is obtained at least at r =5, and the window function is at least 1/h 11 And (4) attenuation. By analogy, can continue to deduce
Figure BDA0003861524950000092
The window function is at least 1/h 13 Attenuation, when a very high attenuation rate has been reached. The constraints for better decay rates can continue to be derived, if desired.
The peak sidelobe levels of the window function can be observed by the amplitude response, which is defined as
Figure BDA0003861524950000093
The peak sidelobe level constraint can therefore be defined as
Figure BDA0003861524950000094
Where sidelibe is a predefined assumption, different sidelibs will result in different feasible fields.
The normalized constraint is defined as
Figure BDA0003861524950000095
The amplitude estimation error is obtained from the windowed DFT of the signal, assuming that the signal x to be estimated 1 (n) is composed of multiple frequency components
Figure BDA0003861524950000096
Figure BDA0003861524950000097
wherein ,
Figure BDA0003861524950000098
denotes the amplitude of the kth frequency component, k × f denotes the frequency of the kth frequency component, f is the fundamental frequency, i.e., the frequency of the first frequency component, Δ t denotes the sampling interval, N denotes the nth sampling, N denotes the total number of samples sampled, α 1k Represents the phase of the k-th frequency component;
wherein order
Figure BDA0003861524950000099
Figure BDA00038615249500000910
Is the conjugate of the two or more different compounds,
Figure BDA00038615249500000911
Δ f is the frequency resolution of the frequency,
Figure BDA0003861524950000101
wherein fs Is the sampling frequency. Can obtain
Figure BDA0003861524950000102
Calculating a signal x 1 (n) windowed DFT is represented as
Figure BDA0003861524950000103
Wherein h is the normalized frequency;
Figure BDA0003861524950000104
is a signal x 1 (n) form of windowed DFT;
Figure BDA0003861524950000105
wherein ,
Figure BDA0003861524950000106
it can be seen that the spectrum of the kth harmonic component is the sum of two spectra, we only focus on the first spectrum, i.e.
Figure BDA0003861524950000107
When in use
Figure BDA0003861524950000108
If the non-integer number is ignored, the maximum spectral line is h = h k To make
Figure BDA0003861524950000109
wherein
Figure BDA00038615249500001010
Can be obtained from the integer part of
Figure BDA00038615249500001011
From equation (5) can be obtained
Figure BDA00038615249500001012
Then for the kth frequency component, the estimated amplitude is
Figure BDA00038615249500001013
And the true amplitude is
Figure BDA00038615249500001014
The difference between the two can obtain the expression of the amplitude estimation error as
Figure BDA0003861524950000111
Further, an amplitude estimation error for the k-th frequency component can be derived as
Figure BDA0003861524950000112
The method can be further simplified into a method by utilizing a sinc (delta) function
Figure BDA0003861524950000113
Thus the optimization objective to minimize the maximum amplitude estimation error can be defined as
Figure BDA0003861524950000114
The technical scheme provided by the invention is further explained by combining the drawings and the embodiment.
The window function design process of the embodiment of the present invention is shown in fig. 1, and the process of the present invention is described below with reference to fig. 1.
1. Firstly, determining the number of window function coefficients to be optimized, and recording the number as an (M + 1) term. In order to obtain a fast fading window function, it is usually necessary to select more than 5 terms to apply more fast fading constraints, but the upper limit is not limited, but is usually better below 10 terms, which is more favorable for obtaining a smaller number of signal cycles for better spectrum analysis.
2. And (2) setting the initial value of the coefficient of the optimized time window function as a random number according to the number of the coefficients determined in the step (1), namely a0= rand (M +1, 1).
3. The initial value sidelobe level constraint is predefined, and the value of the initial value sidelobe level constraint is negative, usually-200 to-80 dB, but not limited to this interval, and needs to be determined according to the number of coefficients (M + 1) selected in step 1, and can be roughly set according to the relation of-20 (M + 1), or can be set by itself.
4. And (3) determining an optimization target, namely an expression for minimizing the maximum amplitude estimation error, namely the formula (27) according to the number of the coefficients determined in the step 1.
5. And (3) determining a normalization constraint condition, namely a formula (16), according to the number of the coefficients determined in the step 1. This condition is a constraint necessary to design the window function.
6. And determining the required fast attenuation constraint condition according to the number of the coefficients determined in the step 1. Fast fading constraintThe number of the window functions needs to be less than or equal to M, and in order to obtain the window functions with the fast attenuation characteristic, the required fast attenuation constraint conditions generally need to be more than or equal to 3, namely, the equations (8) to (10), namely, at least the window functions are ensured to be in 1/h 7 And (4) attenuation.
7. Determining the upper and lower boundary constraints of the coefficients, and requiring each coefficient to be more than or equal to 0 and less than or equal to a m ≤1,0≤m≤M。
8. And (4) carrying out optimization solution on the window function by using a nonlinear optimization function fmincon in matlab, so that parameters of the window function need to be set. Wherein the Algorithm Algorithm selects sqp, i.e. the sequential quadratic programming Algorithm, and the remaining parameters may be self-adjusting or adopt default values.
9. And after the parameter setting is finished, starting optimization solution.
10. If the optimization exit condition exitflag = -2, it means that a feasible solution cannot be found, at this time, the value of the peak side lobe level constraint sidelobe is automatically adjusted, the sidelobe level constraint sidelobe +1 is made to set sidelobe = sidelobe +1, the feasible domain of the solution is gradually expanded, and then the step 9 is returned to, and the optimization solution is carried out again. If exitflag ≠ -2, proceed to the next step, enter step 11.
11. If the optimization exit condition exitflag is not equal to 1 and the exitflag is not equal to 2, the iteration times exceed the setting or the function is not converged, and at this moment, the parameters need to be reset, namely, the step 8 is skipped.
12. And if the optimization exit condition exitflag =1 or exitflag =2 indicates that a feasible solution is found, and outputting the window function coefficient obtained by optimization.
13. And the characteristics of the window function obtained by calculation comprise peak side lobe level and maximum amplitude estimation error, and the attenuation rate is fitted to judge whether the window function obtained by optimization meets the fast attenuation constraint condition.
Taking 6 coefficients as an example, a new window function with fast attenuation characteristic is obtained by optimization according to the above process steps.
For the amplitude spectrum estimation problem of a signal with sampling time approximately being an integer multiple of the signal period (approximately 8 periods) and having a high dynamic range amplitude, 6 coefficients are selected, an initial value sidelobe level constraint of sidelobe = -120dB is defined in advance, an optimization target is determined, and a maximum amplitude estimation error is minimized to be
Figure BDA0003861524950000121
The normalized constraint is
Figure BDA0003861524950000122
To obtain a fast fading window function, 4 fading conditions are chosen, i.e.
Figure BDA0003861524950000123
Figure BDA0003861524950000124
Figure BDA0003861524950000125
Figure BDA0003861524950000131
Determining upper and lower bounds constraints for 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 of the constraint violation, tolCon, 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 computations maxfenevals is set to 200000; the maximum number of iterations maxter allowed is set to 1500. The remaining parameters assume default values.
Through optimization solution, a 6-term coefficient window function with rapid attenuation characteristic is finally obtained, and the coefficient is
Figure BDA0003861524950000132
The peak side lobe level of the new window function is-95 dB, the maximum amplitude estimation error is about 4.78%, and the new window function has the characteristic of quick attenuation and the maximum amplitude estimation error is 1/h 9 The amplitude response is shown in fig. 2.
The verification using a simulated signal with high dynamic range amplitude does not fall within the scope of the solution of the present invention, but since the validity of the new window obtained by the window function design method of the present invention needs to be verified using it, it is necessary to explain.
For signals with high dynamic range amplitudes, the space gravitational wave when considering only the gravitational field of the sun is used to detect the mission "tai chi" inter-satellite distance signal. The satellite-to-satellite distance signal at the moment is a periodic signal, but the period of the satellite-to-satellite distance signal is a non-integer, but the sampling interval can be set to be a decimal number so that the sampling time is exactly integral multiple of the signal period during simulation, namely coherent sampling is carried out, and therefore the amplitude of the satellite-to-satellite distance signal can be accurately estimated at the moment; in reality, the sampling interval is usually an integer, so that the sampling time is not completely 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-term fast attenuation window optimized by the method of the present invention is used to estimate the amplitude of the signal in the presence of leakage and compare the estimated amplitude with the result of coherent sampling.
The estimation of very good signal amplitude can be obtained 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 fast attenuation window obtained by optimization during incoherent sampling; the amplitude estimation results without windowing are also shown for comparison in the case of incoherent sampling.
The sampling duration is selected as the fundamental frequency of 8-period and inter-satellite distance signals
f 0 ＝3.168748934534878517772678686201385×10 -8 Hz
Before comparisonAmplitude estimation results of 16 multiples, i.e. contrast f k ＝k·f 0 K =1,2, the amplitude at 16. Its amplitude dynamic range is about 8X 10 6 m～4×10 -14 m, i.e. over 20 orders of magnitude.
The amplitude estimation results and differences obtained by the coherent sampling non-windowed and incoherent sampling plus new 6-term fast attenuation windows are shown in fig. 3 and fig. 4, where fig. 3 (a) is an amplitude estimation comparison graph of the coherent sampling non-windowed and incoherent sampling plus new 6-term fast attenuation windows, fig. 3 (b) is an amplitude estimation comparison graph of the coherent sampling non-windowed and incoherent sampling plus new 6-term fast attenuation windows, fig. 4 (a) is an amplitude estimation difference graph of the coherent sampling non-windowed and incoherent sampling non-windowed and fig. 4 (b) is an amplitude estimation difference graph of the coherent sampling non-windowed and incoherent sampling plus new 6-term fast attenuation windows. As can be seen from fig. 3 (a) and 4 (a), when the incoherent sample is not windowed, the amplitude estimation result is very different from the coherent sample amplitude estimation result, and as can be seen from fig. 3 (b) and 4 (b), when the incoherent sample is added with a new 6-term fast attenuation window, the amplitude estimation result is very small from the coherent sample amplitude estimation result, thereby effectively reducing the influence of leakage.
The percentage of amplitude estimation error versus the 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 non-windowed and non-coherent sample non-windowed, and fig. 5 (b) is a graph of the percentage of amplitude estimation difference for coherent sample non-windowed and non-coherent sample plus a new 6-term fast fading window; it can be seen that the new 6-term fast attenuation window obtained by using the method of the present application has a relative estimation error within 0.44% for the amplitude of a signal with a high dynamic range (the amplitude range spans about 20 orders of magnitude) in incoherent sampling, and the influence of leakage is well reduced.
As can be seen from the above detailed description of the present invention, the window function design method of the present invention uses the minimum maximum amplitude estimation error as an optimization target, and uses the peak side lobe level and the attenuation rate as constraint conditions, so that a fast attenuation window function with the minimum amplitude estimation error can be obtained on the premise of preferentially satisfying the desired side lobe characteristic constraint. The window function designed by the method is used for estimating the amplitude spectrum of the signal, and particularly for the signal with high dynamic range, the relative error of the estimated amplitude spectrum and a true value can be within 0.44% for the estimation of the amplitude spanning more than 20 magnitude levels.
Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and are not limited. Although the present invention has been described in detail with reference to the embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (10)

1. A peak side lobe constraint fast attenuation window function design method based on multi-objective optimization is disclosed, wherein a window function designed by the method is used for estimating an amplitude spectrum of a signal; the method comprises the following steps:
determining an optimization target for minimizing the maximum amplitude estimation error, a window function coefficient normalization constraint condition, a peak side lobe level constraint condition, a fast attenuation constraint condition and upper and lower limit constraint conditions of the window function coefficient according to the number of the window function coefficients to be designed;
and obtaining a window function coefficient through optimization solution according to the determined optimization target and each constraint condition, thereby completing the design of the window function.
2. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization as claimed in claim 1, wherein the window function designed by the method comprises a combined cosine window; the discrete form of the combined cosine window function g (n) is:
Figure FDA0003861524940000011
wherein M +1 is the number of window function coefficients, N is the number of samples of signal samples, a m Is the m +1 th coefficient, m =0,1,2, \ 8230And M and n represent the nth sampling.
3. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization according to claim 2, wherein the number of the fast attenuation constraint conditions is as follows: m or less and 3 or more; the fast decay constraints all satisfy the following equation:
Figure FDA0003861524940000012
under fast fading constraints, a window function of
Figure FDA0003861524940000013
Attenuation, h is the normalized frequency.
4. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization as claimed in claim 3, wherein the fast attenuation constraint condition is derived by discrete Fourier transform of a window function, and specifically comprises:
the window function g (n) is converted by the euler formula into:
Figure FDA0003861524940000014
wherein the variables are
Figure FDA0003861524940000015
Discrete Fourier transform of the transformed window function g (n)
Figure FDA0003861524940000021
The expression is as follows:
Figure FDA0003861524940000022
wherein the function
Figure FDA0003861524940000023
Function(s)
Figure FDA0003861524940000024
Further simplification results in:
Figure FDA0003861524940000025
will be 1/(h) 2 -m 2 ) The power series expansion yields:
Figure FDA0003861524940000026
and determining the decay rate of the side lobe according to the first non-zero term of the r series.
5. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization according to claim 4, wherein the expression for minimizing the maximum amplitude estimation error is as follows:
Figure FDA0003861524940000027
where k denotes the kth frequency component, variable
Figure FDA0003861524940000028
h k The position of the k-th frequency component maximum spectral line,
Figure FDA0003861524940000029
[·]representing an integer part taking function; sine (·) represents a sine function,abs (-) represents the absolute value function, max (-) represents the maximum value function, min (-) represents the minimum value function; a is a 0 Representing the first coefficient of the window function.
6. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization according to claim 5, wherein the objective of minimizing the maximum amplitude estimation error is obtained according to the windowed discrete Fourier transform of the signal to be estimated, and specifically comprises the following steps:
setting a signal x to be estimated 1 (n) is composed of a plurality of frequency components, and the expression is:
Figure FDA0003861524940000031
wherein ,x1k (n) is the kth frequency component of the signal to be estimated, and the expression is:
Figure FDA0003861524940000032
wherein ,
Figure FDA0003861524940000033
representing the amplitude of the kth frequency component, kf representing the frequency of the kth frequency component, f being the fundamental frequency, i.e. the frequency of the first frequency component, at representing the sampling interval, α 1k Represents the phase of the k-th frequency component;
calculating the signal x 1 (n) windowed discrete Fourier transform
Figure FDA0003861524940000034
The calculation formula is:
Figure FDA0003861524940000035
wherein the function
Figure FDA0003861524940000036
h is the normalized frequency;
Figure FDA0003861524940000037
is the spectrum of the kth harmonic component, is the sum of two spectra, and the expression is:
Figure FDA0003861524940000038
wherein ,
Figure FDA0003861524940000039
discrete Fourier transform of window function g (n) converted by Euler formula;
only the first spectrum is of interest, namely:
Figure FDA00038615249400000310
wherein ,
Figure FDA00038615249400000311
in order to be a phasor, the method comprises the following steps of,
Figure FDA00038615249400000312
Δ f is the frequency resolution of the frequency,
Figure FDA00038615249400000313
wherein fs Is the sampling frequency;
when in use
Figure FDA00038615249400000314
When it is a non-integer, h = h k At maximum spectral line
Figure FDA00038615249400000315
The expression is as follows:
Figure FDA00038615249400000316
wherein ,
Figure FDA00038615249400000317
subtracting the true amplitude to obtain an amplitude estimation error of:
Figure FDA0003861524940000041
wherein ,|δk |≤0.5
Further reduction to using sinc (·) function
Figure FDA0003861524940000042
The optimization objective to minimize the maximum amplitude estimation error is:
Figure FDA0003861524940000043
7. the method for designing peak sidelobe constraint fast attenuation window functions based on multi-objective optimization according to claim 6, characterized in that the expression of the normalized constraint condition is as follows:
Figure FDA0003861524940000044
8. the multi-objective optimization-based peak sidelobe constraint rapid attenuation window function design method according to claim 7, characterized in that the peak sidelobe level constraint condition is observed through amplitude response;
the amplitude response is defined as:
Figure FDA0003861524940000045
the peak sidelobe level constraint is defined as:
Figure FDA0003861524940000046
where h > M +1, sideliobe is a predefined value of the peak sidelobe level constraint.
9. The method for designing the peak sidelobe constraint fast attenuation window function based on the multi-objective optimization according to claim 8, wherein the upper and lower bound constraint conditions of the window function coefficient are as follows: for each coefficient, it is required:
0≤a m ≤1,0≤m≤M
wherein M +1 is the number of coefficients, a m Represents the m +1 th coefficient.
10. The multi-objective optimization-based peak sidelobe constraint fast attenuation window function design method as claimed in claim 9, wherein the method uses fmincon function to perform optimization solution to obtain window function coefficients, and specifically comprises the following steps:
s1, setting the number of window function coefficients, initial values of the window function coefficients and initial values of sideliobe;
s2, setting parameters required by the fmincon function, wherein the Algorithm selects the sqp Algorithm;
s3, inputting the initial value of the set window function coefficient, the sideliobe value, the parameters required by the fmincon function, the optimization objective function and each constraint condition into the fmincon function for optimization solution;
s4, judging a return value exitflag of the fmincon function: if exitflag = -2, the feasible solution cannot be found, at the moment, a peak side lobe level constraint value sildelobe is automatically adjusted, sildelobe +1 is assigned to the sildelobe, the feasible domain of the solution is gradually expanded, the step S3 is carried out, and the optimization solution is carried out again; if exitflag is not equal to-2, the step S5 is carried out;
s5, further judging exitflag: if the exitflag is not equal to 1 and the exitflag is not equal to 2, the iteration times exceed the setting or the function is not converged, and at the moment, the parameters need to be reset, namely, the step S2 is skipped; and if the exitflag =1 or the exitflag =2, the feasible solution is found, and the optimized window function coefficient is output.
CN202211164323.5A 2022-09-23 2022-09-23 Peak sidelobe constraint rapid decay window function design method based on multi-objective optimization Active CN115510787B (en)

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 true CN115510787A (en) 2022-12-23
CN115510787B 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 (7)

* Cited by examiner, † Cited by third party
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
CN109598093A (en) * 2018-12-29 2019-04-09 北京化工大学 Earthquake vector wave field numerical method and system based on fitting window function
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
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
CN114116370A (en) * 2021-08-31 2022-03-01 西南电子技术研究所（中国电子科技集团公司第十研究所） Method for optimizing operation health state monitoring points of complex electronic system
CN114966201A (en) * 2022-05-05 2022-08-30 威胜电气有限公司 Harmonic voltage and current measurement and characteristic analysis method in typical scene

Patent Citations (7)

* Cited by examiner, † Cited by third party
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
CN109598093A (en) * 2018-12-29 2019-04-09 北京化工大学 Earthquake vector wave field numerical method and system based on fitting window function
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
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
CN114116370A (en) * 2021-08-31 2022-03-01 西南电子技术研究所（中国电子科技集团公司第十研究所） Method for optimizing operation health state monitoring points of complex electronic system
CN114966201A (en) * 2022-05-05 2022-08-30 威胜电气有限公司 Harmonic voltage and current measurement and characteristic analysis method in typical scene

Also Published As

Publication number Publication date
CN115510787B (en) 2023-04-21

Similar Documents

Publication Publication Date Title
US5357257A (en) Apparatus and method for equalizing channels in a multi-channel communication system
US7161515B2 (en) Calibration system and method for a linearity corrector using filter products
US7609759B2 (en) Method and system of nonlinear signal processing
Vanbeylen et al. Blind maximum-likelihood identification of Wiener systems
EP1128293A2 (en) Modelling non-linear devices
Tóth et al. Model structure learning: A support vector machine approach for LPV linear-regression models
CN111950759A (en) Short-term wind speed prediction method based on two-stage decomposition, LSTM and AT
CN111553513A (en) Medium-and-long-term runoff prediction method based on quadratic decomposition and echo state network
US7457756B1 (en) Method of generating time-frequency signal representation preserving phase information
US6892155B2 (en) Method for the rapid estimation of figures of merit for multiple devices based on nonlinear modeling
Petri Frequency-domain testing of waveform digitizers
CN113189561A (en) Sea clutter parameter estimation method, system, equipment and storage medium
CN111159888B (en) Covariance matrix sparse iteration time delay estimation method based on cross-correlation function
CN115062668A (en) Harmonic parameter detection method and system based on RAdam optimization width learning
CN114842867A (en) DFT-based audio sinusoidal signal frequency estimation method and system
CN115510787A (en) Peak sidelobe constraint rapid attenuation window function design method based on multi-objective optimization
CN111917676B (en) Linear frequency modulation interference cancellation method
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
Liepin’sh An algorithm for evaluation a discrete Fourier transform for incomplete data
CN106549652A (en) Filter coefficient update in time-domain filtering
Neumayer et al. Bandwidth Based Design Methodology for Wiener Filters for Online Signal Denoising
CN113804981B (en) Time-frequency joint optimization multi-source multi-channel signal separation method
CN101741776A (en) Method and device for eliminating interference signals
CN114330420B (en) Data-driven radar communication aliasing signal separation method and device

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