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

Abstract

The invention discloses a self-adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion. The method comprises the following steps: 1, acquiring a pulse signal sampling data sequence to be processed; 2, initializing an initial window length set and time frame parameters of adaptive time-frequency analysis; 3, updating the adaptive window length set according to the initial window length set and the time frame parameters; 4, calculating the kurtosis 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 of the current time frame; 5, correcting the initial adaptive time-frequency analysis window length according to a window length change continuity criterion to obtain the optimal time-frequency analysis window length of the current time frame; 6, obtaining the self-adaptive short-time power spectrum and the instantaneous discretization frequency of the current time frame, and obtaining a contraction window length set according to the change trend of the instantaneous discretization frequency; 7, judging whether all time frames are processed or not, if not, returning to the step 3, and otherwise, turning to the step 8; and 8, reconstructing an adaptive time-frequency analysis result according to the optimal time-frequency analysis window length sequence.

Description

Technical Field

[0001] The present invention belongs to the technical field of signal processing, and in particular relates to an adaptive time-frequency analysis method based on kurtosis maximization and continuity characteristics. Background Art

[0002] Since the spectrum of non-stationary signals is mostly time-varying, the classic Fourier transform cannot fully reflect the changes of the frequency components of the signal over time. Therefore, in order to reveal the global information of the signal in the time and frequency dimensions, it is required to perform joint processing in the time and frequency domains to obtain more information. Therefore, time-frequency analysis is proposed to solve the analysis problem of non-stationary signals.

[0003] Scholars in the field of signal processing usually divide time-frequency analysis into linear time-frequency analysis and quadratic time-frequency analysis. The representative classical algorithm of quadratic time-frequency analysis is Wigner-Ville distribution (WVD), which has high-resolution time-frequency analysis performance. However, for multi-component signals or nonlinear frequency modulation signals, this algorithm and its subsequent improved algorithms are troubled by cross terms, and the amount of calculation is significantly higher than that of linear time-frequency analysis methods, so it is difficult to deploy and apply in real-time processing systems. The representative classical algorithm of linear time-frequency analysis is short-time Fourier transform (STFT), which is easy to calculate, has clear physical meaning, and has no cross terms, making it widely used in radar or sonar signal processing. 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. Therefore, the short-time analysis window length selection method used to balance the time and frequency domain resolution is the core of the adaptive STFT method. Usually, the selection of the short-time analysis window length of STFT depends on empirical debugging.

[0004] In order to design time-frequency analysis methods more flexibly, an adaptive linear time-frequency analysis method was proposed to solve the problem of selecting the short-time analysis window length. The existing adaptive time-frequency analysis method has overcome the disadvantage that the window length of STFT does not change with the signal characteristics to a certain extent, and achieved an effective compromise between time and frequency resolution, but the redundancy and computational complexity have been significantly increased. Summary of the invention

[0005] Purpose of the invention: The purpose of the present invention is to provide an adaptive time-frequency analysis method based on kurtosis maximization and continuity characteristics, so as to solve the technical problems that conventional time-frequency analysis methods for pulse signals are difficult to achieve time-frequency domain resolution compromise, the adaptive time-frequency analysis method has redundancy, large amount of calculation, and is difficult to be deployed and used in real-time signal processing systems.

[0006] Technical solution, in order to solve the above technical problems, the present invention proposes an adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion, the method comprising the following steps:

[0007] Step 1: Get the 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 Ω and time frame parameters of the adaptive time-frequency analysis;

[0009] 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 ;

[0010] Step 4: Calculate the adaptive window length set Ω for the current time frame l The kurtosis h(l,m) corresponding to each window length in the current time frame is obtained according to the kurtosis maximization criterion. adp (l);

[0011] Step 5: According to the continuity criterion of window length change, the initial adaptive time-frequency analysis window length W of the current time frame is adjusted. adp (l) Make corrections to obtain the optimal time-frequency analysis window length W for the current time frame opt (l);

[0012] 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 the instantaneous discretization frequency And according to the instantaneous discretization frequency The change trend of the shrinking window length set is obtained

[0013] Step 7: Determine whether all time frames have been processed. If not, return to step 3; otherwise, go to step 8.

[0014] Step 8: Based on the optimal time-frequency analysis window length W opt (l), l = 1, 2, ..., L, reconstruct the adaptive time-frequency analysis results Where L is the total number of time frames.

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

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

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

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

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

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

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

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

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

[0024] Let the initial window length set Ω for 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] Initialize the sliding window step S to S = int[H min / 4], then initialize the total number of time frames L for adaptive time-frequency analysis of the pulse signal to be analyzed to: L = int[N / S] - 1, where int[] represents rounding down the value inside [];

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

[0028] Furthermore, the specific method for step 3 is as follows:

[0029] Step 3.1, sliding the time frame of the signal, that is, setting the time frame sequence index l=l+1;

[0030] Step 3.2: The adaptive window length set Ω of the lth time frame l Updated to:

[0031]

[0032] Based on the updated adaptive window length set Ω of the lth time frame l for:

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

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

[0035] Furthermore, the specific method of step 4 is as follows:

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

[0037] Step 4.2: Define m as the adaptive window length set Ω of the lth time frame l The window length index within is initialized to m=0;

[0038] Step 4.3: Update the adaptive window length set Ω of the lth time frame l The window length index in the lth time frame is set to m=m+1, and the adaptive window length set Ω is taken. l The length of the mth window W(l,m) is:

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

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

[0041]

[0042] Among them, kl =0,1,...,W max (l) / 2-1, represents the discretized frequency index of the lth time frame, j represents the imaginary unit, that is,

[0043] The value range of n is an integer satisfying n1(l,m)≤n≤n2(l,m), where the starting position n1(l,m) and the 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] Among them, min{·} means taking the minimum value between the two, and max{·} means taking the maximum value between the two;

[0047] 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 ):

[0048]

[0049] Among them, | | represents the modulus operation;

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

[0051]

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

[0053] Step 4.7: Normalized short-time power spectrum corresponding to the mth window length of the current lth time frame Calculate the kurtosis h(l,m) corresponding to the mth window length of the current lth time frame:

[0054]

[0055] in, is the normalized short-time power spectrum corresponding to the mth window length of the lth time frame at the discretized frequency k l The mean value on this dimension is:

[0056]

[0057] Step 4.8: Determine whether m≥M l Is it true? If so, go to step 4.9; otherwise, return to step 4.3;

[0058] Step 4.9: Based on the M of the current lth time frame l The window length corresponding to M l The kurtosis h(l,m) is obtained to obtain the initial adaptive time-frequency analysis window length W of the current lth time frame adp (l):

[0059]

[0060] in, It represents W(l,m) that takes the maximum value in the brackets.

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

[0062] Step 5.1, determine whether l = 1. If yes, set the optimal time-frequency analysis window length of the lth time frame equal to the initial adaptive time-frequency analysis window length of the lth time frame, that is, W opt (l) = W adp (l), go to step 6, otherwise, go to step 5.2;

[0063] Step 5.2: Calculate the optimal time-frequency analysis window length W for the first time frame opt (l):

[0064]

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

[0066]

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

[0068] Step 6.1: Based on the optimal time-frequency analysis window length W of the current lth time frame opt (l), obtain the adaptive time-varying spectrum y of the current lth time frame opt (l,k l )

[0069]

[0070] Among them, the value range of n satisfies the integer 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 according to the lth time frame respectively. opt(l) The starting position and ending position of the sliding window obtained are calculated as:

[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 first time frame opt (l,k l ), calculate the adaptive short-time power spectrum Y of the lth time frame opt (l,k l ):

[0074]

[0075] Step 6.3: Based on the adaptive short-time power spectrum Y of the first time frame opt (l,k l ), estimate the instantaneous discretization frequency of the lth time frame of the signal to be analyzed

[0076]

[0077] in, Indicates k that takes the maximum value in the brackets l value;

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

[0079] Step 6.5, the directional trend factor α is initialized to: 0.6<α<0.8, α is a real number;

[0080] Step 6.6: Calculate the instantaneous discretization frequency using the Cox-Stuart trend test algorithm The sub-steps are as follows:

[0081] Step 6.6.1, starting from the lth time frame, take l forward begin The instantaneous discretization frequency of time frames, including the lth time frame, constitutes the data set

[0082]

[0083] Pick and Make up a set of data Among them, i represents the trend test data index, is the left value of the trend test data pair, The right value of the trend test data pair i = ll begin +1,ll begin +2,…,lc,i+c=l-c+1,l-c+2,…,l,where c=int[l begin / 2], indicating the number of trend test data pairs;

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

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

[0086]

[0087] Step 6.7, if the test result T is 0, go to step 7, otherwise go to step 6.8;

[0088] Step 6.8: Check T and the first l of the l time frame according to the trend of the l time frame. begin -The optimal time-frequency analysis window length W for 1 time frame opt (l), get the contraction window length set of the lth time frame The expression is:

[0089]

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

[0091]

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

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

[0094] Among them, median{·} represents taking the median value of the values in the sequence within {·}. If the sequence length is odd, the median is taken; if the sequence length is even, the average of the two middle numbers is taken. Since l begin is initialized as an even number in Step 2, l begin - 1 is an odd number, and W mid (l) is the median of the sequence;

[0095] Step 6.9: Let the adaptive window length set of the l-th time frame be equal to the contracted window length set of the l-th time frame, that is

[0096] Furthermore, the specific method of Step 7 is as follows:

[0097] Judge whether l < L. If it holds, return to Step 3.1; otherwise, go to Step 8.

[0098] Furthermore, the specific method of Step 8 is as follows:

[0099] Step 8.1: Normalize the spectral length L norm Initialize it as:

[0100]

[0101] Among them, represents the maximum value in the dimension of the time frame sequence index l;

[0102] And re-initialize the time frame sequence index l = 0;

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

[0104]

[0105] Among them, the value range of n satisfies the integer of n'1(l, m) ≤ n ≤ n'2(l, m), and the definitions of n'1(l, m) and n'2(l, m) are as described in Step 6.1; k represents the discretized frequency index obtained according to the normalized spectral length L norm ;

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

[0107]

[0108] Step 8.4: Determine whether l=L is established. If so, go to step 8.5; otherwise, return to step 8.2;

[0109] Step 8.5: Adapt the normalized spectrum length of the L frame to the short-time power spectrum Y norm (l,k) are arranged in time order, and the normalized spectrum length of L frames is adaptive to the short-time power spectrum Y norm (l,k) performs amplitude normalization to obtain the final time-frequency analysis results

[0110]

[0111] in, Indicates the maximum value in the two dimensions l and k.

[0112] Beneficial effects: Compared with the prior art, the technical solution of the present invention has the following beneficial technical effects:

[0113] 1. The time-frequency analysis method of the present invention utilizes the frequency domain narrowband characteristics of the local signal when obtaining the short-time Fourier transform of the signal, obtains the statistical value describing the time-frequency energy concentration of the signal, namely, the spectrum kurtosis, and adaptively selects the short-time analysis window length to perform short-time Fourier transform on the local signal. In steps 2 and 3, the high-energy concentration spectrum of the signal is adaptively obtained through the unified spectrum kurtosis of each window length in the window length set;

[0114] 2. When obtaining the adaptive short-time analysis window length and short-time Fourier transform of the signal, the time-frequency analysis method of the present invention utilizes the intrinsic connection between the signal analysis window length and the instantaneous discretization frequency of the signal to be analyzed in step 4. By calculating the rate of change of the signal adaptive window length, the continuity criterion of the signal is used to correct the adaptive short-time analysis window length sequence. The rate of change of the window length compared with the previous time frame is calculated by the obtained signal adaptive window length, and then compared with the preset rate of change threshold, and the possible abnormal analysis window length is corrected, so that a relatively ideal short-time Fourier transform result of the signal can be calculated at a lower signal-to-noise ratio;

[0115] 3. When obtaining the adaptive short-time analysis window length and short-time Fourier transform results of several time frames of the signal, the time-frequency analysis method of the present invention utilizes the characteristics of common pulse signals and judges the change trend of the instantaneous discretization 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 invention. In step 5, the normalized short-time power spectrum Y opt (l,k l ) estimates the instantaneous discretization frequency of the signal. In step 6, the Cox-Stuart trend test algorithm is used to test the increase and decrease trend of the signal to guide the shrinking range and shrinking direction of the adaptive time-frequency analysis window length set. In step 7, the ideal time-frequency analysis result Y is finally obtained. norm (l,k). By adaptively shrinking the window length set, the redundancy and large computational complexity of adaptive time-frequency analysis are suppressed, and the robustness of the time-frequency analysis results under low signal-to-noise ratio is improved. The present invention makes full use of the short-term narrowband characteristics of common pulse signals and can obtain adaptive time-frequency analysis results with high time-frequency focus in a relatively robust manner. BRIEF DESCRIPTION OF THE DRAWINGS

[0116] Figure 1 It is a flow chart of the adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion of the present invention;

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

[0118] Figure 3 This is a diagram showing the normalized adaptive time-frequency analysis results of a simulated pulse signal in Example 1;

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

[0120] Figure 5 This is a diagram showing the normalized adaptive time-frequency analysis results of a simulated pulse signal in Example 2; DETAILED DESCRIPTION

[0121] In order to better understand the purpose, structure and function of the present invention, the following is a further detailed description of an adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion of the present invention in conjunction with the accompanying drawings.

[0122] like Figure 1 As shown, the present invention proposes an adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria, the method comprising the following steps:

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

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

[0125] 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 ;

[0126] Step 4: Calculate the adaptive window length set Ω for the current time frame l The kurtosis h(l,m) corresponding to each window length in the current time frame is obtained according to the kurtosis maximization criterion. adp (l);

[0127] Step 5: According to the continuity criterion of window length change, the initial adaptive time-frequency analysis window length W of the current time frame is adjusted. adp (l) Make corrections to obtain the optimal time-frequency analysis window length W for the current time frame opt (l);

[0128] 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 the instantaneous discretization frequency And according to the instantaneous discretization frequency The change trend of the shrinking window length set is obtained

[0129] Step 7: Determine whether all time frames have been processed. If not, return to step 3; otherwise, go to step 8.

[0130] Step 8: Based on the optimal time-frequency analysis window length W opt (l), l = 1, 2, ..., L, reconstruct the adaptive time-frequency analysis results Where L is the total number of time frames.

[0131] Furthermore, the specific method of step 1 is as follows:

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

[0133] Furthermore, the specific method of step 2 is as follows:

[0134] Initialize the 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] The minimum initial window length H for time-frequency analysis of the pulse signal to be analyzed min Initialize to: 4 ≤ H min ≤ 16, and it is an even number;

[0136] The maximum initial window length H for time-frequency analysis of the pulse signal to be analyzed max Initialize to: H max = N;

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

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

[0139] The minimum initial time frame sequence index l at which the adaptive window length set starts to shrink begin Initialize to: N / 8 < l begin < N / 4, and l begin is an even number;

[0140] Let the initial window length set Ω for 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 to S = int[H min / 4], then initialize the total number of time frames L for adaptive time-frequency analysis of the pulse signal to be analyzed to: L = int[N / S] - 1, where int[] represents rounding down the value inside [];

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

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

[0145] Step 3.1, Slide the time frame of the signal, that is, let the time frame sequence index l = l + 1;

[0146] Step 3.2: The adaptive window length set Ω of the lth time frame l Updated to:

[0147]

[0148] Based on the updated adaptive window length set Ω of the lth time frame l for:

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

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

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

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

[0153] Step 4.2: Define m as the adaptive window length set Ω of the lth time frame l The window length index within is initialized to m=0;

[0154] Step 4.3: Update the adaptive window length set Ω of the lth time frame l The window length index in the lth time frame is set to m=m+1, and the adaptive window length set Ω is taken. l The length of the mth window W(l,m) is:

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

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

[0157]

[0158] Among them, k l =0,1,...,Wmax (l) / 2-1, represents the discretized frequency index of the lth time frame, j represents the imaginary unit, that is,

[0159] The value range of n is an integer satisfying n1(l,m)≤n≤n2(l,m), where 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] Among them, min{·} means taking the minimum value between the two, and max{·} means taking the maximum value between the two;

[0163] 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 ):

[0164]

[0165] Among them, | | represents the modulus operation;

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

[0167]

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

[0169] Step 4.7: Normalized short-time power spectrum corresponding to the mth window length of the current lth time frame Calculate the kurtosis h(l,m) corresponding to the mth window length of the current lth time frame:

[0170]

[0171] in, is the normalized short-time power spectrum corresponding to the mth window length of the lth time frame at the discretized frequency k l The mean value on this dimension is:

[0172]

[0173] Step 4.8: Determine whether m≥M l Is it true? If so, go to step 4.9; otherwise, return to step 4.3;

[0174] Step 4.9: Based on the M of the current lth time frame l The window length corresponding to M l The kurtosis h(l,m) is obtained to obtain the initial adaptive time-frequency analysis window length W of the current lth time frame adp (l):

[0175]

[0176] in, It represents W(l,m) that takes the maximum value in the brackets.

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

[0178] Step 5.1, determine whether l = 1. If yes, set the optimal time-frequency analysis window length of the lth time frame equal to the initial adaptive time-frequency analysis window length of the lth time frame, that is, W opt (l) = W adp (l), go to step 6, otherwise, go to step 5.2;

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

[0180]

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

[0182]

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

[0184] Step 6.1: Based on the optimal time-frequency analysis window length W of the current lth time frame opt (l), obtain the adaptive time-varying spectrum y of the current lth time frame opt (l,k l )

[0185]

[0186] Among them, the value range of n satisfies the integer 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 according to the lth time frame respectively. opt(l) The starting position and ending position of the sliding window obtained are calculated as:

[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 first time frame opt (l,k l ), calculate the adaptive short-time power spectrum Y of the lth time frame opt (l,k l ):

[0190]

[0191] Step 6.3: Based on the adaptive short-time power spectrum Y of the first time frame opt (l,k l ), estimate the instantaneous discretization frequency of the lth time frame of the signal to be analyzed

[0192]

[0193] in, Indicates k that takes the maximum value in the brackets l value;

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

[0195] Step 6.5, the directional trend factor α is initialized to: 0.6<α<0.8, α is a real number;

[0196] Step 6.6: Calculate the instantaneous discretization frequency using the Cox-Stuart trend test algorithm The sub-steps are as follows:

[0197] Step 6.6.1, starting from the lth time frame, take l forward begin The instantaneous discretization frequency of time frames, including the lth time frame, constitutes the data set

[0198]

[0199] Pick and Make up a set of data Among them, i represents the trend test data index, is the left value of the trend test data pair, The right value of the trend test data pair i = ll begin +1,ll begin +2,…,lc,i+c=l-c+1,l-c+2,…,l,where c=int[l begin / 2], indicating the number of trend test data pairs;

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

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

[0202]

[0203] Step 6.7, if the test result T is 0, go to step 7, otherwise go to step 6.8;

[0204] Step 6.8: Check T and the first l of the l time frame according to the trend of the l time frame. begin -The optimal time-frequency analysis window length W for 1 time frame opt (l), get the contraction window length set of the lth time frame The expression is:

[0205]

[0206] Among them, W' min (l) and W' max (l) represents the set of contraction window lengths of the lth time frame The minimum and maximum window lengths are calculated as follows:

[0207]

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

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

[0210] where median{·} represents taking the median value of the values in the sequence within {·}. If the sequence length is odd, the median is taken; if the sequence length is even, the average of the two middle numbers is taken. Since l begin is initialized as an even number in step 2, l begin - 1 is odd, and W mid (l) is the median of the sequence;

[0211] Step 6.9. Let the adaptive window length set at the l-th time frame be equal to the shrunk window length set at the l-th time frame, i.e.,

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

[0213] Judge whether l < L. If it holds, return to step 3.1; otherwise, go to step 8.

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

[0215] Step 8.1. Normalize the spectral length L norm Initialize it as:

[0216]

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

[0218] And re-initialize 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 calculate the normalized length adaptive time-frequency spectrum y opt (l, k) of the l-th time frame according to the optimal time-frequency analysis window length W norm (l):

[0220]

[0221] where the value range of n satisfies n'1(l, m) ≤ n ≤ n'2(l, m) as an integer, and the definitions of n'1(l, m) and n'2(l, m) are as described in step 6.1; k represents the discretized frequency index obtained according to the normalized spectral length L norm ;

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

[0223]

[0224] Step 8.4: Determine whether l=L is established. If so, go to step 8.5; otherwise, return to step 8.2;

[0225] Step 8.5: Adapt the normalized spectrum length of the L frame to the short-time power spectrum Y norm (l,k) are arranged in time order, and the normalized spectrum length of L frames is adaptive to the short-time power spectrum Y norm (l,k) performs amplitude normalization to obtain the final time-frequency analysis results

[0226]

[0227] in, Indicates the maximum value in the two dimensions l and k.

[0228] In the embodiment of the present invention, the simulation receiving pulse signal model is:

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

[0230] Among them, w(t) is the observed signal background noise, which is Gaussian white noise; s(t) is the pulse signal, expressed as:

[0231]

[0232] Where 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] Where t0 and τ are the start time and pulse width of the signal respectively.

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

[0236]

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

[0238] In the observed signal background noise, w(t) has a mean of zero and a variance of σ 2 Gaussian white noise with variance σ 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] With sampling frequency f s The noisy combined pulse signal x(t) received in the above simulation is discretely sampled to obtain a noisy combined pulse signal sampling data sequence x(n):

[0240]

[0241] Embodiment 1:

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

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

[0244] In step 3, the time frame of the sliding signal is updated based on the initial window length set Ω and the time frame parameters. l ;

[0245] In step 4, the adaptive window length set Ω of the current time frame is calculated l The kurtosis h(l,m) corresponding to each window length in the current time frame is obtained according to the kurtosis maximization criterion. adp (l); In step 5, according to the continuity criterion of window length change, the initial adaptive time-frequency analysis window length W of the current time frame is 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, according to 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 the instantaneous discretization frequency And according to the instantaneous discretization frequency The change trend of the shrinking window length set is obtained

[0247] In step 7, ensure that all L frames of data are processed by judging the conditions;

[0248] In step 8, according to 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 in , the time-varying window length selected based on kurtosis maximization and corrected by the continuity criterion shows a certain relationship with the discretization frequency of the signal, reflecting the characteristics of the linear frequency modulation signal. Figure 3 As shown, it can be seen that the adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria obtains high-quality time-frequency distribution, and the result presents good time-frequency characteristics of linear frequency modulation pulse signals.

[0249] Embodiment 2:

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

[0251] In step 2, the minimum window length H of the initial window length set is set 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 the value of α is 0.7;

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

[0253] In step 4, the adaptive window length set Ω of the current time frame is calculated l The kurtosis h(l,m) corresponding to each window length in the current time frame is obtained according to the kurtosis maximization criterion. adp (l);

[0254] In step 5, according to the continuity criterion of window length change, the initial adaptive time-frequency analysis window length W of the current time frame is adjusted. 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, according to 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 the instantaneous discretization frequency And according to the instantaneous discretization frequency The change trend of the shrinking window length set is obtained

[0256] In step 7, ensure that all L frames of data are processed by judging the conditions;

[0257] In step 8, according to 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 in the figure, the time-varying window length selected based on kurtosis maximization and corrected by the continuity criterion shows a certain change relationship with the discretization frequency of the signal, reflecting the characteristics of the hyperbolic frequency modulation signal. Figure 5 As shown, it can be seen that the adaptive time-frequency analysis method based on kurtosis maximization and continuity criteria obtains high-quality time-frequency distribution, and the result shows good time-frequency characteristics of hyperbolic frequency modulated pulse signals.

[0258] It is to be understood that the present invention is described by some embodiments, and it is known to those skilled in the art that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teachings of the present invention, these features and embodiments may be modified to adapt to specific circumstances and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the scope of protection of the present invention.

Claims

1. An adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion, characterized in that: The method includes the following steps: Step 1: Get the 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, initialize the initial window length set Ω and time frame parameters of the adaptive time-frequency analysis; 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 parameters. l ; Step 4: Calculate the adaptive window length set Ω for the current time frame l The kurtosis h(l,m) corresponding to each window length in the current time frame is obtained according to the kurtosis maximization criterion. adp (l); Step 5: According to the continuity criterion of window length change, the initial adaptive time-frequency analysis window length W of the current time frame is adjusted. adp (l) Make corrections to obtain the optimal time-frequency analysis window length W for the current time frame opt (l); 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 the instantaneous discretization frequency And according to the instantaneous discretization frequency The change trend of the shrinking window length set is obtained Step 7, determine whether all time frames have been processed. If not, return to Step 3; otherwise, go to Step 8; Step 8: Based on the optimal time-frequency analysis window length W opt (l), l = 1, 2, ..., L, reconstruct the adaptive time-frequency analysis results Where L is the total number of time frames.

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

3. The adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion according to claim 1, characterized in that: The specific method of Step 2 is as follows: Initialize the time-frequency parameters required for the adaptive time-frequency analysis of the pulse signal to be analyzed, specifically including the initialization of the following parameters: The minimum initial window length H for time-frequency analysis of the pulse signal to be analyzed min Initialization: 4≤H min ≤16 and an even number; The maximum initial window length H for time-frequency analysis of the pulse signal to be analyzed max Initialized to: H max =N; The window length search step S for the time-frequency analysis of the pulse signal to be analyzed W Initialized to: 1 W <5, and is an even number, and satisfies (H max -H min ) / S W is an integer;​ The threshold of the change rate of the time-frequency analysis window length of the pulse signal to be analyzed a Initialized to: 0.25 <r a <0.45, and is a real number; The 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 Ω for the time-frequency analysis of the pulse signal to be analyzed be: Ω={H min ,H min +S W ,H min +2S W ,...,H max }; The sliding window step S is initialized to S = int[H min / 4], then the total number of time frames L for adaptive time-frequency analysis of the pulse signal to be analyzed is initialized to: L = int[N / S]-1, where int[] means rounding down within []. The time frame sequence index l for the adaptive time-frequency analysis of the pulse signal to be analyzed is initialized to: l = 0.

4. The adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion according to claim 3, characterized in that: The specific method of Step 3 is as follows: Step 3.1, slide the time frame of the signal, that is, let the time frame sequence index l = l + 1; Step 3.2: The adaptive window length set Ω of the lth time frame l Updated to: Based on the updated adaptive window length set Ω of the lth time frame l for: Ω l ={W min (l),W min (l)+S W ,W min (l)+2S W ,...,W max (l)} Among them, W max (l) is the adaptive window length set Ω of the lth time frame l The maximum window length, W min (l) is the adaptive window length set Ω of the lth time frame l The minimum window length.

5. The adaptive time-frequency analysis method based on kurtosis maximization 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 elements in the middle window length, denoted as M l ; Step 4.2: Define m as the adaptive window length set Ω of the lth time frame l The window length index within is initialized to m=0; Step 4.3: Update the adaptive window length set Ω of the lth time frame l The window length index in the lth time frame is set to m=m+1, and the adaptive window length set Ω is taken. l The length of the mth window W(l,m) is: 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 ): Among them, k l =0,1,...,W max (l)2-1, represents the discretized frequency index of the lth time frame, j represents the imaginary unit, that is, The value range of n satisfies n1(l,m) ≤ n ≤ n2(l,m) for integers, where 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 the two, and max{·} represents taking the maximum value of the 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 modulo operation; Step 4.6: The short-time power spectrum Y corresponding to the m-th window length of the current l-th time frame m (l,k l ) is normalized to obtain the normalized short-time power spectrum corresponding to the m-th window length of the l-th time frame in, Indicates that at the discretized frequency k l The maximum value of this dimension; Step 4.7: Normalized short-time power spectrum corresponding to the mth window length of the current lth time frame Calculate the kurtosis h(l,m) corresponding to the mth window length of the current lth time frame: in, is the normalized short-time power spectrum corresponding to the mth window length of the lth time frame at the discretized frequency k l The mean value on this dimension is: Step 4.8: Determine whether m≥M l Is it true? If so, go to step 4.9; otherwise, return to step 4.3; Step 4.9: Based on the M of the current lth time frame l The window length corresponding to M l The kurtosis h(l,m) is obtained to obtain the initial adaptive time-frequency analysis window length W of the current lth time frame adp (l): in, It represents W(l,m) that takes the maximum value in the brackets.

6. The adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion according to claim 5, characterized in that: The specific method of Step 5 is as follows: Step 5.1, determine whether l = 1. If yes, set the optimal time-frequency analysis window length of the lth time frame equal to the initial adaptive time-frequency analysis window length of the lth time frame, that is, 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 first time frame opt (l): where, r(l) represents the window length change rate of the l-th time frame:

7. The adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion according to claim 6, characterized in that: The specific method of Step 6 is as follows: Step 6.1: Based on the optimal time-frequency analysis window length W of the current lth time frame opt (l), obtain the adaptive time-varying spectrum y of the current lth time frame opt (l,k l ) Among them, the value range of n satisfies the integer 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 according to the lth time frame respectively. opt (l) The starting position and ending position of the sliding window obtained 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]} Step 6.2: Based on the adaptive time-varying spectrum y of the first time frame opt (l,k l ), calculate the adaptive short-time power spectrum Y of the lth time frame opt (l,k l ): Step 6.3: Based on the adaptive short-time power spectrum Y of the first time frame opt (l,k l ), estimate the instantaneous discretization frequency of the lth time frame of the signal to be analyzed in, Indicates k that takes the maximum value in the brackets l value; Step 6.4: Determine if l≥l begin Is it true, that is, whether l has been accumulated begin If the information of the time frame is established, go to step 6.5, otherwise, go to step 6.9; Step 6.5, initialize the direction trend factor α to: 0.6 < α < 0.8, where α is a real number; Step 6.6: Calculate the instantaneous discretization frequency using the Cox-Stuart trend test algorithm The sub-steps are as follows: Step 6.6.1, starting from the lth time frame, take l forward begin The instantaneous discretization frequency of the time frames, including the lth time frame, constitutes the data set Pick and Make up a set of data Among them, i represents the trend test data index, is the left value of the trend test data pair, The right value of the trend test data pair i = ll begin +1,ll begin +2,…,lc,i+c=l-c+1,l-c+2,…,l,where c=int[l begin 2], indicating the number of trend test data pairs; Step 6.6.2: Subtract the two numbers in each data pair in the time series direction and record the trend value And introduce the number of positive trends S + Record D i is the number of positive trends, and the number of negative trends is S - Record D i is a negative quantity; Step 6.6.3: Based on 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 T and the first l of the l time frame according to the trend of the l time frame. begin -The optimal time-frequency analysis window length W for 1 time frame opt (l), get the contraction window length set of the lth time frame The expression is: Among them, W′ min (l) and W′ max (l) represents the set of contraction window lengths of the lth time frame The minimum and maximum window lengths are calculated as follows: W mid (l) represents the first l time frame participating in the Cox-Stuart trend test begin -The median value of the optimal time-frequency analysis window length for 1 time frame, 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{·} means taking the median of the values ​​in the sequence in {·}. If the length of the sequence is an odd number, the median is taken; if the length of the sequence is an even number, the mean of the two middle numbers is taken. begin Initialized to an even number in step 2, l begin -1 is an odd number, W mid (l) is the median of the sequence; Step 6.9: Let the adaptive window length set of the lth time frame be equal to the contracted window length set of the lth time frame, that is, 8. The adaptive time-frequency analysis method based on kurtosis maximization and continuity criterion according to claim 7, characterized in that: The specific method of Step 7 is as follows: Determine whether l < L. If it holds, return to Step 3.1; otherwise, go to Step 8.

9. The adaptive time-frequency analysis method 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: Normalize the spectrum length L norm Initialized to: in, Indicates the maximum value of the dimension of time frame sequence index l; And re-initialize the time frame sequence index l = 0; Step 8.2: Slide the time frame of the signal, that is, let the time frame sequence index l=l+1, and according to the optimal time-frequency analysis window length W of the lth time frame opt (l), calculate the normalized length adaptive time-varying spectrum y of the lth time frame norm (l,k): Wherein, the value range of n satisfies the integer n′1(l,m)≤n≤n′2(l,m), and the definitions of n′1(l,m) and n′2(l,m) are as described in step 6.1; k represents the length of the normalized spectrum L norm The discretized frequency index obtained; Step 8.3: Normalized spectrum length adaptive short-time power spectrum Y of the first time frame norm (l,k) is calculated as: Step 8.4, determine whether l = L holds. If it holds, go to Step 8.5; otherwise, return to Step 8.2; Step 8.5: Adapt the normalized spectrum length of the L frame to the short-time power spectrum Y norm (l,k) are arranged in time order, and the normalized spectrum length of L frames is adaptive to the short-time power spectrum Y norm (l,k) performs amplitude normalization to obtain the final time-frequency analysis results in, Indicates the maximum value in two dimensions: l and k.

Images