CN107622035B - Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing - Google Patents
Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing Download PDFInfo
- Publication number
- CN107622035B CN107622035B CN201710914700.5A CN201710914700A CN107622035B CN 107622035 B CN107622035 B CN 107622035B CN 201710914700 A CN201710914700 A CN 201710914700A CN 107622035 B CN107622035 B CN 107622035B
- Authority
- CN
- China
- Prior art keywords
- frequency
- signal
- polynomial phase
- phase signal
- time
- 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
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
- Complex Calculations (AREA)
Abstract
The invention provides a polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing, which can complete time-frequency decomposition of a polynomial phase signal, wherein each signal component obtained by decomposition is a single component which only corresponds to one frequency point at any moment, then the signal frequency distribution corresponding to the corresponding moment is directly calculated and generated by using each signal component and the instantaneous frequency value of each moment through a Sinc function only retaining main lobe response, the defect that non-single components of a plurality of frequency points corresponding to one moment in the traditional time-frequency transformation have cross terms is overcome, and finally time-frequency distribution which does not have any cross term interference and has better time-frequency joint resolution is output; the method has simple principle and convenient operation, can effectively overcome the adverse effect of cross term interference of the classical time-frequency analysis method and the loss of time-frequency joint resolution, and can effectively improve the quality and benefit of time-frequency analysis of the non-stationary polynomial phase signals.
Description
Technical Field
The invention belongs to the field of signal analysis, and particularly relates to a polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing.
Background
Many natural and artificial signals, such as speech, biomedical signals, waves propagating in dispersive media, mechanical vibrations, animal sounds, music, radar, sonar signals, etc., are typically non-stationary signals characterized by finite duration and time-varying frequency, and are characterized by non-stationary, non-linear, non-uniform, non-structural, non-deterministic, non-integrable, non-reciprocal, amorphous, non-regular, non-continuous, non-smooth, non-periodic, non-symmetric characteristics. The joint time-frequency analysis (abbreviated as time-frequency analysis) is to express a one-dimensional time signal in a form of a two-dimensional time-frequency density function by focusing on the time-varying characteristics of the components of the real signal, and aims to disclose how many frequency components are included in the signal and how each component varies with time.
In 1948, french scholars j.ville introduced the Wigner distribution, proposed in 1932 by american physicist e.p. Wigner who was born to budapest hungary, into the field of signal processing, and came to be called "Wigner-Ville distribution" (WVD). The follow-up scholars are imitated and provide some novel time-frequency distribution. The whole time-frequency analysis history is almost a history of struggling with the deficiency of WVD. According to the essential characteristics of each group, the time-frequency distribution of the shapes and colors can be classified into the following categories: (1) linear time-frequency representation; (2) cohen-like bilinear time-frequency distribution; (3) affine bilinear time frequency distribution; (4) rearranging bilinear time frequency distribution; (5) self-adaptive kernel function class time frequency distribution; (6) and (4) parameterizing time-frequency distribution.
The time-frequency resolution of the linear time-frequency transformation Gabor transformation and STFT is limited by the shape and width of the window function. Wavelet analysis is essentially a time-scale analysis, more suitable for analyzing signals with self-similar structures (e.g. fractal) and abrupt (transient) signals, while the results of wavelet transform are often difficult to interpret from the point of view of characterizing the time-varying structure of a signal.
The essence of Cohen-like bilinear time-frequency distribution is to distribute the energy of a signal (some square form of the signal) in a time-frequency plane, which is based on WVD. However, WVD is not linear, i.e., the WVD for the sum of the two signals is not the sum of the WVD for each signal, with an additional term added. The effect of the cross terms on WVD is strong, and spots are visible. The cross terms are real, are mixed among self-term components and have larger amplitude; in addition, the cross terms are oscillatory in type, with one cross term being generated every two signal components. If the signal has N components, it will generate
And (4) a cross item. The existence of cross terms seriously interferes with the interpretation of the WVD; the time-frequency distribution given by WVD becomes even meaningless when the signal components become complex.
In order to solve the influence of cross-term interference, several different time-frequency analysis methods have been proposed in succession, wherein the Cohen-type kernel time-frequency analysis method generates time-frequency distribution with desired characteristics by designing a two-dimensional kernel (two-dimensional filter). However, the class of time-frequency distribution uses a smoothing method to suppress cross terms at the expense of the time-frequency resolution of the whole distribution.
The description of the local characteristics of the signal in time and frequency by the various distributions (except wavelet transform) described above is realized by time shift and frequency shift transform; in contrast, the distribution of affine classes is achieved by time-shifting and scaling.
The most well-known distribution in affine class is the quadratic of the wavelet transform. Since the basis for this type of distribution is still WVD, WVD is naturally one of its members. In fact, WVD is just a tie connecting Cohen classes and affine classes. The former is based on time-frequency smoothing (time-frequency smoothing) performed on the WVD, and the latter is based on affine smoothing (affine smoothing).
The Cohen-like and affine-like time-frequency distributions can greatly suppress cross term interference by performing time-frequency smoothing and time-scale smoothing on the WVD, as shown in FIG. 4, but still a lot of cross terms remain, and some distributions introduce some new cross terms. To further improve the performance of the two types of distributions, k.kdoea and the like firstly propose the idea of rearranging the time-frequency plane, and then f.auger and p.flandrin develop and perfect the rearrangement method.
In bilinear time-frequency distributions of the Cohen class, the affine class and the rearrangement class, each distribution corresponds to a fixed kernel function, and the kernel function determines the cross term suppression characteristic of the corresponding distribution. It goes without saying that a kernel function is only valid for one class of signals, and hence the bilinear time-frequency distribution in all the three classes lacks adaptability to signals.
The various time-frequency analysis methods described above are non-parametric methods, and they do not a priori assume what model signal the signal is composed of. And a parametric time frequency analysis (parametric time frequency analysis) method constructs a signal model which is best matched with the signal hierarchical structure according to the analysis of the signal hierarchical structure, thereby concentrating the information of the signal, simplifying the representation of the signal and obtaining the time frequency distribution of the signal.
In the linear time-frequency representation (atom decomposition) method, if the selected atoms are similar to the main components of the signal, only linear combination of a few atoms is needed to represent the signal more accurately, and the decomposition result is sparse (sparse). On the contrary, if the character of the atom is far from or different from the main structure of the signal, a large number of atoms, even infinite number of atoms, are required to be assembled into the original signal accurately enough, and the information of the signal is dispersed on too many atoms, which is not favorable for effectively representing the signal. Therefore, when the atom decomposition method is adopted, the combination of atoms must be selected adaptively according to the local structural characteristics of the signal, so as to decompose the signal with as few atoms as possible. In 1994, s.qian and d.chen pioneer the parametric time-frequency analysis (the inventive idea goes back to 1988), and an adaptive expansion algorithm is proposed; mallat and z zhang proposed a sister algorithm with adaptive unfolding- "matchingpursuit" algorithm in 1993. The two algorithms, although named differently, have essentially no irreconcilable outliers and are essentially equivalent.
The essence of the adaptive matching projection tower decomposition algorithm is that the time-frequency energy distribution of atoms is used for approximating the time-frequency energy distribution of an original signal. Since the s.qian and s.mallat use Gabor atoms whose frequencies are not changed, the partitioning of the time-frequency plane by the iterative algorithm is a lattice segmentation. This algorithm works well for frequency components that are invariant, but when the signal to be analyzed is a Chirp signal, the matching is equivalent to zeroth order curve approximation, which inevitably causes many truncations and mixed distortions between the components in the decomposition process. To overcome this drawback, s.mann and s.haykin et al propose almost simultaneously to use the gaussian function of stretching, time shifting, frequency shifting and frequency tilting, a "chirp" as an atom to replace the Gabor atom with constant frequency (s.mann adds another time tilting operation), and to use the inner product method to obtain the so-called "chirp transform".
The Chirplet transformation based on the adaptive matching projection turriform decomposition algorithm is essentially to perform linear approximation on any energy curve on a time-frequency plane by using a group of line segments with any inclination. Obviously, the first order approximation can express chirp-like signals more compactly than the zero order approximation. While scholars have been fond of linearity, nature is often non-linear. When the frequency components of a signal vary non-linearly with time, such as a naturally occurring or artificially generated Doppler signal, characterized by atoms whose frequencies vary linearly (zeroth or first order) with time, the number of atoms must be increased, thereby affecting both the understanding and interpretation of the decomposition results and the data compression capabilities of the decomposition results.
However, the above-mentioned adaptive time-frequency decomposition method is based on a known signal hierarchy or signal model, and is adaptively decomposed in signal parameters, but is non-blind in signal type. New attempts are needed to further improve the above-mentioned parametric time-frequency analysis methods.
Time-frequency analysis has been used with many success in speech recognition, radar signal processing and image processing, seismic signal processing, signal reconstruction, and interference suppression in spread spectrum communications. In general, the application fields of time-frequency analysis are roughly classified into four types: firstly, time-varying spectrum analysis; secondly, some physical quantities are indirectly calculated through time-frequency distribution; thirdly, the time-frequency distribution is used as a carrier of information carried by the signal (without concern whether it can really represent the energy density); fourthly, signal reconstruction, compression, encoding and the like.
The result of the wavelet transform is very confusing from the viewpoint of examining the time-dependent evolution of the frequency components of the signal, although this field is good. The wavelet transform is a multi-resolution analysis method, which takes time and scale as parameters and has different resolutions at different positions of a time-scale plane. The wavelet analysis benefits from the completeness, self-similarity and multi-resolution of wavelet atoms, and can obtain two most important reasons for success, namely that the wavelet analysis has a tower-shaped rapid algorithm and good time-frequency local characteristics; the disadvantage is that once the mother wavelet is not properly selected, the application effect is greatly influenced. From the perspective of signal compression and cross term interference elimination, the parameterized time-frequency analysis method is better, but the calculation of the parameters of the atomic model is not easy. WVD and all other Cohen-like time-frequency distributions can be used for analyzing narrowband signals, but they are not well suited for analyzing multi-component wideband signals as well as radar and sonar signals, and if affine-like distributions are used, the problem of resolution space-variant exists. In addition, various time-frequency analysis techniques are difficult to be distinguished, and the key is the type of signals suitable for the time-frequency analysis techniques. A frequent trial and error experience is that a very fast running spectrum of STFT can be tried first, as shown in fig. 3; if a higher time-frequency resolution is required, a parameterized time-frequency analysis method can be adopted.
However, the existing time-frequency analysis methods are not suitable for phase modulation of polynomial phase signals which can be expressed as finite-term polynomial series.
The patent "a model-driven polynomial phase signal adaptive time-frequency transformation method" provides a feasible new idea and method for the above-mentioned practical problems, but does not specifically provide specific practical operation and operation process of various modern optimization algorithms. The invention provides a new self-adaptive time-frequency decomposition algorithm for a non-stationary polynomial phase signal, which is a specific implementation mode, a practical operation process and a calculation process of a new model-driven polynomial phase signal self-adaptive time-frequency decomposition method, can self-adaptively finish the time-frequency decomposition of the polynomial phase signal, and output time-frequency distribution which has no cross term interference and better time-frequency joint resolution.
Disclosure of Invention
In order to solve the problems, the invention provides a polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing, which utilizes instantaneous frequency values of each signal component and each moment, directly calculates and generates signal frequency distribution corresponding to the corresponding moment by only reserving a Sinc function of main lobe response, overcomes the defect that non-single components of a plurality of frequency points corresponding to one moment in the traditional time-frequency transformation have cross terms, and finally outputs the time-frequency distribution which has no cross term interference and better time-frequency joint resolution.
A polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing comprises the following steps:
step 1: performing multiple decompositions on an original polynomial phase signal s (t) received by a radar by using a simulated annealing algorithm, wherein each decomposition determines an optimal model order and an optimal coefficient set to be determined of a signal component, and specifically:
step 11: generating a reference function h corresponding to the original polynomial phase signal s (t) according to the polynomial phase modelpAnd initializing a reference function hpModel order of (N)11, and calculating the energy E of the original polynomial phase signal s (t)0Randomly generating a reference function hpTwo sets of coefficients to be determined { a }n},n＝0,1,2,...,N1Wherein a set of coefficients { a } to be determinednTaking the solution as an initial solution in a simulated annealing algorithm;
step 12: two coefficient sets to be determined { a }nH corresponding topThe conjugate of the signal is multiplied with an original polynomial phase signal s (t) respectively to obtain a mixed modulation signal x (t), Fourier transform is carried out on the mixed modulation signal x (t) to obtain a transform result X (f), a characteristic value of one frequency domain characteristic in the transform result X (f) is calculated, the characteristic value is used as the fitness of the evaluation transform result X (f), and the fitness is good or bad according to each transform result X (f), wherein the transform result X (f) with the better fitness corresponds to the current optimal solution Gbest, the transform result X (f) with the worse fitness is the secondary optimal solution PreBest, and the current optimal solution Gbest is assigned to the pre;
step 13: according to the pre, according to the set annealing step factor TstepCalculating a new possible solution Nextp, and calculating the fitness fit of the possible solution Nextp by the method of step 12nextThe fitness fit of Nextp is solved if possiblenextReplacing the fitness of the current optimal solution Gbest which is superior to the fitness of the current optimal solution Gbest in the step 12 with the current optimal solution Gbest, and replacing the original optimal solution with a secondary optimal solution PreBest, otherwise, keeping the current optimal solution Gbest and the secondary optimal solution PreBest unchanged;
step 14: randomly generating another new possible solution Nextq and corresponding undetermined coefficient set { a }nCalculating the fitness fit of the new possible solution Nextq corresponding to the transformation result X (f) according to the method of step 12newIf the fitness fit is fitnewBetter than fitness fitnextReplacing Nextq of the new possible solution with pre, otherwise according to the replacement probability pωThe magnitude of the random probability rand is used to decide whether to replace Nextq with Prelocal, if the probability p is replacedωIf the probability is greater than the random probability rand, replacing Nextq with the pre, otherwise, keeping the pre unchanged;
step 15, repeating the step 13 to the step 14 until the cycle number reaches L markov chain values, so as to obtain the current optimal solution Gbest and the suboptimal solution PreBest after iteration;
step 16: calculating the fitness difference value delta fit of the current optimal solution Gbest and the suboptimal solution PreBest obtained in the step 15, and calculating the annealing temperature TtWherein the annealing temperature TtDecays as the number of iterations of step 16 increases, if the absolute value of the fitness difference Δ fit is greater than the set annealing tolerance ξ, and the annealing temperature TtGreater than a set terminal temperature TlowThen steps 13-15 are repeated until the absolute value of the fitness difference Δ fit is not greater than the set annealing tolerance ξ, or the annealing temperature TtNot greater than a set cut-off temperature TlowSo as to obtain the final optimal solution Gtest';
and step 17: order of model N1Sequentially obtaining the maximum possible order _ max from 2, and repeating the steps 11 to 16 to obtain the final optimal solution Gbest' under different model orders;
step 18: selecting a solution with optimal fitness from all optimal solutions Gbest' obtained in the step 17 so as to determine the optimal model order N of the current signal component of the original polynomial phase signal s (t)pAnd NpCorresponding optimal pending coefficient set { an}max；
Step 2: using the optimal model order NpAnd an optimal pending coefficient set { a }n}maxConjugate of determined signal component and original polynomialMultiplying the phase signals s (t) and performing Fourier transform to obtain a frequency spectrum X '(f), and performing inverse Fourier transform after zeroing the intensity complex value at the maximum value of the envelope of the frequency spectrum X' (f) to obtain a time domain signal y (t);
and step 3: using model order of NpAnd an optimal pending coefficient set { a }n}maxMultiplying the determined signal component by a time domain signal y (t) to obtain a residual signal z (t), and finishing the decomposition;
and 4, step 4: calculating the energy E of the residual signal z (t)dWith the energy E of the original polynomial phase signal s (t)0Taking the ratio R, if the ratio R is less than the set threshold gamma or the decomposition times reaches the set upper limit number NmaxStopping decomposition to obtain the optimal model order N of each component signalpAnd an optimal pending coefficient set { a }n}max(ii) a Otherwise, replacing the original polynomial phase signal s (t) in the step 1 with the residual signal z (t) calculated in the step 3 to recalculate the mixed modulation signal x (t), and repeating the steps 1 to 3 until the ratio R is smaller than the set threshold gamma or the decomposition times reaches the set upper limit number Nmax；
And 5: the optimal model order NpAnd an optimal pending coefficient set { a }n}maxThe corresponding signal components are numbered h in the order of decompositioncWherein C is 1, 2.., C, wherein C is the number of adaptive decompositions;
step 6: obtaining each signal component h according to the physical definition of the instantaneous frequency of the single component signalcCorresponding instantaneous frequency curve fc(t) and from all signal components hcInstantaneous frequency curve fc(t) determining the frequency distribution range of the original polynomial phase signal s (t), and finally performing discretization, main lobe response reservation and superposition operations on the frequency distribution range in sequence to obtain the time-frequency joint distribution f (t) of the original polynomial phase signal s (t).
Further, the original polynomial phase signal s (t) corresponds to a reference function hpThe method specifically comprises the following steps:
where order _ max is the maximum possible order of the original polynomial phase signal s (t).
Further, the frequency distribution range of the original polynomial phase signal s (t) is specifically:
all signal components hcInstantaneous frequency curve fcThe instantaneous frequency minimum in (t) is the minimum frequency f of the original polynomial phase signal s (t)iminThe instantaneous frequency maximum is the maximum frequency f of the original polynomial phase signal s (t)imax。
Further, the set of coefficients to be determined { a }nIs specifically a vector
Wherein N is1Representing the model order of the current optimization, a0Setting the position serial number of a sampling point where a peak value obtained by direct Fourier transform of an original polynomial phase signal s (t) is located for a fixed value;
the coefficient to be optimized is generated by a random function Rand with the value range of [ -M, M]Where M is the number of sampling points of the original polynomial phase signal s (t), specifically:
wherein rand (2, N)1) Representing random generation of 2 dimensions N1Distribution of (2) in [0,1 ]]To obtain a model order of N1The initial solution of (a).
Further, the frequency domain feature in step 12 is any one of a spectrum peak, a contrast of a spectrum, or a spectrum shape entropy of a spectrum, where a feature value of each frequency domain feature is a spectrum peak XpSpectral contrast XcAnd the spectral shape entropy X of the frequency spectrumeRespectively as follows:
Xp＝max{abs[X(f)]}
Xe＝sum(-pS·logpS)
wherein
sum (-) represents the summation operation of the vectors, abs (-) is the absolute value;
when selecting spectral peak XpOr spectral contrast XcWhen the fitness is larger, the evaluation result is better, and the spectrum shape entropy X of the selected frequency spectrumeIn the case of fitness, the smaller the value, the better the evaluation result.
Further, the step 13 is according to the set annealing step factor TstepCalculating a new possible solution Nextp, wherein the specific method is as follows:
Nextp＝Prelocal+Tstep·[2M·rand(1,N1)-M]
wherein rand (1, N)1) Representing random generation of 1 dimension N1Distribution of (2) in [0,1 ]]M is the number of samples of the original polynomial phase signal s (t).
Further, the replacement probability p of step 14ωThe specific calculation method comprises the following steps:
pω＝exp[-(fitnext-fitnew)]。
further, the annealing temperature T in step 16tThe specific calculation method comprises the following steps:
Tt＝KtT0
wherein K is a temperature decay parameter of the set simulated annealing, T0T is the number of repetitions of step 16 for a set simulated annealing initiation temperature, and is at least 1.
Further, the step 6 of sequentially discretizing the frequency distribution range specifically includes:
fis＝[fiminfimin+Δf fimin+2Δf fimin+3Δf ... fimax]T
wherein f isisFor the instantaneous frequency value of the discrete frequency vector, Δ f is the set frequency resolution, specifically:
wherein f isiminIs the minimum frequency, f, of the original polynomial phase signal s (t)imaxThe maximum frequency of the original polynomial phase signal s (t), K is the set discrete frequency dimension, and i is the imaginary unit.
Further, the specific calculation method for reserving the main lobe response for the frequency distribution range in step 6 is as follows:
wherein A iscFor the c-th signal component hcI is the imaginary unit, deltaf is the set frequency resolution, f is the intensity complex value at the maximum of the spectral envelope of (1)isAs a discrete frequency vector fsInstantaneous frequency value of fic(t) is the respective signal component hcCorresponding instantaneous frequency curve fc(t) instantaneous frequency value.
Has the advantages that:
the invention provides a polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing, which can complete time-frequency decomposition of a polynomial phase signal, wherein each signal component obtained by decomposition is a single component which only corresponds to one frequency point at any moment, then the signal frequency distribution corresponding to the corresponding moment is directly calculated and generated by using each signal component and the instantaneous frequency value of each moment through a Sinc function only retaining main lobe response, the defect that non-single components of a plurality of frequency points corresponding to one moment in the traditional time-frequency transformation have cross terms is overcome, and finally time-frequency distribution which does not have any cross term interference and has better time-frequency joint resolution is output;
the method has simple principle and convenient operation, can effectively overcome the adverse effect of cross term interference of the classical time-frequency analysis method and the loss of time-frequency joint resolution, and can effectively improve the quality and benefit of time-frequency analysis of the non-stationary polynomial phase signals.
Drawings
FIG. 1 is a flow chart of the adaptive time-frequency transform algorithm based on simulated annealing according to the present invention;
FIG. 2 is a time domain waveform of a multi-component polynomial phase signal according to the present invention;
FIG. 3 is a time-frequency diagram of a classical STFT analysis in the prior art;
FIG. 4 is a typical WVD time-frequency diagram of the prior art;
FIG. 5 is a time-frequency diagram obtained by the adaptive time-frequency transform algorithm based on simulated annealing according to the present invention.
Detailed Description
The invention is described in detail below by way of example with reference to the accompanying drawings.
The method comprises the steps of firstly estimating the model order of each component of the polynomial phase signal and corresponding phase parameters of each order through simulated annealing, then extracting corresponding polynomial phase signal components by using the cleaning thought, removing the signal components from the original signal to obtain a residual signal, then repeatedly using the simulated annealing to implement self-adaptive model order determination and optimal model parameter optimization on the residual signal in an iterative cleaning mode, gradually extracting each signal component, and repeating the steps until the energy of the residual signal is lower than a preset threshold or the number of the extracted signal components reaches a preset maximum value. Then, constructing a phase-time history by using phase parameters of each component signal obtained by simulated annealing decomposition, deriving the phase history according to the physical definition of the instantaneous frequency of the single component signal, directly obtaining the frequency-time history of each signal component in an analytic mode, namely, an instant frequency change curve, determining the maximum and minimum frequencies of the whole signal according to the time-frequency change curve of each component of the signal, determining the frequency change range of the whole signal by using the maximum and minimum frequencies, discretizing the frequency range according to the frequency resolution requirement required by application to obtain a discrete frequency vector of a required dimension, and finally, directly calculating and generating the signal frequency distribution corresponding to the corresponding moment by using the instantaneous frequency values of each component and each moment and only keeping a Sinc function of main lobe response; repeating the steps until the frequency distribution of all the signal components at all the moments is completely generated, storing the signal components in time sequence and accumulating the signal components one by one to obtain the time-frequency joint distribution of the final polynomial phase signals; the method is realized by utilizing a simulated annealing process, has simple principle and convenient operation, can effectively overcome the adverse effect of cross term interference and the loss of time-frequency joint resolution of a classical time-frequency analysis method, and improves the quality and benefit of time-frequency analysis of non-stable polynomial phase signals.
As shown in fig. 1, a polynomial phase signal adaptive time-frequency transform method based on simulated annealing includes the following steps:
step 1: performing multiple decompositions on an original polynomial phase signal s (t) received by a radar by using a simulated annealing algorithm, wherein each decomposition determines an optimal model order and an optimal coefficient set to be determined of a signal component, and specifically:
step 11: generating a reference function h corresponding to the original polynomial phase signal s (t) according to the polynomial phase modelpAnd initializing a reference function hpModel order of (N)11, and calculating the energy E of the original polynomial phase signal s (t)0Randomly generating a reference function hpTwo sets of coefficients to be determined { a }n},n＝0,1,2,...,N1Wherein a set of coefficients { a } to be determinednTaking the solution as an initial solution in a simulated annealing algorithm;
step 12: two coefficient sets to be determined { a }nH corresponding topRespectively multiplying the conjugate of the signal with an original polynomial phase signal s (t) to obtain a mixed modulation signal x (t), performing Fourier transform on the mixed modulation signal x (t) to obtain a transform result X (f), calculating a characteristic value of one frequency domain characteristic of the transform result X (f), taking the characteristic value as the fitness of the evaluation transform result X (f), and obtaining the fitness of each transform result X (f), wherein the initial solution corresponding to the transform result X (f) with better fitness is the current optimal solution Gbest and the fitness is betterThe corresponding initial solution of the poor transformation result X (f) is the secondary optimal solution PreBest, and simultaneously the current optimal solution Gbest is assigned to the pre;
step 13: according to the pre, according to the set annealing step factor TstepCalculating a new possible solution Nextp, and calculating the fitness fit of the possible solution Nextp by the method of step 12nextThe fitness index fit of Nextp if possiblenextReplacing the fitness of the current optimal solution Gbest which is superior to the fitness of the current optimal solution Gbest in the step 12 with the current optimal solution Gbest, and replacing the original optimal solution with a secondary optimal solution PreBest, otherwise, keeping the current optimal solution Gbest and the secondary optimal solution PreBest unchanged;
step 14: carrying out the Metropolis process: another new possible solution Nextq is randomly generated and its corresponding set of pending coefficients { a }nCalculating the fitness fit of the new possible solution Nextq corresponding to the transformation result X (f) according to the method of step 12newIf the fitness fit is fitnewBetter than fitness fitnextReplacing Nextq of the new possible solution with pre, otherwise according to the replacement probability pωThe magnitude of the random probability rand is used to decide whether to replace Nextq with Prelocal, if the probability p is replacedωIf the probability is greater than the random probability rand, replacing Nextq with the pre, otherwise, keeping the pre unchanged;
step 15, repeating the step 13 to the step 14 until the cycle number reaches L markov chain values, so as to obtain the current optimal solution Gbest and the suboptimal solution PreBest after iteration;
step 16: calculating the fitness difference value delta fit of the current optimal solution Gbest and the suboptimal solution PreBest obtained in the step 15, and calculating the annealing temperature TtWherein the annealing temperature TtDecays as the number of iterations of step 16 increases, if the absolute value of the fitness difference Δ fit is greater than the set annealing tolerance ξ, and the annealing temperature TtGreater than a set terminal temperature TlowThen steps 13-15 are repeated until the absolute value of the fitness difference Δ fit is not greater than the set annealing tolerance ξ, or the annealing temperature TtNot greater than a set cut-off temperature TlowThereby obtainingThe final optimal solution Gbest';
and step 17: order of model N1Sequentially obtaining the maximum possible order _ max from 2, and repeating the steps 11 to 16 to obtain the final optimal solution Gbest' under different model orders;
step 18: selecting a solution with optimal fitness from all optimal solutions Gbest' obtained in the step 17 so as to determine the optimal model order N of the current signal component of the original polynomial phase signal s (t)pAnd NpCorresponding optimal pending coefficient set { an}max；
Step 2: using the optimal model order NpAnd an optimal pending coefficient set { a }n}maxMultiplying the determined conjugate of the signal component by the original polynomial phase signal s (t) and performing fourier transform to obtain a frequency spectrum X '(f), and performing inverse fourier transform after zeroing the intensity complex value at the maximum value of the envelope of the frequency spectrum X' (f) to obtain a time domain signal y (t), as shown in fig. 2;
and step 3: using model order of NpAnd an optimal pending coefficient set { a }n}maxMultiplying the determined signal component by a time domain signal y (t) to obtain a residual signal z (t), and finishing the decomposition;
and 4, step 4: calculating the energy E of the residual signal z (t)dWith the energy E of the original polynomial phase signal s (t)0Taking the ratio R, if the ratio R is less than the set threshold gamma or the decomposition times reaches the set upper limit number NmaxStopping decomposition to obtain the optimal model order N of each component signalpAnd an optimal pending coefficient set { a }n}max(ii) a Otherwise, replacing the original polynomial phase signal s (t) in the step 1 with the residual signal z (t) calculated in the step 3 to recalculate the mixed modulation signal x (t), and repeating the steps 1 to 3 until the ratio R is smaller than the set threshold gamma or the decomposition times reaches the set upper limit number Nmax；
And 5: the optimal model order NpAnd an optimal pending coefficient set { a }n}maxThe corresponding signal components are numbered h in the order of decompositioncWherein C is 1, 2.., C, wherein C is the number of adaptive decompositions;
step 6: obtaining each signal component h according to the physical definition of the instantaneous frequency of the single component signalcCorresponding instantaneous frequency curve fc(t) and from all signal components hcInstantaneous frequency curve fc(t) determining the frequency distribution range of the original polynomial phase signal s (t), and finally performing discretization, main lobe response reservation and superposition operations on the frequency distribution range in sequence to obtain the time-frequency joint distribution f (t) of the original polynomial phase signal s (t), as shown in fig. 5.
Further, the original polynomial phase signal s (t) corresponds to a reference function hpThe method specifically comprises the following steps:
wherein, the set of coefficients to be determined { a }n},n＝0,1,2,...N1Undetermined coefficient of (a)nThe initial values are randomly generated and order _ max is the maximum possible order of the original polynomial phase signal s (t).
Further, the frequency distribution range of the original polynomial phase signal s (t) is specifically:
all signal components hcInstantaneous frequency curve fcThe instantaneous frequency minimum in (t) is the minimum frequency f of the original polynomial phase signal s (t)iminThe instantaneous frequency maximum is the maximum frequency f of the original polynomial phase signal s (t)imax。
Further, the set of coefficients to be determined { a }nIs specifically a vector
Wherein N is1Representing the model order of the current optimization, a0Setting the position serial number of a sampling point where a peak value obtained by direct Fourier transform of an original polynomial phase signal s (t) is located for a fixed value;
the coefficient to be optimized is generated by a random function Rand and is subjected to value rangeIs surrounded by [ -M, M [)]Where M is the number of sampling points of the original polynomial phase signal s (t), specifically:
wherein rand (2, N)1) Representing random generation of 2 dimensions N1Distribution of (2) in [0,1 ]]To obtain a model order of N1The initial solution of (a).
Further, the frequency domain feature in step 12 is any one of a spectrum peak, a contrast of a spectrum, or a spectrum shape entropy of a spectrum, where a feature value of each frequency domain feature is a spectrum peak XpSpectral contrast XcAnd the spectral shape entropy X of the frequency spectrumeRespectively as follows:
Xp＝max{abs[X(f)]}
Xe＝sum(-pS·logpS)
wherein
sum (-) represents the summation operation of the vectors, abs (-) is the absolute value;
when selecting spectral peak XpOr spectral contrast XcWhen the fitness is larger, the evaluation result is better, and the spectrum shape entropy X of the selected frequency spectrumeIn the case of fitness, the smaller the value, the better the evaluation result.
Further, the step 13 is according to the set annealing step factor TstepCalculating a new possible solution Nextp, wherein the specific method is as follows:
Nextp＝Prelocal+Tstep·[2M·rand(1,N1)-M]
wherein rand (1, N)1) Representing random generation of 1 dimension N1Distribution of (2) in [0,1 ]]M is the originalThe number of samples of the polynomial phase signal s (t).
Further, the replacement probability p of step 14ωThe specific calculation method comprises the following steps:
pω＝exp[-(fitnext-fitnew)]。
further, the annealing temperature T in step 16tThe specific calculation method comprises the following steps:
Tt＝KtT0
wherein K is a temperature decay parameter of the set simulated annealing, T0T is the number of repetitions of step 16 for a set simulated annealing initiation temperature, and is at least 1.
Further, the step 6 of sequentially discretizing the frequency distribution range specifically includes:
fis＝[fiminfimin+Δf fimin+2Δf fimin+3Δf ... fimax]T
wherein f isisFor the instantaneous frequency value of the discrete frequency vector, Δ f is the set frequency resolution, specifically:
wherein f isiminIs the minimum frequency, f, of the original polynomial phase signal s (t)imaxThe maximum frequency of the original polynomial phase signal s (t), K is the set discrete frequency dimension, and i is the imaginary unit.
Further, the specific calculation method for reserving the main lobe response for the frequency distribution range in step 6 is as follows:
wherein the former factor represents the direct frequency distribution of the Sinc function; the latter factor represents an operation that retains only the main lobe response, AcFor the c-th signal component hcThe intensity at the maximum of the spectral envelope of (a) is taken to be a complex value,i is the imaginary unit, Δ f is the set frequency resolution, fisAs a discrete frequency vector fsInstantaneous frequency value of fic(t) is the respective signal component hcCorresponding instantaneous frequency curve fc(t) instantaneous frequency value.
The present invention may be embodied in other specific forms without departing from the spirit or essential attributes thereof, and it should be understood that various changes and modifications can be effected therein by one skilled in the art without departing from the spirit and scope of the invention as defined in the appended claims.
Claims (10)
1. A polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing is characterized by comprising the following steps:
step 1: performing multiple decompositions on an original polynomial phase signal s (t) received by a radar by using a simulated annealing algorithm, wherein each decomposition determines an optimal model order and an optimal coefficient set to be determined of a signal component, and specifically:
step 11: generating a reference function h corresponding to the original polynomial phase signal s (t) according to the polynomial phase modelpAnd initializing a reference function hpModel order of (N)11, and calculating the energy E of the original polynomial phase signal s (t)0Randomly generating a reference function hpTwo sets of coefficients to be determined { a }n},n＝0,1,2,...,N1Wherein a set of coefficients { a } to be determinednTaking the solution as an initial solution in a simulated annealing algorithm;
step 12: two coefficient sets to be determined { a }nH corresponding topThe conjugate of the first and second signals is multiplied by an original polynomial phase signal s (t) respectively to obtain a mixed modulation signal x (t), Fourier transformation is carried out on the mixed modulation signal x (t) to obtain a transformation result X (f), a characteristic value of one frequency domain characteristic in the transformation result X (f) is calculated, the characteristic value is used as the fitness of the evaluation transformation result X (f), and the fitness is good or bad according to each transformation result X (f), wherein the transformation result X (f) with the optimal fitness corresponds to the current optimal solution Gbest, and the transformation result with the suboptimal fitness is a secondary optimal solution Prebast, and simultaneously assigning the current optimal solution Gbest to a preset;
step 13: according to the pre, according to the set annealing step factor TstepCalculating a new possible solution Nextp, and calculating the fitness fit of the possible solution Nextp by the method of step 12nextThe fitness fit of Nextp is solved if possiblenextReplacing the fitness of the current optimal solution Gbest which is superior to the fitness of the current optimal solution Gbest in the step 12 with the current optimal solution Gbest, and replacing the original optimal solution with a secondary optimal solution PreBest, otherwise, keeping the current optimal solution Gbest and the secondary optimal solution PreBest unchanged;
step 14: another new possible solution Nextq is randomly generated and its corresponding set of pending coefficients { a }nCalculating the fitness fit of the new possible solution Nextq corresponding to the transformation result X (f) according to the method of step 12newIf the fitness fit is fitnewBetter than fitness fitnextReplacing Nextq of the new possible solution with pre, otherwise according to the replacement probability pωThe magnitude of the random probability rand is used to decide whether to replace Nextq with Prelocal, if the probability p is replacedωIf the probability is greater than the random probability rand, replacing Nextq with the pre, otherwise, keeping the pre unchanged;
step 15, repeating the step 13 to the step 14 until the cycle number reaches L markov chain values, so as to obtain the current optimal solution Gbest and the suboptimal solution PreBest after iteration;
step 16: calculating the fitness difference value delta fit of the current optimal solution Gbest and the suboptimal solution PreBest obtained in the step 15, and calculating the annealing temperature TtWherein the annealing temperature TtDecays as the number of iterations of step 16 increases, if the absolute value of the fitness difference Δ fit is greater than the set annealing tolerance ξ, and the annealing temperature TtGreater than a set terminal temperature TlowThen steps 13-15 are repeated until the absolute value of the fitness difference Δ fit is not greater than the set annealing tolerance ξ, or the annealing temperature TtNot greater than a set cut-off temperature TlowSo as to obtain the final optimal solution Gtest';
and step 17: order of model N1Sequentially obtaining the maximum possible order _ max from 2, and repeating the steps 11 to 16 to obtain the final optimal solution Gbest' under different model orders;
step 18: selecting a solution with optimal fitness from all optimal solutions Gbest' obtained in the step 17 so as to determine the optimal model order N of the current signal component of the original polynomial phase signal s (t)pAnd NpCorresponding optimal pending coefficient set { an}max；
Step 2: using the optimal model order NpAnd an optimal pending coefficient set { a }n}maxMultiplying the determined conjugate of the signal component by the original polynomial phase signal s (t) and performing Fourier transform to obtain a frequency spectrum X '(f), and performing inverse Fourier transform after zeroing the intensity complex value at the maximum value of the envelope of the frequency spectrum X' (f) so as to obtain a time domain signal y (t);
and step 3: using model order of NpAnd an optimal pending coefficient set { a }n}maxMultiplying the determined signal component by a time domain signal y (t) to obtain a residual signal z (t), and finishing the decomposition;
and 4, step 4: calculating the energy E of the residual signal z (t)dWith the energy E of the original polynomial phase signal s (t)0Taking the ratio R, if the ratio R is less than the set threshold gamma or the decomposition times reaches the set upper limit number NmaxStopping decomposition to obtain the optimal model order N of each component signalpAnd an optimal pending coefficient set { a }n}max(ii) a Otherwise, replacing the original polynomial phase signal s (t) in the step 1 with the residual signal z (t) calculated in the step 3 to recalculate the mixed modulation signal x (t), and repeating the steps 1 to 3 until the ratio R is smaller than the set threshold gamma or the decomposition times reaches the set upper limit number Nmax；
And 5: the optimal model order NpAnd an optimal pending coefficient set { a }n}maxThe corresponding signal components are numbered h in the order of decompositioncWherein C is 1, 2.., C, wherein C is the number of adaptive decompositions;
step 6: obtaining each signal component h according to the physical definition of the instantaneous frequency of the single component signalcCorresponding instantaneous frequency curve fc(t) and from all signal components hcInstantaneous frequency curve fc(t) determining the frequency distribution range of the original polynomial phase signal s (t), and finally performing discretization, main lobe response reservation and superposition operations on the frequency distribution range in sequence to obtain the time-frequency joint distribution f (t) of the original polynomial phase signal s (t).
2. The simulated annealing-based polynomial phase signal adaptive time-frequency transform method as claimed in claim 1, wherein said original polynomial phase signal s (t) corresponds to a reference function hpThe method specifically comprises the following steps:
where order _ max is the maximum possible order of the original polynomial phase signal s (t).
3. The simulated annealing-based polynomial phase signal adaptive time-frequency transform method of claim 1, wherein the frequency distribution range of the original polynomial phase signal s (t) is specifically:
all signal components hcInstantaneous frequency curve fcThe instantaneous frequency minimum in (t) is the minimum frequency f of the original polynomial phase signal s (t)iminThe instantaneous frequency maximum is the maximum frequency f of the original polynomial phase signal s (t)imax。
4. The simulated annealing-based polynomial phase signal adaptive time-frequency transform method as claimed in claim 1, wherein said set of coefficients to be determined { a } isnIs specifically a vector
Wherein N is1Representing the model order of the current optimization, a0For a fixed value, set as the peak obtained by direct Fourier transform of the original polynomial phase signal s (t)The position serial number of the sampling point of the value;
the coefficient to be optimized is generated by a random function Rand with the value range of [ -M, M]Where M is the number of sampling points of the original polynomial phase signal s (t), specifically:
wherein rand (2, N)1) Representing random generation of 2 dimensions N1Distribution of (2) in [0,1 ]]To obtain a model order of N1The initial solution of (a).
5. The simulated annealing-based polynomial phase signal adaptive time-frequency transformation method as claimed in claim 1, wherein the frequency domain feature of step 12 is any one of a spectrum peak, a spectrum contrast or a spectrum shape entropy of a spectrum, wherein a feature value of each frequency domain feature is a spectrum peak XpSpectral contrast XcAnd the spectral shape entropy X of the frequency spectrumeRespectively as follows:
Xp＝max{abs[X(f)]}
Xe＝sum(-pS·logpS)
wherein
sum (-) represents the summation operation of the vectors, abs (-) is the absolute value;
when selecting spectral peak XpOr spectral contrast XcWhen the fitness is larger, the evaluation result is better, and the spectrum shape entropy X of the selected frequency spectrumeWhen the fitness is small, the evaluation is madeThe better the result.
6. The simulated annealing-based polynomial phase signal adaptive time-frequency transform method as claimed in claim 1 wherein step 13 is performed according to a set annealing step factor TstepCalculating a new possible solution Nextp, wherein the specific method is as follows:
Nextp＝Prelocal+Tstep·[2M·rand(1,N1)-M]
wherein rand (1, N)1) Representing random generation of 1 dimension N1Distribution of (2) in [0,1 ]]M is the number of samples of the original polynomial phase signal s (t).
7. The method for adaptive time-frequency transform of polynomial phase signals based on simulated annealing as claimed in claim 1 wherein said alternative probability p of step 14ωThe specific calculation method comprises the following steps:
pω＝exp[-(fitnext-fitnew)]。
8. the simulated annealing-based polynomial phase signal adaptive time-frequency transform method of claim 1, wherein said annealing temperature T of step 16tThe specific calculation method comprises the following steps:
Tt＝KtT0
wherein K is a temperature decay parameter of the set simulated annealing, T0T is the number of repetitions of step 16 for a set simulated annealing initiation temperature, and is at least 1.
9. The simulated annealing-based polynomial phase signal adaptive time-frequency transformation method as claimed in claim 1, wherein said discretizing the frequency distribution range in sequence in step 6 specifically comprises:
fis＝[fiminfimin+Δf fimin+2Δf fimin+3Δf...fimax]T
wherein f isisFor the instantaneous frequency value of the discrete frequency vector, Δ f is the set frequency resolution, specifically:
wherein f isiminIs the minimum frequency, f, of the original polynomial phase signal s (t)imaxThe maximum frequency of the original polynomial phase signal s (t), K is the set discrete frequency dimension, and i is the imaginary unit.
10. The adaptive time-frequency transform method for polynomial phase signals based on simulated annealing according to claim 1, wherein the specific calculation method for preserving the main lobe response of the frequency distribution range in step 6 is:
wherein A iscFor the c-th signal component hcI is the imaginary unit, deltaf is the set frequency resolution, f is the intensity complex value at the maximum of the spectral envelope of (1)isAs a discrete frequency vector fsInstantaneous frequency value of fic(t) is the respective signal component hcCorresponding instantaneous frequency curve fc(t) instantaneous frequency value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710914700.5A CN107622035B (en) | 2017-09-30 | 2017-09-30 | Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710914700.5A CN107622035B (en) | 2017-09-30 | 2017-09-30 | Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107622035A CN107622035A (en) | 2018-01-23 |
| CN107622035B true CN107622035B (en) | 2020-07-17 |
Family
ID=61091513
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710914700.5A Active CN107622035B (en) | 2017-09-30 | 2017-09-30 | Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107622035B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112130138B (en) * | 2020-07-06 | 2022-10-14 | 电子科技大学 | Ionized layer phase pollution correction method for over-the-horizon radar |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102023963A (en) * | 2010-12-13 | 2011-04-20 | 西安电子科技大学 | High-speed multi-mode time domain and frequency domain transform method |
| CN104731762A (en) * | 2015-03-25 | 2015-06-24 | 西安电子科技大学 | Cubic phase signal parameter estimation method based on cyclic shift |
| CN105606894A (en) * | 2016-01-28 | 2016-05-25 | 南京信息工程大学 | Instantaneous frequency estimation method based on simulated annealing (SA) |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9454511B2 (en) * | 2011-05-04 | 2016-09-27 | American University | Windowing methods and systems for use in time-frequency analysis |
-
2017
- 2017-09-30 CN CN201710914700.5A patent/CN107622035B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102023963A (en) * | 2010-12-13 | 2011-04-20 | 西安电子科技大学 | High-speed multi-mode time domain and frequency domain transform method |
| CN104731762A (en) * | 2015-03-25 | 2015-06-24 | 西安电子科技大学 | Cubic phase signal parameter estimation method based on cyclic shift |
| CN105606894A (en) * | 2016-01-28 | 2016-05-25 | 南京信息工程大学 | Instantaneous frequency estimation method based on simulated annealing (SA) |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107622035A (en) | 2018-01-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107729289B (en) | Polynomial phase signal self-adaptive time-frequency transformation method based on genetic optimization | |
| Gokhale et al. | Time domain signal analysis using wavelet packet decomposition approach | |
| Venema et al. | A stochastic iterative amplitude adjusted Fourier transform algorithm with improved accuracy | |
| CN110705762B (en) | Ubiquitous power Internet of things perception data missing repairing method based on matrix filling | |
| CN111723701B (en) | Underwater target identification method | |
| CN113094993B (en) | Modulation signal denoising method based on self-coding neural network | |
| CN107729288B (en) | Polynomial phase signal time-frequency transformation method based on particle swarm optimization | |
| CN109086555A (en) | Using the simulation spectrum curve emulation mode of Monte Carlo method | |
| CN107622036B (en) | Polynomial phase signal self-adaptive time-frequency transformation method based on ant colony optimization | |
| CN107783939B (en) | Model-driven polynomial phase signal self-adaptive time-frequency transformation method | |
| CN107622035B (en) | Polynomial phase signal self-adaptive time-frequency transformation method based on simulated annealing | |
| KR102645315B1 (en) | Instrument synthesis sound generation system | |
| CN110428848B (en) | Speech enhancement method based on public space speech model prediction | |
| Popescu et al. | New procedure for change detection operating on Rényi entropy with application in seismic signals processing | |
| CN113935246A (en) | Signal robust sparse time-frequency analysis method, terminal equipment and storage medium | |
| US20070124137A1 (en) | Highly optimized nonlinear least squares method for sinusoidal sound modelling | |
| CN116522073A (en) | Adaptive signal decomposition method of fractional Fourier transform domain, storage medium and computer | |
| Bagwe et al. | Automated Radar Signal Analysis Based on Deep Learning | |
| CN113655534B (en) | Nuclear magnetic resonance FID signal noise suppression method based on multi-linear singular value tensor decomposition | |
| CN115859048A (en) | Noise processing method and device for partial discharge signal | |
| Sturm et al. | Dark energy in sparse atomic estimations | |
| CN115910091A (en) | Method and device for separating generated voice by introducing fundamental frequency clues | |
| Zivanovic et al. | Nonparametric spectrum interpolation methods: A comparative study | |
| CN112054803A (en) | Communication signal sorting method based on compressed sensing | |
| de Oliveira et al. | Implementation of an electrical signal compression system using sparce representation |
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 |