An adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion

Abstract

The application discloses a kind of based on kurtosis maximization and continuity criterion's adaptive time-frequency analysis method, comprising:1, obtaining the pulse signal sample data sequence to be processed;2, initialization adaptive time-frequency analysis initial window length set and time frame parameter;3, according to initial window length set and time frame parameter update adaptive window length set 4, calculate the kurtosis corresponding to each window length in current time frame adaptive window length set, obtain the initial adaptive time-frequency analysis window length of current time frame;5. According to window length change continuity criterion, the initial adaptive time-frequency analysis window length is corrected, and the best time-frequency analysis window length of current time frame is obtained;6. Obtain the adaptive short-time power spectrum and instantaneous discretization frequency of current time frame, and obtain the contraction window length set according to the change trend of instantaneous discretization frequency;7. Determine whether to process all time frames, if not, return to step 3, otherwise, go to step 8;8. According to the best time-frequency analysis window length sequence reconstructs adaptive time-frequency analysis result.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of signal processing, and particularly relates to an adaptive time-frequency analysis method based on kurtosis maximization and continuity characteristics. BACKGROUND

[0002] Since the spectrum of a non-stationary signal is mostly time-varying, the classic Fourier transform cannot fully reflect the change of the frequency component of the signal over time. Therefore, in order to reveal the global information of the signal in the time and frequency dimensions, joint processing in the time-frequency domain is required to obtain more information, and thus the time-frequency analysis is proposed to solve the analysis problem of the non-stationary signal.

[0003] Scholars in the field of signal processing usually divide the time-frequency analysis into linear time-frequency analysis and quadratic time-frequency analysis. The representative classical algorithm of the quadratic time-frequency analysis is the Wigner-Ville distribution (WVD), which has high-resolution time-frequency analysis performance, but for multi-component signals or non-linear frequency modulation signals, the algorithm and its subsequent improved algorithms are plagued by cross terms, and the computational amount is significantly higher than that of the linear time-frequency analysis method, so it is difficult to be deployed and applied in real-time processing systems. The representative classical algorithm of the linear time-frequency analysis is the short-time Fourier transform (STFT), which is simple to calculate, has clear physical meaning, and has no cross terms, so it is widely used in the processing of radar or sonar signals. The problem is that the resolution of the time-frequency analysis result is limited by the uncertainty principle, and it is impossible to improve the time and frequency domain resolution at the same time, so the method of selecting the length of the short-time analysis window for trade-off between time and frequency domain resolution is the core of the adaptive STFT method. In general, the length of the short-time analysis window of the STFT is selected depending on the experience-based debugging.

[0004] In order to more flexibly design the time-frequency analysis method, the adaptive linear time-frequency analysis method is proposed to solve the problem of selecting the length of the short-time analysis window. The existing adaptive time-frequency analysis method to some extent overcomes the shortcoming that the length of the window of the STFT does not change with the characteristics of the signal, and effectively trades off the time and frequency resolution, but the redundancy and computational amount are significantly improved. SUMMARY

[0005] The purpose of the application is to provide an adaptive time-frequency analysis method based on kurtosis maximization and continuity characteristics, so as to solve the technical problem that the conventional time-frequency analysis method of the pulse signal is difficult to realize the trade-off between the time and frequency domain resolution, the adaptive time-frequency analysis method has redundancy and large computational amount, and is difficult to be deployed and used in the real-time signal processing system.

[0006] To solve the above technical problems, the application provides an adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria, which comprises the following steps:

[0007] Step 1, obtaining a signal sampling data sequence x(n) to be processed, n = 1, 2, …, N, n is the sampling point index of x(n), N = 2 a , a is an integer greater than 10;

[0008] Step 2, initializing the initial window length set Ω of adaptive time-frequency analysis and the time frame parameter;

[0009] Step 3, sliding the time frame of the signal and updating the adaptive window length set Ω according to the initial window length set Ω and the time frame parameter l ;

[0010] Step 4, calculating the kurtosis h(l, m) corresponding to each window length in the adaptive window length set Ω l of the current time frame, and obtaining the initial adaptive time-frequency analysis window length W adp (l) of the current time frame according to the maximum kurtosis criterion;

[0011] Step 5, correcting the initial adaptive time-frequency analysis window length W adp (l) of the current time frame according to the continuity criterion of window length change, and obtaining the optimal time-frequency analysis window length W opt (l) of the current time frame;

[0012] Step 6, obtaining the adaptive short-time power spectrum Y opt (l, k opt ) and the instantaneous discrete frequency of the current time frame according to the optimal time-frequency analysis window length W l (l) of the current time frame, and obtaining the contraction window length set according to the change trend of the instantaneous discrete frequency

[0013] Step 7, judging whether all time frames are processed, if not, returning to step 3, otherwise going to step 8;

[0014] Step 8, reconstructing the adaptive time-frequency analysis result according to the optimal time-frequency analysis window length W opt (l), l = 1, 2, …, L; Wherein, L is the total number of time frames.

[0015] Further, the specific method of step 1 is as follows:

[0016] Receiving real-time collection data of N sampling points from a sensor as the data sequence x(n) to be processed, or extracting N sampling point data containing the entire pulse signal to be analyzed from a memory as the data sequence x(n) to be processed, N = 2 a , a is an integer greater than 10.

[0017] Further, the specific method of step 2 is as follows:

[0018] The time-frequency parameters required for adaptive time-frequency analysis of the pulse signal to be analyzed are initialized, including the initialization of the following parameters:

[0019] The minimum initial window length H of the time-frequency analysis of the pulse signal to be analyzed min is initialized as: 4≤H min ≤16, and is an even number;

[0020] The maximum initial window length H of the time-frequency analysis of the pulse signal to be analyzed max is initialized as: H max =N;

[0021] The window length search step S of the time-frequency analysis of the pulse signal to be analyzed W is initialized as: 1<S W <5, and is an even number, and satisfies (H max -H min ) / S W is an integer;

[0022] The window length change rate threshold r of the time-frequency analysis of the pulse signal to be analyzed a is initialized as: 0.25<r a <0.45, and is a real number;

[0023] The minimum initial time frame sequence index l of the adaptive window length set starting to shrink begin is initialized as: N / 8<l begin <N / 4, and l begin is an even number;

[0024] Let the initial window length set Ω of the time-frequency analysis of the pulse signal to be analyzed be:

[0025] Ω={H min ,H min +S W ,H min +2S W ,...,H max};

[0026] The sliding window step S is initialized as S=int[H min / 4], and the total time frame number L of the adaptive time-frequency analysis of the pulse signal to be analyzed is initialized as: L=int[N / S]-1, wherein int[] represents the downward integer of [] inside;

[0027] The time frame sequence index l of the adaptive time-frequency analysis of the pulse signal to be analyzed is initialized as: l=0.

[0028] Further, the specific method of step 3 is as follows:

[0029] Step 3.1, slide the time frame of the signal, i.e. let the time frame sequence index l = l + 1;

[0030] Step 3.2, update the adaptive window length set Ω l to:

[0031]

[0032] update the adaptive window length set Ω l of the lth time frame based on the above formula to:

[0033] Ω l = {W min (l), W min (l) + S W , W min (l) + 2S W ,..., W max (l)}

[0034] wherein W max (l) is the maximum window length of the adaptive window length set Ω l of the lth time frame, and W min (l) is the minimum window length of the adaptive window length set Ω l of the lth time frame.

[0035] Further, the specific method of Step 4 is as follows:

[0036] Step 4.1, count the number of window length elements in the adaptive window length set Ω l of the lth time frame, denoted as M l ;

[0037] Step 4.2, define m as the window length index within the adaptive window length set Ω l of the lth time frame, and initialize m = 0;

[0038] Step 4.3, update the window length index within the adaptive window length set Ω l of the lth time frame, i.e. let m = m + 1, and take the mth window length W(l, m) within the adaptive window length set Ω l of the lth time frame:

[0039] W(l, m) = W min (l) + (m - 1) · S W ;

[0040] Step 4.4, calculate the time-varying spectrum y m (l, k l ) corresponding to the mth window length of the lth time frame:

[0041]

[0042] wherein kl =0,1,...,W max (l) / 2-1 represents the discretized frequency index of the frame at time l, where j represents the imaginary unit, i.e.

[0043] The value of n is an integer satisfying n1(l,m)≤n≤n2(l,m), where the starting position n1(l,m) and ending position n2(l,m) of the sliding window are respectively:

[0044] n1(l,m)=max{1,int[(l-1)S+1-W(l,m) / 2]}

[0045] n2(l,m)=min{N,int[(l-1)S+1+W(l,m) / 2]}

[0046] Where min{·} represents taking the minimum value of the two, and max{·} represents taking the maximum value of the two;

[0047] Step 4.5: Calculate the short-time power spectrum Y corresponding to the m-th window length of the current time frame l. m (l,k l ):

[0048]

[0049] Where | represents the modulo operation;

[0050] Step 4.6: Calculate the short-time power spectrum Y corresponding to the m-th window length of the current time frame l. m (l,k l Normalization is performed to obtain the normalized short-time power spectrum corresponding to the m-th window length of the l-th frame.

[0051]

[0052] in, Indicates at the discretized frequency k l The maximum value of this dimension;

[0053] Step 4.7: Based on the normalized short-time power spectrum corresponding to the m-th window length of the current l-th time frame. The kurtosis h(l,m) corresponding to the m-th window length in the l-th frame is calculated as follows:

[0054]

[0055] in, The normalized short-time power spectrum corresponding to the m-th window length in the l-th frame at the discretized frequency k l The mean in this dimension:

[0056]

[0057] Step 4.8, judge whether m≥M is true, if true, go to step 4.9, otherwise, return to step 4.3; l

[0058] Step 4.9, according to the M l peak degrees h(l,m) corresponding to the M l window length of the current l time frame, obtain the initial adaptive time-frequency analysis window length W adp (l) of the current l time frame:

[0059]

[0060] wherein, W(l,m) is the maximum value in the parentheses.

[0061] Further, the specific method of step 5 is as follows:

[0062] Step 5.1, judge whether l=1 is true, if true, let the l time frame best time-frequency analysis window length equal to the l time frame initial adaptive time-frequency analysis window length, that is, W opt (l) = W adp (l), go to step 6, otherwise, go to step 5.2;

[0063] Step 5.2, calculate the l time frame best time-frequency analysis window length W opt (l):

[0064]

[0065] wherein, r(l) represents the window length change rate of the l time frame:

[0066]

[0067] Further, the specific method of step 6 is as follows:

[0068] Step 6.1, according to the l time frame best time-frequency analysis window length W opt (l), obtain the adaptive time-varying spectrum y opt (l,k l )

[0069]

[0070] wherein, the value range of n satisfies n'1(l,m)≤n≤n'2(l,m) integer, n'1(l,m) and n'2(l,m) respectively represent according to the l time frame best time-frequency analysis window length W opt ​(l) The starting and ending positions of the obtained sliding window are calculated as follows:

[0071] n'1(l,m)=max{1,int[(l-1)S+1-W opt (l) / 2]}

[0072] n'2(l,m)=min{N,int[(l-1)S+1+W opt (l) / 2]}

[0073] Step 6.2: Based on the adaptive time-varying spectrum y of the l-th time frame opt (l,k l ), calculate the adaptive short-time power spectrum Y of the l-th frame. opt (l,k l ):

[0074]

[0075] Step 6.3: Based on the adaptive short-time power spectrum Y of the l-th time frame opt (l,k l Estimate the instantaneous discretization frequency of the signal to be analyzed in the l-th time frame.

[0076]

[0077] in, This represents the value of k that makes the expression inside the parentheses reach its maximum value. l value;

[0078] Step 6.4: Determine if l ≥ l begin Whether it is true, that is, whether l has been accumulated. begin If the information of each time frame is true, proceed to step 6.5; otherwise, proceed to step 6.9.

[0079] Step 6.5: Initialize the directional trend factor α as follows: 0.6 < α < 0.8, where α is a real number;

[0080] Step 6.6: Calculate the instantaneous discretized frequency using the Cox-Stuart trend test algorithm. The changing trend, with the following sub-steps:

[0081] Step 6.6.1: Starting from the l-th time frame, take l frames backwards. begin The instantaneous discretized frequencies of each time frame, including the l-th time frame, constitute the dataset.

[0082]

[0083] Pick and Form a set of data Where i represents the trend test data index. The left value of the data pair is used to test the trend. For trend testing, the right value i = ll begin +1,ll begin +2,…,lc,i+c=l-c+1,l-c+2,…,l, where c=int[l begin [ / 2] indicates the number of trend test data pairs;

[0084] Step 6.6.2: Subtract the two numbers in each data pair along the time series direction and record the trend value. And introduce the number of positive trends S + Record D i The number of positive trends, S - Record D i The quantity is negative;

[0085] Step 6.6.3: Based on the number S of positive trends + and the number of negative trends S - Calculate the test result T:

[0086]

[0087] Step 6.7: If the test result T is 0, proceed to step 7; otherwise, proceed to step 6.8.

[0088] Step 6.8: Based on the trend of the l-th frame, examine T and the previous l of the l-th frame. begin -1 time frame optimal time-frequency analysis window length W opt (l), obtain the long set of the contracted window for the l-th time frame. The expression is:

[0089]

[0090] Among them, W' min (l) and W' max (l) represents the set of the contracted window length of the l-th time frame. The minimum and maximum window lengths are calculated as follows:

[0091]

[0092] W mid (l) represents the first l frames of the l-th time frame participating in the Cox-Stuart trend test. begin -1 time frame median of the optimal time-frequency analysis window length, W mid (l) is calculated as follows:

[0093] Wmid (l) = median{W opt (l-l begin +1), W opt (l-l begin +2),..., W opt (l-1)}

[0094] where median{•} denotes the median of the values in {•}, and if the length of {•} is odd, the median is taken, and if the length of {•} is even, the average of the two middle values is taken, since l begin is initialized to even in step 2, l begin -1 is odd, W mid (l) is the median of the sequence;

[0095] Step 6.9, let the set of adaptive window lengths of the l-th time frame equal to the set of contracted window lengths of the l-th time frame, i.e.

[0096] Further, the specific method of step 7 is as follows:

[0097] If l < L is true, return to step 3.1, otherwise, go to step 8.

[0098] Further, the specific method of step 8 is as follows:

[0099] Step 8.1, normalize the length L norm of the spectrum to:

[0100]

[0101] where denotes the maximum value in the dimension of the time frame sequence index l;

[0102] and reinitialize the time frame sequence index l = 0;

[0103] Step 8.2, slide the time frame of the signal, i.e. let the time frame sequence index l = l + 1, and according to the optimal time-frequency analysis window length W opt (l) of the l-th time frame, calculate the normalized length adaptive time-varying spectrum y norm (l, k) of the l-th time frame:

[0104]

[0105] where n takes an integer value satisfying n'1(l, m)≤n≤n'2(l, m), and n'1(l, m) and n'2(l, m) are defined as in step 6.1; k denotes the discrete frequency index obtained according to the normalized length L norm of the spectrum;

[0106] Step 8.3, normalizing the spectrum length of the Ith frame of the short-time power spectrum Y norm (l,k) is calculated as:

[0107]

[0108] Step 8.4, judging whether l=L is true, if true, going to step 8.5; otherwise, returning to step 8.2;

[0109] Step 8.5, arranging the L frames of the normalized spectrum length of the short-time power spectrum Y norm (l,k) in time sequence, performing amplitude normalization on the L frames of the normalized spectrum length of the short-time power spectrum Y norm (l,k) to obtain the final time-frequency analysis result

[0110]

[0111] wherein, represents the maximum value in the two dimensions of l and k.

[0112] Beneficial effects: compared with the prior art, the technical scheme of the present application has the following beneficial technical effects:

[0113] 1. In the time-frequency analysis method of the present application, the frequency domain narrowband characteristics of the local signal are utilized when the short-time Fourier transform of the signal is obtained, the statistical value describing the time-frequency energy aggregation of the signal, i.e. the spectral kurtosis, is obtained, and the length of the short-time analysis window is adaptively selected to perform the short-time Fourier transform on the local signal. In steps 2 and 3, the high-energy aggregation spectrum of the signal is adaptively obtained through the normalized spectral kurtosis of each window length in the window length set.

[0114] 2. In the time-frequency analysis method of the present application, the internal relationship between the signal analysis window length and the instantaneous discretized frequency of the signal to be analyzed is utilized when the adaptive short-time analysis window length and the short-time Fourier transform of the signal are obtained. In step 4, the continuity criterion of the signal is used to correct the adaptive short-time analysis window length sequence by calculating the change rate of the signal adaptive window length. The change rate of the window length compared with the previous time frame is calculated through the obtained signal adaptive window length, and then compared with the preset change rate threshold, and the abnormal analysis window length that may exist is corrected. The short-time Fourier transform result of the signal that is more ideal can be calculated at a lower signal-to-noise ratio.

[0115] 3、The time-frequency analysis method of the present application uses the characteristics of common pulse signals when calculating the adaptive short-time analysis window length of each time frame of the signal and the short-time Fourier transform result, judges the change trend of the instantaneous discrete frequency of the signal based on the historical information of the adaptive short-time analysis window length, thereby realizing the distribution range of the self-shrinking window length set and reducing the redundancy and computational complexity of the method of the present application. In step 5, the normalized short-time power spectrum Y opt (l,k l ) is used to estimate the instantaneous discrete frequency of the signal, in step 6, the Cox-stuart trend test algorithm is used to test the increasing or decreasing trend of the signal, so as to guide the shrinkage range and shrinkage direction of the adaptive time-frequency analysis window length set, and in step 7, the ideal time-frequency analysis result Y norm (l,k) is finally obtained. Through the adaptive shrinkage of the window length set, the redundancy and large amount of calculation of the adaptive time-frequency analysis are suppressed, and the robustness of the time-frequency analysis result under low signal-to-noise ratio is improved. The present application fully utilizes the short-time narrowband characteristics of common pulse signals, and can obtain adaptive time-frequency analysis results with high time-frequency focusing stability. BRIEF DESCRIPTION OF DRAWINGS

[0116] Figure 1 Fig. 1 is a flowchart of the adaptive time-frequency analysis method based on the kurtosis maximization and continuity criterion of the present application;

[0117] Figure 2 Fig. 2 is the adaptive time-frequency analysis window length sequence obtained in Example 1;

[0118] Figure 3 Fig. 3 is a normalized adaptive time-frequency analysis result graph of the simulation pulse signal in Example 1;

[0119] Figure 4 Fig. 4 is the adaptive time-frequency analysis window length sequence obtained in Example 2;

[0120] Figure 5 Fig. 5 is a normalized adaptive time-frequency analysis result graph of the simulation pulse signal in Example 2; DETAILED DESCRIPTION

[0121] In order to better understand the purpose, structure and function of the present application, the adaptive time-frequency analysis method based on the kurtosis maximization and continuity criterion of the present application will be further described in detail below in combination with the drawings.

[0122] As shown in Figure 1 , the present application proposes an adaptive time-frequency analysis method based on the kurtosis maximization and continuity criterion, which comprises the following steps:

[0123] Step 1, obtaining the signal sample data sequence x(n) to be processed, n=1, 2, …, N, n is the sampling point index of x(n), and N=2a , a is an integer greater than 10;

[0124] Step 2, initializing the initial window length set Ω of adaptive time-frequency analysis and time frame parameters;

[0125] Step 3, sliding the time frame of the signal and updating the adaptive window length set Ω according to the initial window length set Ω and the time frame parameters l ;

[0126] Step 4, calculating the kurtosis h(l, m) corresponding to each window length in the adaptive window length set Ω of the current time frame, and obtaining the initial adaptive time-frequency analysis window length W l (l) of the current time frame according to the maximum kurtosis criterion adp ;

[0127] Step 5, correcting the initial adaptive time-frequency analysis window length W adp (l) of the current time frame according to the continuity criterion of window length change, and obtaining the optimal time-frequency analysis window length W opt (l) of the current time frame;

[0128] Step 6, obtaining the adaptive short-time power spectrum Y opt (l, k opt ) and the instantaneous discrete frequency of the current time frame according to the optimal time-frequency analysis window length W l (l) of the current time frame, and obtaining the contraction window length set according to the change trend of the instantaneous discrete frequency

[0129] Step 7, judging whether all time frames have been processed, if not, returning to Step 3, otherwise going to Step 8;

[0130] Step 8, reconstructing the adaptive time-frequency analysis result according to the optimal time-frequency analysis window length W opt (l), l = 1, 2, …, L , where L is the total number of time frames.

[0131] Further, the specific method of Step 1 is as follows:

[0132] Receiving real-time collection data of N sampling points from the sensor as the data sequence x(n) to be processed, or extracting N sampling point data containing the entire pulse signal to be analyzed from the memory as the data sequence x(n) to be processed, N = 2 a , a is an integer greater than 10.

[0133] Further, the specific method of Step 2 is as follows:

[0134] Initialize time-frequency parameters required for adaptive time-frequency analysis of the pulse signal to be analyzed, specifically including the initialization of the following parameters:

[0135] Minimum initial window length H of time-frequency analysis of the pulse signal to be analyzed min initialized as: 4≤H min ≤16, and is an even number;

[0136] Maximum initial window length H of time-frequency analysis of the pulse signal to be analyzed max initialized as: H max =N;

[0137] Window length search step S of time-frequency analysis of the pulse signal to be analyzed W initialized as: 1<S W <5, and is an even number, and satisfies (H max -H min ) / S W is an integer;

[0138] Window length change rate threshold r of time-frequency analysis of the pulse signal to be analyzed a initialized as: 0.25<r a <0.45, and is a real number;

[0139] Minimum initial time frame sequence index l of the start of the adaptive window length set contraction begin initialized as: N / 8<l begin <N / 4, and l begin is an even number;

[0140] Let the initial window length set Ω of time-frequency analysis of the pulse signal to be analyzed be:

[0141] Ω={H min ,H min +S W ,H min +2S W ,...,H max};

[0142] Initialize the sliding window step S as S=int[H min / 4], and then initialize the total time frame number L of adaptive time-frequency analysis of the pulse signal to be analyzed as L=int[N / S]-1, where int[] represents the downward integer of [] inside;

[0143] Initialize the time frame sequence index l of adaptive time-frequency analysis of the pulse signal to be analyzed as l=0.

[0144] Further, the specific method of step 3 is as follows:

[0145] Step 3.1, slide the time frame of the signal, i.e. let the time frame sequence index l=l+1;

[0146] Step 3.2, adaptive window length set Ω of the lth time frame l is updated as:

[0147]

[0148] The updated adaptive window length set Ω of the lth time frame based on the above formula is: l is:

[0149] Ω l = {W min (l), W min (l) + S W , W min (l) + 2S W ,..., W max (l)}

[0150] wherein W max (l) is the maximum window length of the adaptive window length set Ω l of the lth time frame, and W min (l) is the minimum window length of the adaptive window length set Ω l of the lth time frame.

[0151] Further, the specific method of step 4 is as follows:

[0152] Step 4.1, count the number of window length elements in the adaptive window length set Ω l of the lth time frame, denoted as M l ;

[0153] Step 4.2, define m as the window length index in the adaptive window length set Ω l of the lth time frame, and initialize m = 0;

[0154] Step 4.3, update the window length index in the adaptive window length set Ω l of the lth time frame, i.e. m = m + 1, and take the mth window length W(l, m) in the adaptive window length set Ω l of the lth time frame:

[0155] W(l, m) = W min (l) + (m - 1) · S W ;

[0156] Step 4.4, calculate the time-varying spectrum y m (l, k l ) corresponding to the mth window length of the lth time frame:

[0157]

[0158] wherein k l = 0, 1,..., Wmax (l) / 2-1, represents the discrete frequency index of the l-th time frame, j represents the imaginary unit, i.e.

[0159] The value range of n satisfies n1(l,m)≤n≤n2(l,m) integer, wherein the starting position n1(l,m) and the ending position n2(l,m) of the sliding window are respectively:

[0160] n1(l,m)=max{1,int[(l-1)S+1-W(l,m) / 2]}

[0161] n2(l,m)=min{N,int[(l-1)S+1+W(l,m) / 2]}

[0162] Wherein, min{·} represents taking the minimum value of two, max{·} represents taking the maximum value of two;

[0163] Step 4.5, calculating the short-time power spectrum Y m (l,k l ) corresponding to the m-th window length of the l-th time frame:

[0164]

[0165] Wherein, | | represents the modulus operation;

[0166] Step 4.6, normalizing the short-time power spectrum Y m (l,k l ) corresponding to the m-th window length of the l-th time frame, to obtain the normalized short-time power spectrum

[0167]

[0168] Wherein, represents the maximum value in the dimension of the discrete frequency k l ;

[0169] Step 4.7, according to the normalized short-time power spectrum (l,k l ) corresponding to the m-th window length of the l-th time frame, the kurtosis h(l,m) corresponding to the m-th window length of the l-th time frame is calculated:

[0170]

[0171] Wherein, is the mean value of the normalized short-time power spectrum l (l,k l ) corresponding to the m-th window length of the l-th time frame in the dimension of the discrete frequency k l ;

[0172]

[0173] Step 4.8, judge whether m≥M l is true, if true, go to step 4.9, otherwise, return to step 4.3;

[0174] Step 4.9, according to the M l peak degrees h(l,m) corresponding to the M l window length of the current l time frame, obtain the initial adaptive time-frequency analysis window length W adp (l) of the current l time frame:

[0175]

[0176] wherein, W(l,m) is the maximum value in the parentheses.

[0177] Further, the specific method of step 5 is as follows:

[0178] Step 5.1, judge whether l=1 is true, if true, let the l time frame optimal time-frequency analysis window length equal to the l time frame initial adaptive time-frequency analysis window length, that is, W opt (l)=W adp (l), go to step 6, otherwise, go to step 5.2;

[0179] Step 5.2, calculate the l time frame optimal time-frequency analysis window length W opt (l):

[0180]

[0181] wherein, r(l) represents the window length change rate of the l time frame:

[0182]

[0183] Further, the specific method of step 6 is as follows:

[0184] Step 6.1, according to the l time frame optimal time-frequency analysis window length W opt (l), obtain the adaptive time-varying spectrum y opt (l,k l )

[0185]

[0186] wherein, the value range of n satisfies n'1(l,m)≤n≤n'2(l,m) integer, n'1(l,m) and n'2(l,m) respectively represent the l time frame optimal time-frequency analysis window length W opt(l) The starting and ending positions of the obtained sliding window are calculated as follows:

[0187] n'1(l,m)=max{1,int[(l-1)S+1-W opt (l) / 2]}

[0188] n'2(l,m)=min{N,int[(l-1)S+1+W opt (l) / 2]}

[0189] Step 6.2: Based on the adaptive time-varying spectrum y of the l-th time frame opt (l,k l ), calculate the adaptive short-time power spectrum Y of the l-th frame. opt (l,k l ):

[0190]

[0191] Step 6.3: Based on the adaptive short-time power spectrum Y of the l-th time frame opt (l,k l Estimate the instantaneous discretization frequency of the signal to be analyzed in the l-th time frame.

[0192]

[0193] in, This represents the value of k that makes the expression inside the parentheses reach its maximum value. l value;

[0194] Step 6.4: Determine if l ≥ l begin Whether it is true, that is, whether l has been accumulated. begin If the information of each time frame is true, proceed to step 6.5; otherwise, proceed to step 6.9.

[0195] Step 6.5: Initialize the directional trend factor α as follows: 0.6 < α < 0.8, where α is a real number;

[0196] Step 6.6: Calculate the instantaneous discretized frequency using the Cox-Stuart trend test algorithm. The changing trend, with the following sub-steps:

[0197] Step 6.6.1: Starting from the l-th time frame, take l frames backwards. begin The instantaneous discretized frequencies of each time frame, including the l-th time frame, constitute the dataset.

[0198]

[0199] Pick and Form a set of data Where i represents the trend test data index. The left value of the data pair is used to test the trend. For trend testing, the right value i = ll begin +1,ll begin +2,…,lc,i+c=l-c+1,l-c+2,…,l, where c=int[l begin [ / 2] indicates the number of trend test data pairs;

[0200] Step 6.6.2: Subtract the two numbers in each data pair along the time series direction and record the trend value. And introduce the number of positive trends S + Record D i The number of positive trends, S - Record D i The quantity is negative;

[0201] Step 6.6.3: Based on the number S of positive trends + and the number of negative trends S - Calculate the test result T:

[0202]

[0203] Step 6.7: If the test result T is 0, proceed to step 7; otherwise, proceed to step 6.8.

[0204] Step 6.8: Based on the trend of the l-th frame, check T and the previous l of the l-th frame. begin -1 time frame optimal time-frequency analysis window length W opt (l), obtain the long set of the contracted window for the l-th time frame. The expression is:

[0205]

[0206] Among them, W' min (l) and W' max (l) represents the set of the contracted window length of the l-th frame. The minimum and maximum window lengths are calculated as follows:

[0207]

[0208] W mid (l) represents the first l frames of the l-th time frame participating in the Cox-Stuart trend test. begin -1 time frame median of the optimal time-frequency analysis window length, W mid (l) is calculated as follows:

[0209] Wmid (l) = median{W opt (l-l begin +1), W opt (l-l begin +2),..., W opt (l-1)}

[0210] where median{•} denotes taking the median of the values in {•}, and if the length of {•} is odd, the median is taken, and if the length of {•} is even, the average of the two middle values is taken, since l begin is initialized to even in step 2, l begin -1 is odd, W mid (l) is the median of the sequence;

[0211] Step 6.9, let the set of adaptive window lengths of the l-th time frame equal to the set of contracted window lengths of the l-th time frame, i.e.

[0212] Further, the specific method of step 7 is as follows:

[0213] If l < L is true, return to step 3.1, otherwise, go to step 8.

[0214] Further, the specific method of step 8 is as follows:

[0215] Step 8.1, normalize the length L norm of the spectrum to:

[0216]

[0217] where denotes the maximum value in the dimension of the time frame sequence index l;

[0218] and reinitialize the time frame sequence index l = 0;

[0219] Step 8.2, slide the time frame of the signal, i.e. let the time frame sequence index l = l + 1, and according to the optimal time-frequency analysis window length W opt (l) of the l-th time frame, calculate the normalized length adaptive time-varying spectrum y norm (l, k) of the l-th time frame:

[0220]

[0221] where n takes an integer value satisfying n'1(l, m)≤n≤n'2(l, m), and n'1(l, m) and n'2(l, m) are defined as in step 6.1; k denotes the discrete frequency index obtained according to the normalized length L norm of the spectrum;

[0222] Step 8.3, Normalized spectral length adaptive short-time power spectrum Y of the l-th frame norm (l,k) is calculated as:

[0223]

[0224] Step 8.4, judging whether l=L is true, if true, going to step 8.5; otherwise, returning to step 8.2;

[0225] Step 8.5, amplitude-normalizing the normalized spectral length adaptive short-time power spectrum Y of the L frames norm (l,k) in time sequence, amplitude-normalizing the normalized spectral length adaptive short-time power spectrum Y of the L frames norm (l,k) to obtain the final time-frequency analysis result

[0226]

[0227] wherein, represents the maximum value in the two dimensions of l and k.

[0228] In the embodiment of the present application, the simulation received pulse signal model is:

[0229] x(t) = s(t) + w(t)

[0230] wherein, w(t) is the background noise of the observation signal, which is Gaussian white noise; s(t) is the pulse signal, which is expressed as:

[0231]

[0232] wherein, A is the amplitude of the signal s(t), Π(t) is the duration of s(t), is the instantaneous phase of s(t), and Π(t) is expressed as:

[0233]

[0234] wherein, t0 and τ are the starting time and pulse width of the signal respectively.

[0235] For different types of frequency modulation signals, the instantaneous phase of the pulse signal can be expressed as:

[0236]

[0237] wherein, represents the initial phase, f1 represents the starting frequency of the pulse signal, f2 represents the cutoff frequency of the pulse signal, λ is the frequency modulation rate of the LFM signal, μ=(f2-f1) / τ represents the periodic slope of the HFM signal, A m represents the modulation index of the sinusoidal wave, fm denotes the modulation frequency of the sinusoidal wave, denotes the initial phase of the sinusoidal wave.

[0238] In the observation signal background noise, w(t) is a Gaussian white noise with zero mean and variance σ 2 2 The size depends on the signal-to-noise ratio SNR: SNR = 10log 10 [A 2 / 2σ 2 ], A is the pulse signal amplitude.

[0239] The above simulation received noisy combined pulse signal x(t) is discretely sampled at a sampling frequency f s , and the noisy combined pulse signal sampling data sequence x(n) is obtained as:

[0240]

[0241] Embodiment 1:

[0242] The simulation signal parameters are set as: signal amplitude A = 1, initial phase pulse width T = 1 s, linear frequency modulation pulse signal starting time t0 = 0.25 s, linear frequency modulation pulse signal pulse width τ = 0.5 s, linear frequency modulation pulse signal starting frequency f1 = 3600 Hz, linear frequency modulation pulse signal terminal frequency f2 = 900 Hz, sampling frequency f s = 10000 Hz, observation data sequence points L = 10000, signal-to-noise ratio SNR = -5 dB.

[0243] In step 2, the initial window length set minimum window length H min is set to 16, H max is set to 10000, S W is set to 4, r a is set to 0.3, l begin is set to 170, and α is set to 0.7.

[0244] In step 3, the time frame of the signal is slid and the adaptive window length set Ω l is updated according to the initial window length set Ω and the time frame parameters.

[0245] In step 4, the kurtosis h(l, m) corresponding to each window length in the current time frame adaptive window length set Ω l is calculated, and the initial adaptive time-frequency analysis window length W adp (l) of the current time frame is obtained according to the maximum kurtosis criterion; in step 5, the continuity criterion of window length change is used to adjust the initial adaptive time-frequency analysis window length W adp ​(l) Make corrections to obtain the optimal time-frequency analysis window length W for the current time frame. opt (l);

[0246] In step 6, based on the optimal time-frequency analysis window length W of the current time frame... opt (l) Obtain the adaptive short-time power spectrum Y of the current time frame. opt (l,k l and instantaneous discretized frequency And based on instantaneous discretization frequency To determine the changing trend and obtain the long set of the contraction window.

[0247] In step 7, the judgment conditions are used to ensure that all L-frame data is processed;

[0248] In step 8, based on the optimal time-frequency analysis window length sequence W opt (l), l=1,2,…,L, reconstruct the adaptive time-frequency analysis results. like Figure 2 As shown, the time-varying window length, selected based on kurtosis maximization and corrected by the continuity criterion, exhibits a certain relationship with the discretized frequency of the signal, reflecting the characteristics of a linear frequency modulated signal. For example... Figure 3 As shown, it can be seen that the adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria obtained a high-quality time-frequency distribution, and the results showed good time-frequency characteristics of linear frequency modulated pulse signals.

[0249] Example 2:

[0250] The simulation signal parameters are set as follows: signal amplitude A = 1, initial phase... The pulse width T = 1.4s, the start time t0 = 0.2s, the pulse width τ = 0.9s, the start frequency f1 = 100Hz, the end frequency f2 = 300Hz, and the sampling frequency f... s =1000Hz, observed data sequence points L=1400, signal-to-noise ratio SNR=-5dB.

[0251] In step 2, the initial window length set is set to the minimum window length H. min The value is 8, H max The value is 1400, S W The value is 2, r a The value is 0.3, l begin The value is 170, and α is 0.7;

[0252] In step 3, the time frame of the sliding signal is determined, and the adaptive window length set Ω is updated according to the initial window length set Ω and the time frame parameters. l ;

[0253] In step 4, the adaptive window length set Ω for the current time frame is calculated. l The kurtosis h(l,m) corresponding to each window length is used to obtain the initial adaptive time-frequency analysis window length W for the current frame based on the kurtosis maximization criterion. adp (l);

[0254] In step 5, based on the continuity criterion of window length variation, the initial adaptive time-frequency analysis window length W for the current time frame is determined. adp (l) Make corrections to obtain the optimal time-frequency analysis window length W for the current time frame. opt (l);

[0255] In step 6, based on the optimal time-frequency analysis window length W of the current time frame... opt (l) Obtain the adaptive short-time power spectrum Y of the current time frame. opt (l,k l and instantaneous discretized frequency And based on instantaneous discretization frequency To determine the changing trend and obtain the long set of the contraction window.

[0256] In step 7, the judgment conditions are used to ensure that all L-frame data is processed;

[0257] In step 8, based on the optimal time-frequency analysis window length sequence W opt (l), l=1,2,…,L, reconstruct the adaptive time-frequency analysis results. like Figure 4 As shown, the time-varying window length, selected based on kurtosis maximization and corrected by the continuity criterion, exhibits a certain relationship with the discretized frequency of the signal, reflecting the characteristics of hyperbolic frequency modulated signals. Figure 5 As shown, the adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria yielded a high-quality time-frequency distribution, and the results exhibited good time-frequency characteristics of hyperbolic frequency-modulated pulse signals.

[0258] It is understood that the present invention has been described through some embodiments, and those skilled in the art will recognize that various changes or equivalent substitutions can be made to these features and embodiments without departing from the spirit and scope of the invention. Furthermore, under the teachings of the present invention, these features and embodiments can be modified to adapt to specific situations and materials without departing from the spirit and scope of the invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the protection scope of the present invention.

Claims

1. A method of adaptive time-frequency analysis based on kurtosis maximization and continuity criterion, characterized in that, The method comprises the following steps: Step 1, obtaining a signal sampling data sequence x(n) to be processed, n = 1, 2, …, N, n is the sampling point index of x(n), N = 2 a , a is an integer greater than 10; Step 2, initializing an initial window length set Omega of adaptive time-frequency analysis and time frame parameters; Step 3, slide the time frame of the signal and update the adaptive window length set Ω according to the initial window length set Ω and the time frame parameter l ; Step 4, calculate the current time frame adaptive window length set Ω l The kurtosis h(l, m) corresponding to each window length is obtained, and the initial adaptive time-frequency analysis window length W of the current time frame is obtained according to the kurtosis maximization criterion adp (l); Step 5. Obtain the optimal time-frequency analysis window length W for the current time frame based on the continuity criterion of window length variation adp (l) performing a modification to obtain the optimal time-frequency analysis window length W for the current time frame opt (l); Step 6, obtaining the best time-frequency analysis window length W of the current time frame according to the current time frame opt (l), obtaining the adaptive short-time power spectrum Y of the current time frame opt (l,k l ) and the instantaneous discrete frequency and obtaining the set of contraction window lengths according to the variation trend of the instantaneous discrete frequency Step 7, judging whether all time frames are processed, if not, returning to step 3, otherwise, going to step 8; Step 8, determining the best time-frequency analysis window length W according to the reconstructed adaptive time-frequency analysis result opt (l), l = 1, 2, …, L, reconstructed adaptive time-frequency analysis result where L is the total number of time frames.

2. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 1, characterized in that, The specific method of step 1 is as follows: Receive real-time acquisition data of N sampling points from the sensor as the data sequence x(n) to be processed, or extract N sampling point data containing the entire pulse signal to be analyzed from the memory as the data sequence x(n) to be processed, N=2 a a is an integer greater than 10.

3. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 1, characterized in that, The specific method of step 2 is as follows: The time-frequency parameters required for adaptive time-frequency analysis of the pulse signal to be analyzed are initialized, specifically including the initialization of the following parameters: Minimum initial window length H for time-frequency analysis of the pulse signal to be analyzed min initialized to: 4 ≤ H min ≤ 16, and is even; The maximum initial window length H of time-frequency analysis of the pulse signal to be analyzed max Initialization: H max = N; Window length search step S for time-frequency analysis of the pulse signal to be analyzed W initialized to: 1 < S W <5, and even, and satisfies (H max - H min ) / S W is an integer; Threshold r of the rate of change of the time-frequency analysis window length when analyzing the pulse signal a initialized to: 0.25 < r a <0.45, and is a real number; minimum initial time frame sequence index l at which the adaptive window length set starts to shrink begin initialized as: N / 8 < l begin < N / 4, and l begin is even; Let the initial window length set Omega of the pulse signal to be analyzed for time-frequency analysis be: Ω = {H min ,H min +S W ,H min +2S W ,...,H max} The sliding window step S is initialized as S = int[H min / 4], and the total time frame number L of the pulse signal adaptive time-frequency analysis to be analyzed is initialized as: L = int[N / S] - 1, wherein int[] represents the downward integer of [] inside. The time frame sequence index l of the adaptive time-frequency analysis of the pulse signal to be analyzed is initialized as: l=0.

4. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 3, characterized in that, The specific method of step 3 is as follows: Step 3.1, sliding the time frame of the signal, that is, setting the time frame sequence index l=l+1; Step 3.2 Set of adaptive window lengths Ω for the 1st time frame l is updated to: The updated adaptive window length set Ω of the 1th time frame based on the above formula l is: Ω l = {W min (l), W min (l) + S W , W min (l) + 2S W ,..., W max (l)} where W max (l) is the maximum window length of the adaptive window set Ω l (l) for the lth time frame, and W min (l) is the minimum window length of the adaptive window set Ω l (l) for the lth time frame.

5. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 3, characterized in that, The specific method of step 4 is as follows: Step 4.1, count the adaptive window length set Ω of the lth time frame l the number of medium window length elements, denoted as M l ; Step 4.2, define m as the 1st frame adaptive window length set Ω l the window length index in Ω, and initialize m = 0; Step 4.3, update the set of adaptive window lengths Ω l the window length index m within the set of adaptive window lengths Ω l the mth window length W(l,m) within the set of adaptive window lengths Ω W(l,m) = W min (l) + (m - 1) - S W ; Step 4.4, calculate the time-varying spectrum y corresponding to the mth window length of the lth time frame m (l, k l ): where k l = 0, 1,..., W max (l) 2 - 1, denotes the discrete frequency index of the lth time frame, j denotes the imaginary unit, i.e. The value range of n satisfies n1(l,m)≤n≤n2(l,m), wherein the starting position n1(l,m) and the ending position n2(l,m) of the sliding window are respectively: n1(l,m)=max{1,int[(l-1)S+1-W(l,m)2]} n2(l,m)=min{N,int[(l-1)S+1+W(l,m)2]} Where, min{·} represents taking the minimum value of two, and max{·} represents taking the maximum value of two; Step 4.5, calculate the short-time power spectrum Y corresponding to the mth window length of the current lth time frame m (l, k l ): Where, | | represents the modulus operation; Step 4.6, normalizing the short-time power spectrum Y m (l,k l ) corresponding to the mth window length of the current ith time frame to obtain the normalized short-time power spectrum wherein, represents the discretized frequency k l the maximum value of this dimension; Step 4.7, obtaining the normalized short-time power spectrum corresponding to the mth window length of the current lth time frame calculating the kurtosis h(l,m) corresponding to the mth window length of the current lth time frame wherein, is the normalized short-time power spectrum corresponding to the lth time frame and mth window length at the discrete frequency k l Mean value in this dimension: Step 4.8, determine if m > M l is true, if true, go to step 4.9, otherwise, go back to step 4.3; Step 4.9: Based on the M of the current l-th time frame l M corresponding to the length of each window l Given a kurtosis h(l,m), obtain the initial adaptive time-frequency analysis window length W for the current l-th frame. adp (l): wherein denotes W(l,m) such that the maximum value within the brackets is taken.

6. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 5, characterized in that, The specific method of step 5 is as follows: Step 5.1, judge whether l = 1 is true, if true, let the best time-frequency analysis window length of the lth time frame equal the initial adaptive time-frequency analysis window length of the lth time frame, i.e. W opt (l) = W adp (l), go to step 6, otherwise, go to step 5.2; Step 5.2, Calculate the optimal time-frequency analysis window length W for the lth time frame opt (l): Where, r(l) represents the window length change rate of the lth time frame:

7. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 6, characterized in that, The specific method of step 6 is as follows: Step 6.1, determine the optimal time-frequency analysis window length W according to the current time frame l opt (l), obtain the adaptive time-varying spectrum y opt (l, k l ) where n is an integer satisfying n'1(l,m)≤n≤n'2(l,m), n'1(l,m) and n'2(l,m) represent the optimal time-frequency analysis window length W opt (l) the start position and the end position of the sliding window are calculated as n'1(l,m) = max{1, int[(l - 1)S + 1 - W opt (l)2]} n'2(l,m) = min{N, int[(l - 1)S + 1 + W opt (l)2]} n'2(l,m) = min{N, int[(l - 1)S + 1 + W <000007 Step 6.

2. Compute the adaptive time-varying spectrum y opt (l,k l ) of the lth time frame from the adaptive short-time power spectrum Y opt (l,k l ) of the lth time frame. Step 6.

3. Estimate the instantaneous discrete frequency f of the l-th time frame of the signal to be analyzed from the adaptive short-time power spectrum Y opt (l,k l ) of the l-th time frame of the signal to be analyzed wherein denotes k such that the maximum value within the brackets is obtained l value; Step 6.4: Determine if l ≥ l begin Whether it is true, that is, whether l has been accumulated. begin If the information of each time frame is true, proceed to step 6.5; otherwise, proceed to step 6.

9. Step 6.5, the direction trend factor a is initialized as: 0.6< a <0.8, and a is a real number; Step 6.6 Calculate the instantaneous discretized frequency using the Cox-Stuart trend test algorithm The trend of the change is calculated as follows in sub-steps: Step 6.6.1, from the 1st time frame, take the instantaneous discretized frequency of the l+1 time frames, and include the 1st time frame, to make a data set begin ​ Take and comprise a set of data wherein, i represents the trend test data index, is the left value of the trend test data pair, is the right value of the trend test data pair i=l-l begin +1,l-l begin +2,…,l-c,i+c=l-c+1,l-c+2,…,l, wherein, c=int[l begin 2] represents the number of trend test data pairs; Step 6.6.2, difference the two numbers in each data pair in the time direction, record the trend value And introduce the number of positive trends S + Record D i The number of positive, the number of negative trends S - Record D i The number of negative Step 6.6.

3. Calculate the number of positive trends S + and the number of negative trends S - Calculate the test result T: Step 6.7, if the test result T is 0, go to step 7, otherwise go to step 6.8; Step 6.

8. Check the trend of the 1st time frame T and the first l begin - the optimal time-frequency analysis window length W of the 1st time frame opt (l), obtain the set of contraction window lengths of the 1st time frame The expression is: where W' (l) is the set of window lengths of the lth time frame min (l) and W' (l) represent the set of window lengths of the lth time frame max (l) represents the set of window lengths of the lth time frame The minimum window length and the maximum window length of W' (l) are calculated as follows: W mid (l) denotes the first l begin - the median value of the optimal time-frequency analysis window length, W mid (l) is calculated as: W mid (l) = median{W opt (l-l begin +1),W opt (l-l begin +2),…,W opt (l-1)} where median{•} denotes taking the median of the values within •, taking the median if the length of • is odd, and taking the average of the middle two values if the length of • is even, since l begin is initialized to even in step 2, l begin -1 is odd, W mid (l) is the median of the sequence; Step 6.

9. Let the set of adaptive window lengths of the ith time frame equal the set of contracted window lengths of the ith time frame, i.e.

8. The adaptive time-frequency analysis method based on the maximization of kurtosis and continuity criterion according to claim 7, characterized in that, The specific method of step 7 is as follows: Judging whether l < L, if yes, returning to step 3.1, otherwise, going to step 8.

9. The method of adaptive time-frequency analysis based on kurtosis maximization and continuity criterion according to claim 8, characterized in that, The specific method of step 8 is as follows: Step 8.1, Normalizing the spectrum length L norm initialized to: wherein denotes the maximum value in the dimension of the time frame sequence index l; And reinitialize the time frame sequence index l=0; Step 8.2, slide the signal's time frame, i.e. let the time frame sequence index l = l + 1, and determine the optimal time-frequency analysis window length W according to the l-th time frame opt (l), compute the normalized length-adaptive time-varying spectrum y norm (l,k) : where n is an integer satisfying n'1(l,m)≤n≤n'2(l,m), n'1(l,m) and n'2(l,m) are defined as in step 6.1; k represents the length of the normalized spectrum L norm the resulting discretized frequency index; Step 8.

3. Normalized spectrum length adaptive short-time power spectrum Y of the 1st time frame norm (l,k) is calculated as: Step 8.4, judging whether l=L is true, if yes, go to step 8.5; otherwise, return to step 8.2; Step 8.5, amplitude-normalize the L frame of normalized spectrum length adaptive short-time power spectrum Y norm (l,k) in time sequence, to obtain the final time-frequency analysis result norm (l,k) in time sequence, to obtain the final time-frequency analysis result wherein denotes the maximum value in both dimensions l, k.

Images