Harmonic family detection method based on multi-window joint detection and harmonic characteristic evaluation

Abstract

The invention discloses a harmonic family detection method based on multi-window joint detection and harmonic characteristic evaluation. The method comprises the following steps: 1, obtaining a data sequence; 2, initializing parameters; 3, obtaining a normalized amplitude spectrum; 4, acquiring a joint detection window center vector, a line spectrum frequency vector and a line spectrum energy vector; 5, updating a line spectrum frequency vector and a line spectrum energy vector; 6, obtaining an effective harmonic order vector; 7, calculating the largest common divisor; 8, updating the extracted fundamental frequency; ninthly, judging whether a fundamental frequency search ending condition is met or not, if the condition is met, entering the tenth step, and otherwise, returning to the fourth step; tenth, judging whether a time frame traversal ending condition is met or not, if the condition is met, entering the eleventh step, otherwise, returning to the third step; and eleventhly, outputting the fundamental frequency extracted by each time frame. According to the method, effective detection of the harmonic families can be realized under the condition that one or more harmonic families exist.

Description

Technical Field

[0001] This invention belongs to the field of signal processing, specifically relating to a harmonic family detection method based on multi-window joint detection and harmonic feature evaluation. Background Technology

[0002] Signal characteristics are crucial for signal detection and analysis, encompassing features in the time domain, frequency domain, time-frequency domain, and fuzzy domain. Aerodynamic noise from UAV propellers, cavitation noise from ship propellers, and vital signs all exhibit harmonic characteristics in the frequency domain. This means their spectra contain multiple lines with multiples in a proportional relationship; these lines are called harmonic families, and the fundamental frequency of these families has practical physical significance. Detecting these harmonic families allows for UAV detection; extracting the fundamental frequency of these harmonic families can suppress the impact of ship propeller cavitation noise on received signals during combat, and is also significant for monitoring vital signs.

[0003] The current methods for detecting harmonic families mainly include: (1) methods based on traditional signal processing. Methods based on traditional signal processing mainly obtain the frequency domain information of the signal by methods such as fast Fourier transform, and determine the presence or absence of harmonic families by judging whether the line spectrum extracted from the frequency domain is a multiple relationship. The accuracy is high only under high signal-to-noise ratio conditions, and it is not applicable in complex situations such as the presence of line spectrum interference and multiple harmonic families; (2) methods based on machine learning and deep learning. Methods based on machine learning and deep learning usually use feature vectors such as Mel feature coefficients to represent the signal, and train feature vectors through a network to realize the binary classification problem of harmonic family detection. However, a large amount of data is required as the training set of the model, and the actual physical meaning is difficult to extract. Summary of the Invention

[0004] Technical Problem: The purpose of this invention is to provide a method for detecting harmonic families. This method can detect the presence and number of harmonic families, extract the fundamental frequency information of each harmonic family, and can also achieve effective detection of harmonic families under conditions of single line spectrum interference, multiple harmonic families, or no fundamental frequency.

[0005] Technical solution: To solve the above technical problems, such as Figure 1 As shown, this invention proposes a harmonic family detection method based on multi-window joint detection and harmonic feature evaluation, comprising the following steps:

[0006] (1) Obtain the data sequence x(n) to be processed;

[0007] (2) For each parameter of the detected harmonic family, including the fundamental frequency extracted from the p-th frame of the m-th time frame to be analyzed. Initialize the frame index m to be analyzed;

[0008] (3) Obtain the normalized amplitude spectrum of the m-th time frame to be analyzed in the data sequence x(n) to be processed. The fundamental frequency index k to be analyzed m and the extracted fundamental frequency index p;

[0009] (4) Obtain the kth frame of the m-th time frame to be analyzed. m The joint detection window center vector of the fundamental frequencies to be analyzed Joint detection of line spectrum frequency vector Joint detection of line spectrum energy vector

[0010] (5) Normalized amplitude spectrum based on the m-th time frame to be analyzed and the kth time frame to be analyzed in the m-th time frame m The joint detection window center vector of the fundamental frequencies to be analyzed Update the joint detection line spectrum frequency vector and joint detection of line spectrum energy vector

[0011] (6) Obtain the kth frame of the m-th time frame to be analyzed. m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed And update the joint detection line spectrum frequency vector and joint detection of line spectrum energy vector

[0012] (7) Based on the m-th time frame to be analyzed, the k-th time frame... m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Calculate the greatest common divisor

[0013] (8) Based on the m-th time frame to be analyzed, the k-th time frame... m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed Joint detection of line spectrum energy vector and the greatest common divisor Calculate the overall evaluation index of the harmonic family And update the fundamental frequency extracted from the p-th frame of the m-th time frame to be analyzed.

[0014] (9) Determine whether the fundamental frequency search termination condition is met. If the condition is met, proceed to step (10); otherwise, return to step (4).

[0015] (10) Determine whether the time frame traversal end condition is met. If the condition is met, proceed to step (11); otherwise, return to step (3).

[0016] (11) Output the fundamental frequency extracted from each time frame.

[0017] Further, in step (1), the following method is used to obtain the data sequence x(n) to be processed, which specifically includes the following steps:

[0018] Receive the real-time acquisition data of N sampling points from the sensor as the data sequence x(n) to be processed, n = 0, 1,..., N - 1, or extract the data of N sampling points from the memory as the data sequence x(n) to be processed, n = 0, 1,..., N - 1, where n is the index of the data sequence x(n) to be processed, and N is the number of sampling points corresponding to the pulse width length of the data sequence x(n) to be processed, and the value is N = 2 κ , where κ is a positive integer greater than or equal to 10.

[0019] Further, in step (2), the following method is used to initialize the parameters of each harmonic family to be detected, which specifically includes the following steps:

[0020] (2-1) Initialize the parameters for obtaining the normalized amplitude spectrum of the time frame to be analyzed, which specifically includes the initialization of the following parameters:

[0021] ① The length N of the sliding time window w Initialize to: 4 < N w < N / 4 and is a positive even number;

[0022] ② The step N of the sliding time window s Initialize to: 1 ≤ N s ≤ round{N w / 4}, where round{·} is the rounding function;

[0023] ③ The total number of time frames M total Initialize to: M total = floor{(N - N w ) / N s}-1, where floor{·} is the floor function;

[0024] (2-2) Initialize the parameters for obtaining the fundamental frequency search range, which specifically includes the initialization of the following parameters:

[0025] ① The lower limit of the fundamental frequency search Initialize to: a positive real number;

[0026] ② The upper limit of the fundamental frequency search Initialize to: a positive real number, where f s is the sampling frequency;

[0027] ③ The fundamental frequency search width f liberInitialize to: 10≤f liber Positive real numbers ≤50;

[0028] ④ Number of frames m when harmonic families have been detected detec Initialized to: m detec =0;

[0029] ⑤ Frame number threshold M for detecting harmonic families detec Initialize as: 1≤M detec ≤M total Positive integers;

[0030] (2-3) Initialize the parameters of the joint detection window, specifically including the initialization of the following parameters:

[0031] ① The harmonic order H of the joint detection is initialized to a positive integer of 5 ≤ H ≤ 10;

[0032] ② Joint detection window length vector W d Initialize to: Where h is the harmonic index to be analyzed. Let h be the detection window length for the h-th harmonic to be analyzed, satisfying To meet Positive integers;

[0033] (2-4) Initialize the parameters for calculating the local signal-to-noise ratio, specifically including the initialization of the following parameters:

[0034] ① Calculate the signal window length W for the local signal-to-noise ratio sig Initialize to: in, Let h be the signal window length of the h-th harmonic to be analyzed, satisfying To meet Positive integers;

[0035] ② Calculate the background noise window length W for the local signal-to-noise ratio noi Initialize to: in, Let h be the background noise window length for the h-th harmonic to be analyzed, satisfying To meet Positive integers;

[0036] ③ The parameter α1, which determines the relationship between background noise and the maximum values ​​on both sides, is initialized to a positive real number: 0 < α1 ≤ 2;

[0037] ④ The parameter α2, which determines the relationship between background noise and the minimum values ​​on both sides, is initialized to a positive real number of 1 ≤ α2 < 5;

[0038] (2-5) Initialize the parameters of the decision line spectrum, specifically including the initialization of the following parameters:

[0039] ① The line spectrum amplitude threshold A0 is initialized as a non-negative real number: 0 ≤ A0 < 1;

[0040] ② The local signal-to-noise ratio threshold ξ0 of the line spectrum is initialized as a non-negative real number: 0≤ξ0≤10;

[0041] ③ The limit number of line spectrum search attempts Y0 is initialized to a positive integer of 1 ≤ Y0 ≤ 10;

[0042] (2-6) Initialize the parameters for calculating the greatest common divisor, specifically including the initialization of the following parameters:

[0043] ① The frequency tolerance η1 for determining whether frequencies are the same is initialized as a non-negative real number: 0 ≤ η1 ≤ 20;

[0044] ② The quality factor calculation threshold η2 is initialized to a positive real number of 10 ≤ η2 ≤ 500;

[0045] ③ The quality factor calculation parameter η3 when it is less than the threshold is initialized to a positive real number of 20≤η3≤100;

[0046] ④ The quality factor calculation parameter η4 when it is greater than the threshold is initialized to a positive real number of η3≤η4≤200;

[0047] (2-7) Initialize the parameters for evaluating the harmonic family, specifically including the initialization of the following parameters:

[0048] ① The harmonic order stability weighting factor ω1 is initialized as a non-negative real number: 0≤ω1≤1;

[0049] ② The weighting factor ω2 of the greatest common divisor of harmonics is initialized as: 0≤ω2≤1-ω1, a non-negative real number;

[0050] ③ The harmonic energy weighting factor ω3 is initialized as: ω3 = 1 - ω1 - ω2, a non-negative real number;

[0051] ④ The overall evaluation threshold δ0 of the harmonic family is initialized as a positive real number where 0 < δ0 < 1;

[0052] (2-8) Initialize the parameters of the fundamental frequency extracted from each time frame, specifically including the initialization of the following parameters:

[0053] ① The maximum number of extractable base frequencies P in a single time frame is initialized to a positive integer of 1 ≤ P ≤ 10;

[0054] ②The fundamental frequency extracted from each time frame Initialize to: m is the frame index to be analyzed, m = 1, 2, ..., Mtotal M total The total number of time frames initialized in step (2-1), p is the extracted base frequency index, p = 1, 2, ..., P;

[0055] (2-9) Initialize other parameters, specifically including the initialization of the following parameters:

[0056] ① The frame index m to be analyzed is initialized to m = 1.

[0057] Furthermore, in step (3), the normalized amplitude spectrum of the m-th time frame to be analyzed in the data sequence x(n) to be processed is obtained using the following method. Specifically, the following steps are included:

[0058] (3-1) Obtain the segmented data sequence s of the m-th time frame to be analyzed. m (n seg ):

[0059] s m (n seg )=x((m-1)N s +n seg ),n seg =0,1,…,N w -1

[0060] Where m is the time frame index to be analyzed initialized in step (2-9), and n seg For segmented data sequence s m (n seg The index of N s and N w These are the sliding time window step and sliding time window length initialized in step (2-1), respectively;

[0061] (3-2) Obtain the m-th time frame segment data sequence s to be analyzed m (n seg The amplitude spectrum S m (l):

[0062]

[0063] Where l is the frequency index and j is the imaginary unit, i.e. |·| is the modulo function;

[0064] (3-3) Obtain the m-th time frame segment data sequence s to be analyzed. m (n seg Normalized amplitude spectrum

[0065]

[0066] Where max{·} is the function for finding the maximum value;

[0067] (3-4) Obtain the lower limit of the base frequency search and baseband search upper limit Corresponding frequency index and And calculate the upper limit of the fundamental frequency index of the m-th time frame to be analyzed.

[0068]

[0069] in, and These are the lower and upper limits of the fundamental frequency search initialized in step (2-2), respectively, Δf = f s / N w To normalize the resolution of the amplitude spectrum, f s Here, `round{·}` is the sampling frequency, and `round{·}` is the rounding function.

[0070] (3-5) Initialize the fundamental frequency index k of the m-th time frame to be analyzed. m And update the extracted baseband index p:

[0071] k m =1

[0072] p=1

[0073] Where, k m Used to iterate through each fundamental frequency between the lower and upper bounds of the fundamental frequency search. Corresponding to the kth m The frequency index of the base frequency.

[0074] Furthermore, in step (4), the k-th frame of the m-th time frame to be analyzed is obtained using the following method. m The joint detection window center vector of the fundamental frequencies to be analyzed Joint detection of line spectrum frequency vector Joint detection of line spectrum energy vector Specifically, the following steps are included:

[0075] (4-1) Obtain the k-th frame of the m-th time frame to be analyzed m The joint detection window center vector of the fundamental frequencies to be analyzed

[0076]

[0077] in, For the m-th time frame to be analyzed, the k-th time frame is... m The frequency index corresponding to the center of the detection window of the h-th harmonic to be analyzed for the fundamental frequency to be analyzed satisfies H is the joint detection harmonic order initialized in step (2-3), and The lower limit of the fundamental frequency search obtained in step (3-4) The corresponding frequency index;

[0078] (4-2) Initialize the k-th frame to be analyzed for the m-th time frame. m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed

[0079]

[0080] in, For the m-th time frame to be analyzed, the k-th time frame is... m The line spectrum frequency of the h-th harmonic to be analyzed from the fundamental frequency to be analyzed satisfies

[0081] (4-3) Initialize the kth frame of the m-th time frame to be analyzed m The joint detection line spectrum energy vector of the fundamental frequencies to be analyzed

[0082]

[0083] in, For the m-th time frame to be analyzed, the k-th time frame is... m The line spectrum energy of the h-th harmonic to be analyzed at the fundamental frequency satisfies

[0084] (4-4) Update the harmonic index h to be analyzed:

[0085] h = 1.

[0086] Furthermore, in step (5), the joint detection line spectrum frequency vector is updated using the following method. and joint detection of line spectrum energy vector Specifically, the following steps are included:

[0087] (5-1) Determine the k-th time frame to be analyzed for the m-th time frame. m The frequency index corresponding to the center of the detection window of the h-th harmonic to be analyzed for the fundamental frequency to be analyzed. Does the following overflow condition apply?

[0088]

[0089] in, For the m-th time frame to be analyzed obtained in step (4-1), the k-th time frame is... m The frequency index corresponding to the center of the detection window of the h-th harmonic to be analyzed for the fundamental frequency to be analyzed, N.w / 2 represents the normalized amplitude spectrum of the m-th time frame to be analyzed, obtained in step (3-3). The length; if the condition is met, proceed to step (6-1); otherwise, proceed to step (5-2);

[0090] (5-2) Obtain the k-th frame of the m-th time frame to be analyzed. m The search interval of the detection window for the h-th harmonic to be analyzed from the fundamental frequency to be analyzed.

[0091]

[0092] in, Let be the detection window length of the h-th harmonic to be analyzed initialized in step (2-3), max{·} and min{·} are the maximum value function and minimum value function respectively, and round{·} is the rounding function;

[0093] (5-3) Initialize the kth frame of the m-th time frame to be analyzed m Set of detection window cut-out points for the h-th harmonic to be analyzed from the fundamental frequency to be analyzed.

[0094]

[0095] Where Φ represents the empty set;

[0096] (5-4) Initialize the k-th frame to be analyzed for the m-th time frame. m Number of times the detection window spectrum search is performed for the h-th harmonic of the fundamental frequency to be analyzed.

[0097]

[0098] (5-5) Update the k-th frame of the m-th time frame to be analyzed. m The search interval of the detection window for the h-th harmonic to be analyzed from the fundamental frequency to be analyzed.

[0099]

[0100] Where - denotes the complement operation. Indicates by All that does not belong to A set of elements;

[0101] (5-6) Search for the k-th frame to be analyzed in the m-th time frame. m Normalized amplitude spectrum within the detection window of the h-th harmonic to be analyzed for the fundamental frequency. maximum value and its corresponding frequency index

[0102]

[0103] in, and They represent in Searching for normalized amplitude spectrum within range The maximum value and its corresponding frequency index;

[0104] (5-7)Judgment Does the following line spectrum determination amplitude condition meet:

[0105]

[0106] Where A0 is the line spectrum amplitude threshold initialized in step (2-5); if the condition is met, proceed to step (5-8); otherwise, proceed to step (5-22).

[0107] (5-8) Calculate the maximum index signal range

[0108]

[0109] Among them, l sig For signal interval frequency index, Let h be the signal window length of the h-th harmonic to be analyzed initialized in step (2-4), max{·} and min{·} are the maximum value function and minimum value function respectively, and round{·} is the rounding function;

[0110] (5-9) Calculate the maximum index average signal amplitude

[0111]

[0112] Where mean{·} is the mean function. It means that for all satisfying Normalized amplitude spectrum Take the average;

[0113] (5-10) Calculate the maximum index Left background noise range

[0114]

[0115] Among them, l noi This is a frequency index for the background noise range. and These are the signal window length and background noise window length of the h-th harmonic to be analyzed, initialized in steps (2-4), respectively;

[0116] (5-11) Calculate the maximum index Right side background noise range

[0117]

[0118] (5-12) Calculate the maximum index Background noise range

[0119]

[0120] Where ∪ represents the union of the left and right sets;

[0121] (5-13) Initialize the set of background noise interval removal points

[0122]

[0123] Where Φ represents the empty set;

[0124] (5-14) Update the background noise interval frequency index l noi In the background noise range The corresponding index v noi :

[0125] v noi =1;

[0126] (5-15) Update the background noise interval frequency index l noi :

[0127]

[0128] in, Indicates taking The vth noi One element;

[0129] (5-16) Determine the frequency index of the background noise interval l noi The corresponding normalized amplitude spectrum Does the following local energy condition apply?

[0130]

[0131] Where α1 and α2 are the parameters initialized in step (2-4) for judging the relationship between a certain background noise and the maximum and minimum values ​​on both sides, respectively, and N w / 2 represents the normalized amplitude spectrum of the m-th time frame to be analyzed, obtained in step (3-3). The length of the background noise interval; if the condition is met, then let the set of background noise interval removal points be... [l noi ] indicates a single element l noi If the set is formed, proceed to step (5-17); otherwise, proceed to step (5-17).

[0132] (5-17) Determine the frequency index of the background noise interval l noi In the background noise range The corresponding index v noi Does the following traversal termination condition apply?

[0133]

[0134] Where length{·} represents the length of the vector; if the condition is true, proceed to step (5-18); otherwise, let v noi =v noi +1, and return to step (5-15);

[0135] (5-18) Update the background noise range

[0136]

[0137] Where - denotes the complement operation. Indicates by All that does not belong to A set of elements;

[0138] (5-19) Calculate the maximum index mean amplitude of background noise

[0139]

[0140] Where mean{·} is the mean function. It means that for all satisfying Normalized amplitude spectrum Take the average;

[0141] (5-20) Calculate the maximum index Local signal-to-noise ratio

[0142]

[0143] in, and These are the mean signal amplitude and the mean background noise amplitude calculated in steps (5-9) and (5-19), respectively.

[0144] (5-21) Determine the local signal-to-noise ratio Does the following line spectrum decision local signal-to-noise ratio condition meet:

[0145]

[0146] Where ξ0 is the local signal-to-noise ratio threshold of the line spectrum initialized in step (2-5); if the condition is met, proceed to step (5-24); otherwise, proceed to step (5-22).

[0147] (5-22) Update the detection window cut-out point set

[0148]

[0149] in, and They represent taking The first and last elements;

[0150] (5-23) Determine the set of removal points in the detection window and number of line spectrum searches Does the following line spectrum search cutoff condition apply?

[0151] or

[0152] Where Y0 is the limit number of line spectrum searches initialized in step (2-5), and == indicates whether the elements contained in the left and right sets are completely identical; if the condition is true, proceed to step (5-26); otherwise, let Set the number of line spectrum searches for the detection window initialized in step (5-4), and return to step (5-5);

[0153] (5-24) Update the k-th frame of the m-th time frame to be analyzed. m The line spectrum frequency of the h-th harmonic to be analyzed from the fundamental frequency to be analyzed. and line spectrum energy

[0154]

[0155] Where, Δf=f s / N w To normalize the resolution of the amplitude spectrum, f s Where N is the sampling frequency. w The sliding time window length initialized in step (2-1), and These are the normalized amplitude spectra within the detection window calculated in steps (5-6). The maximum value and its corresponding frequency index;

[0156] (5-25) Update the k-th frame of the m-th time frame to be analyzed. m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed and line spectrum energy vector The h-th element and

[0157]

[0158]

[0159] (5-26) Determine whether the harmonic index h to be analyzed satisfies the following traversal termination condition:

[0160] h≥H

[0161] Where H is the number of harmonics detected in the joint detection initialized in step (2-3); if the condition is met, proceed to step (6-1); otherwise, let h = h + 1 and return to step (5-1).

[0162] Furthermore, in step (6), the k-th frame to be analyzed is obtained using the following method. m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed And update the joint detection line spectrum frequency vector and joint detection of line spectrum energy vector Specifically, the following steps are included:

[0163] (6-1) Obtain the k-th frame of the m-th time frame to be analyzed. m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed

[0164]

[0165] in, For the m-th time frame to be analyzed, the k-th time frame is... m The line spectrum frequency of the h-th harmonic to be analyzed from the fundamental frequency to be analyzed, where H is the joint detection harmonic order initialized in step (2-3). Indicates that m and k are fixed. m At that time, all satisfied The set consisting of the harmonic orders h;

[0166] (6-2) Update the k-th frame of the m-th time frame to be analyzed. m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed

[0167]

[0168] in, Indicates all that satisfy of The set that constitutes;

[0169] (6-3) Update the k-th frame of the m-th time frame to be analyzed. m The joint detection line spectrum energy vector of the fundamental frequencies to be analyzed

[0170]

[0171] in, Indicates all that satisfy of The set that constitutes;

[0172] (6-4) Calculation length

[0173]

[0174] Here, length{·} represents the length of the vector.

[0175] Furthermore, in step (7), the greatest common divisor is calculated using the following method. Specifically, the following steps are included:

[0176] (7-1)Judgement length Does the following condition for not extracting valid line spectra meet:

[0177]

[0178] in, For the calculation in step (6-4) The length of is equal to the length of the left and right numbers; if the condition is true, then let the length of the m-th frame to be analyzed be k-th. m The greatest common divisor extracted from the fundamental frequencies to be analyzed Proceed to step (8-1); otherwise, proceed to step (7-2).

[0179] (7-2)Judgement length Does the following condition, which allows for the extraction of only one valid spectral line, apply?

[0180]

[0181] Here, == indicates whether the two numbers on the left and right are the same; if the condition is true, then let the k-th number of the m-th frame to be analyzed be the same. mThe greatest common divisor extracted from the fundamental frequencies to be analyzed Indicates taking If the first element is found, proceed to step (8-1); otherwise, proceed to step (7-3).

[0182] (7-3) Find Difference frequency matrix

[0183]

[0184] in, and They represent taking The a-th and b-th elements, where a and b are the difference frequency matrices respectively. x-axis index and y-axis index,

[0185]

[0186] (7-4) Find Difference frequency vector

[0187]

[0188] Where sort{·} represents the ascending sort operation;

[0189] (7-5) Initialize the k-th frame to be analyzed for the m-th time frame. m The set of redundant difference frequencies after analyzing the fundamental frequency. Set of frequency occurrences of sum and difference

[0190]

[0191] Where Φ is the empty set;

[0192] (7-6) Obtain the end index of the difference frequency search and traversal index

[0193]

[0194]

[0195] Where length{·} represents the length of the vector;

[0196] (7-7) Obtain the same difference frequency vector R same :

[0197]

[0198] Where, r same The difference frequency vector elements, Indicates taking The There are 1 element, where η1 is the frequency tolerance initialized in step (2-6) to determine whether the frequencies are the same;

[0199] (7-8) Update the k-th frame of the m-th time frame to be analyzed. m The set of redundant difference frequencies after analyzing the fundamental frequency. Set of frequency occurrences of sum and difference

[0200]

[0201] Where ∪ represents the union of the left and right sets, and mean{·} is the mean function. Indicates to With a single element mean{R same The set composed of [mean{R] same Do the union of sets. Indicates to With a single element length{R same The set consisting of [length{R] same Do the union of sets;

[0202] (7-9) Update the difference frequency search traversal index

[0203]

[0204] (7-10)Judgment Does the following traversal termination condition apply?

[0205]

[0206] in, This is the end index for the difference frequency search. The end index of the difference frequency search obtained in step (7-6) is used; if the condition is met, proceed to step (7-11); otherwise, return to step (7-7).

[0207] (7-11) Obtain the k-th frame of the m-th time frame to be analyzed. m The set of quality factors of the fundamental frequency to be analyzed

[0208]

[0209] in, Indicates length is A vector of all zeros;

[0210] (7-12) Obtain the end index of the difference frequency search after redundancy removal and traversal index

[0211]

[0212] (7-13) Calculate the first The quality factor d of the difference frequency after redundancy removal liber :

[0213]

[0214] in, express The There are 10 elements, where η2, η3 and η4 are the quality factor calculation threshold, the quality factor calculation parameter when the value is less than the threshold and the quality factor calculation parameter when the value is greater than the threshold, respectively, initialized in step (2-6).

[0215] (7-14) Obtain the end index of the line spectrum frequency search and traversal index

[0216]

[0217]

[0218] in, For the m-th time frame to be analyzed calculated in step (6-4), the k-th time frame is... m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Length;

[0219] (7-15) Calculate the first The frequency of the first line spectrum and the first The ratio λ of the difference frequencies after redundancy removal:

[0220]

[0221] in, and They represent The elements and The One element;

[0222] (7-16) Update the set of quality factors The element

[0223]

[0224] in, Let d be the set of quality factors obtained in step (7-11). liber The first one calculated in step (7-13) The quality factor decision threshold of the difference frequency after redundancy removal, |·| is the modulus function, and round{·} is the rounding function;

[0225] (7-17) Determine the frequency of the line spectrum and search the traversal index. Does the following traversal termination condition apply?

[0226]

[0227] in, It is the end index of the line spectrum frequency search obtained in step (7-14); if the condition is met, proceed to step (7-18); otherwise, let And return to step (7-15);

[0228] (7-18) Determine the difference frequency search traversal index after redundancy removal. Does the following traversal termination condition apply?

[0229]

[0230] in, It is the end index of the difference frequency search after redundancy removal obtained in step (7-12); if the condition is met, proceed to step (7-19); otherwise, let And return to step (7-13);

[0231] (7-19) Extracting the greatest common divisor

[0232]

[0233] in, Indicates that m and k are fixed. m At that time, search for the set of quality factors The index corresponding to the maximum value.

[0234] Furthermore, in step (8), the overall evaluation index of the harmonic family is calculated using the following method. And update the fundamental frequency extracted from the p-th frame of the m-th time frame to be analyzed. Specifically, the following steps are included:

[0235] (8-1) Obtain the k-th frame of the m-th time frame to be analyzed. m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed difference sequence

[0236]

[0237] Where ∪ represents taking the union of the left and right sets, and dif{·} represents calculating the vector difference sequence. This indicates that the set [0] consisting of a single element 0 is first compared with... Perform a union, then differ the union. The elements in are The difference between adjacent elements;

[0238] (8-2) Calculate the stability evaluation index of harmonic order

[0239]

[0240] Where H is the joint detection harmonic order initialized in step (2-3), length{·} represents the length of the vector, and sum{·} represents the sum of the elements of the vector. Indicates to Decrease all elements by 1;

[0241] (8-3) Obtain the k-th frame of the m-th time frame to be analyzed. m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed difference sequence

[0242]

[0243] in, This indicates that the set [0] consisting of a single element 0 is first compared with... Perform a union, then differ the union. The elements in are The difference between adjacent elements;

[0244] (8-4) Calculate the greatest common divisor of harmonics and evaluate the index

[0245]

[0246] in, The greatest common divisor is the one calculated in step (7), mean{·} is the mean function, and floor{·} is the floor function. Indicates to Divide all elements by And round down;

[0247] (8-5) Calculate the harmonic energy evaluation index

[0248]

[0249] in, For the m-th time frame to be analyzed, the k-th time frame is... m The joint detection line spectrum energy vector of the fundamental frequency to be analyzed;

[0250] (8-6) Calculate the overall evaluation index of the harmonic family.

[0251]

[0252] Wherein, ω1, ω2 and ω3 are the harmonic order stability weight factor, harmonic greatest common divisor weight factor and harmonic energy weight factor initialized in step (2-7), respectively;

[0253] (8-7) Update the fundamental frequency index p extracted from the m-th time frame to be analyzed and the fundamental frequency extracted from the p-th time frame.

[0254]

[0255] Wherein, δ0 is the total evaluation threshold of the harmonic family initialized in step (2-7).

[0256] Furthermore, in step (9), the following method is used to determine whether the fundamental frequency search termination condition is met, specifically including the following steps:

[0257] Determine the fundamental frequency index k to be analyzed m Does the extracted fundamental frequency index p satisfy the following fundamental frequency search termination condition:

[0258] Or p > P

[0259] in, Let k be the upper limit of the fundamental frequency index of the m-th time frame to be analyzed calculated in step (3-4), and P be the maximum number of fundamental frequencies that can be extracted from a single time frame initialized in step (2-8). If the condition is met, proceed to step (10-1); otherwise, let k m =k m +1, and return to step (4-1).

[0260] Furthermore, in step (10), the following method is used to determine whether the time frame traversal end condition is met, specifically including the following steps:

[0261] (10-1) Update the number of frames m when the detected harmonic families are updated. detec :

[0262]

[0263] Where, m detec The number of time frames for the detected harmonic families initialized in step (2-2);

[0264] (10-2) Determine whether the time frame index m to be analyzed satisfies the following time frame traversal termination condition:

[0265] m≥M total

[0266] Among them, M total The total number of time frames initialized in step (2-1); if the condition is met, proceed to step (11); otherwise, proceed to step (10-3);

[0267] (10-3) Update the lower limit of the baseband search and baseband search upper limit

[0268]

[0269] in, f liber and M detec These are the lower limit, upper limit, and width limit of the fundamental frequency search, and the threshold for the number of frames to detect harmonic families initialized in step (2-2), respectively. m0 is the index of the detected frames, satisfying m0 = 1, 2, ..., m, p is the index of the extracted fundamental frequencies, satisfying p = 1, 2, ..., P, and P is the maximum number of fundamental frequencies that can be extracted in a single frame initialized in step (2-8). and Indicates in Under the conditions, search respectively The results corresponding to the minimum and maximum values;

[0270] (10-4) Update the time frame index m to be analyzed:

[0271] m = m + 1

[0272] Then return to step (3-1).

[0273] Furthermore, in step (11), the fundamental frequency extracted from each time frame is output using the following method.

[0274] and

[0275] Where m is the time frame index to be analyzed, M total The total number of time frames initialized in step (2-1), p is the extracted baseband index, and P is the maximum number of basebands that can be extracted in a single time frame initialized in step (2-8). The p-th fundamental frequency extracted from the m-th time frame to be analyzed is initialized in step (2-8) and updated in step (8-7). If its value is -1, it means that the value is invalid, that is, the fundamental frequency has not been extracted. and This indicates that m = 1, 2, ..., M total p = 1, 2, ..., P and Under the conditions

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

[0277] 1. This invention makes full use of the harmonic characteristics of the harmonic family and designs a joint detection window to achieve joint detection of each harmonic, as shown in steps 2 and 4. The designed joint detection window includes a window length and center that vary with the harmonic order. While ensuring that the extracted line spectrum has harmonic characteristics, it can also adapt to situations where there is a certain deviation between each harmonic frequency and an integer multiple of the fundamental frequency.

[0278] 2. This invention makes full use of the effective harmonic order, greatest common divisor and energy of the line spectrum for joint detection, and comprehensively evaluates the extracted harmonic family features. As shown in step 8, the three designed evaluation indicators can reduce the impact of noise on line spectrum detection and improve the accuracy of harmonic family detection and fundamental frequency extraction.

[0279] 3. This invention makes full use of the fundamental frequency extracted from the analyzed time frame, and realizes the self-shrinking of the fundamental frequency search range. As shown in step 10, the self-shrinking of the fundamental frequency search range reduces the amount of computation and improves the accuracy of harmonic family detection and fundamental frequency extraction. Attached Figure Description

[0280] Figure 1 This is a schematic flowchart of the method of the present invention;

[0281] Figure 2 The normalized amplitude spectra for each time frame in Example 1 are shown.

[0282] Figure 3 The fundamental frequency maps extracted from each time frame of Example 1;

[0283] Figure 4 The normalized amplitude spectra for each time frame in Example 2 are shown below.

[0284] Figure 5 This is the fundamental frequency map extracted from each time frame of Example 2. Detailed Implementation

[0285] To better understand the purpose, structure, and function of this invention, the invention will be further described below with reference to the accompanying drawings.

[0286] In the embodiments of the present invention, the simulated received signal model x(t) is:

[0287]

[0288] Among them, U count U represents the number of harmonic families. count For harmonic family index, H real h represents the number of harmonics in a single harmonic family. real For signal harmonic indexing; and They are respectively the uth count The amplitude vector, frequency vector, and initial phase vector of each harmonic family have a length of H. real ; and They are respectively the uth count The h-th harmonic family real The amplitude, frequency, and initial phase are given; τ0 is the signal start time, τ is the signal pulse width; and ο(t) is a signal with a mean of 0 and a variance of σ. 2 Gaussian white noise, variance σ 2 The magnitude depends on the signal-to-noise ratio (SNR): SNR = 10log 10 {A max / (2σ 2 )},A max It represents the maximum amplitude of all harmonics in all harmonic families.

[0289] With sampling frequency f s Discrete sampling is performed on the received signal x(t) from the above simulation to obtain the signal sampling data sequence x(n):

[0290]

[0291] Wherein, the signal starting point n0 = round{τ0f s}, the number of signal sampling points N = round{τf s}, round{·} is the rounding function.

[0292] Example 1:

[0293] The simulation signal parameters are set as follows: number of harmonic families U count =1, harmonic family index u count =1, the number of harmonics H of a single harmonic family real =5, signal harmonic index h real =1,2,…,5; Amplitude vector A of the first harmonic family 1 =[1,1,1,1,1], frequency vector f 1= [235,470,705,940,1175] and the initial phase vector Signal start time τ0 = 0s, signal pulse width τ = 5s, signal-to-noise ratio SNR = -5dB, sampling frequency f s =4kHz.

[0294] The following section estimates the parameters of the simulated signal:

[0295] According to step (1), extract the data of 20,000 sampling points containing the entire pulse signal from the memory as the data sequence x(n) to be processed;

[0296] Based on step (2), set the sliding time window length N. w =512, sliding time window step N s =128, total number of time frames M total =151; Lower bound of fundamental frequency search Baseband search limit Baseband search grace period f liber =30, number of frames m when harmonic families have been detected detec =0, frame number threshold M when detecting harmonic families detec =10; Joint detection harmonic order H=8, joint detection window length vector Calculate the signal window length for the local signal-to-noise ratio Background noise window length for calculating local signal-to-noise ratio

[40] Parameter α1 = 1 for judging the relationship between a certain background noise and the maximum values ​​on both sides; parameter α2 = 2.2 for judging the relationship between a certain background noise and the minimum values ​​on both sides; line spectrum amplitude threshold A0 = 0.1; line spectrum local signal-to-noise ratio threshold ξ0 = 5; line spectrum search limit number Y0 = 2; frequency tolerance η1 = 10 for judging whether the frequencies are the same; quality factor calculation threshold η2 = 100; quality factor calculation parameter η3 = 50 when less than the threshold; quality factor calculation parameter η4 = 100 when greater than the threshold; harmonic order stability weight factor ω1 = 0.7; harmonic greatest common divisor weight factor ω2 = 0.2; harmonic energy weight factor ω3 = 0.1; harmonic family total evaluation threshold δ0 = 0.65; maximum number of fundamental frequencies that can be extracted in a single time frame P = 5; time frame index to be analyzed m = 1;

[0297] Based on step (3), for the data sequence x(n), obtain the normalized amplitude spectrum of the m-th time frame segment data sequence to be analyzed. Get the lower limit frequency index of the baseband search and upper limit frequency index And calculate the upper limit of the fundamental frequency index of the m-th time frame to be analyzed. Initialize the fundamental frequency index k of the m-th time frame to be analyzed m=1, and update the extracted fundamental frequency index p=1; the normalized amplitude spectrum of each time frame is as follows Figure 2 As shown;

[0298] According to step (4), for the m-th time frame to be analyzed, the k-th time frame... m One fundamental frequency index to be analyzed, and the center vector of the joint detection window is obtained. Initialize the joint detection line spectrum frequency vector Joint detection of line spectrum energy vector And update the harmonic index to be analyzed, h = 1;

[0299] According to step (5), for the m-th time frame to be analyzed, the k-th time frame... m The frequency index corresponding to the center of the detection window of the h-th harmonic to be analyzed for the fundamental frequency to be analyzed. and the length of the detection window Get the search range of the detection window Line spectrum frequency and line spectrum energy And update the joint detection line spectrum frequency vector and line spectrum energy vector

[0300] According to step (6), for the m-th time frame to be analyzed, the k-th time frame... m The line spectrum frequency of the h-th harmonic to be analyzed from the fundamental frequency to be analyzed. Obtain the effective harmonic order vector of joint detection And update the joint detection line spectrum frequency vector Line spectrum energy vector and its length

[0301] According to step (7), for the m-th time frame to be analyzed, the k-th time frame... m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Calculate the greatest common divisor

[0302] According to step (8), for the kth frame of the m-th time frame to be analyzed... m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed Greatest common divisor and line spectrum energy vector Calculate the overall evaluation index of the harmonic family And update the extracted baseband.

[0303] According to step (9), for the fundamental frequency index k of the m-th time frame to be analyzed... m Based on the extracted fundamental frequency index p, determine whether the fundamental frequency search termination condition is met. If p > P, and the condition is true, proceed to step (10); otherwise, let k m =k m +1, and return to step (4);

[0304] Based on step (10), update the number of time frames m of the detected harmonic families. detec For the time frame index m to be analyzed, determine whether the time frame traversal termination condition m≥M is met. total If the condition is met, proceed to step (11); otherwise, update the lower limit of the baseband search. and baseband search upper limit Let m = m + 1, and return to step (3);

[0305] Based on step (11), output the fundamental frequency extracted from each time frame. and like Figure 3 As shown.

[0306] Example 2:

[0307] The simulation signal parameters are set as follows: number of harmonic families U count =2, harmonic family index u count =1,2, the number of harmonics H of a single harmonic family real =10, signal harmonic index h real =1,2,…,10; Amplitude vector of the first harmonic family Frequency vector f 1 =[320,640,960,1280,1600,1920,2240,2560,2880,3200], Initial phase vector The amplitude vector of the second harmonic family Frequency vector f 2 =[370,740,1110,1480,1850,2220,2590,2960,3330,3700], Initial phase vector Signal start time τ0 = 0s, signal pulse width τ = 3s, signal-to-noise ratio SNR = 3dB, sampling frequency f s =10kHz.

[0308] The following section estimates the parameters of the simulated signal:

[0309] According to step (1), extract the data of 30,000 sampling points containing the entire pulse signal from the memory as the data sequence x(n) to be processed;

[0310] Based on step (2), set the sliding time window length N. w =1280, sliding time window step N s=320, total number of time frames M total =88; Lower bound of baseband search Baseband search limit Baseband search grace period f liber =30, number of frames m when harmonic families have been detected detec =0, frame number threshold M when detecting harmonic families detec =10; Joint detection harmonic order H=8, joint detection window length vector Calculate the signal window length for the local signal-to-noise ratio Background noise window length for calculating local signal-to-noise ratio The parameters for determining the relationship between a certain background noise and its two maximum values ​​are α1 = 1, and the parameters for determining the relationship between a certain background noise and its two minimum values ​​are α2 = 2.2; the line spectrum amplitude threshold is A0 = 0.1, the line spectrum local signal-to-noise ratio threshold is ξ0 = 5, and the line spectrum search limit number is Y0 = 2; the frequency tolerance for determining whether frequencies are the same is η1 = 10, the quality factor calculation threshold is η2 = 100, the quality factor calculation parameter η3 = 50 when the threshold is less than the threshold, and the quality factor calculation parameter η4 = 100 when the threshold is greater than the threshold; the harmonic order stability weighting factor is ω1 = 0.7, the harmonic greatest common divisor weighting factor is ω2 = 0.2, the harmonic energy weighting factor is ω3 = 0.1, and the harmonic family total evaluation threshold is δ0 = 0.7; the maximum number of fundamental frequencies that can be extracted in a single time frame is P = 5; the time frame index to be analyzed is m = 1;

[0311] Based on step (3), for the data sequence x(n), obtain the normalized amplitude spectrum of the m-th time frame segment data sequence to be analyzed. Get the lower limit frequency index of the baseband search and upper limit frequency index And calculate the upper limit of the fundamental frequency index of the m-th time frame to be analyzed. Initialize the fundamental frequency index k of the m-th time frame to be analyzed m =1, and update the extracted fundamental frequency index p=1; the normalized amplitude spectrum of each time frame is as follows Figure 4 As shown;

[0312] According to step (4), for the m-th time frame to be analyzed, the k-th time frame... m One fundamental frequency index to be analyzed, and the center vector of the joint detection window is obtained. Initialize the joint detection line spectrum frequency vector Joint detection of line spectrum energy vector And update the harmonic index to be analyzed, h = 1;

[0313] According to step (5), for the m-th time frame to be analyzed, the k-th time frame... m The frequency index corresponding to the center of the detection window of the h-th harmonic to be analyzed for the fundamental frequency to be analyzed. and the length of the detection window Get the search range of the detection window Line spectrum frequency and line spectrum energy And update the joint detection line spectrum frequency vector and line spectrum energy vector

[0314] According to step (6), for the m-th time frame to be analyzed, the k-th time frame... m The line spectrum frequency of the h-th harmonic to be analyzed from the fundamental frequency to be analyzed. Obtain the effective harmonic order vector of joint detection And update the joint detection line spectrum frequency vector Line spectrum energy vector and its length

[0315] According to step (7), for the m-th time frame to be analyzed, the k-th time frame... m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Calculate the greatest common divisor

[0316] According to step (8), for the kth frame of the m-th time frame to be analyzed... m The joint detection effective harmonic order vector of the fundamental frequency to be analyzed Greatest common divisor and line spectrum energy vector Calculate the overall evaluation index of the harmonic family And update the extracted baseband.

[0317] According to step (9), for the fundamental frequency index k of the m-th time frame to be analyzed... m Based on the extracted fundamental frequency index p, determine whether the fundamental frequency search termination condition is met. If p > P, and the condition is true, proceed to step (10); otherwise, let k m =k m +1, and return to step (4);

[0318] Based on step (10), update the number of time frames m of the detected harmonic families. detec For the time frame index m to be analyzed, determine whether the time frame traversal termination condition m≥M is met. total If the condition is met, proceed to step (11); otherwise, update the lower limit of the baseband search. and baseband search upper limit Let m = m + 1, and return to step (3);

[0319] Based on step (11), output the fundamental frequency extracted from each time frame. and like Figure 5 As shown.

[0320] 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 harmonic family detection method based on multi-window joint detection and harmonic feature evaluation, characterized in that: The following steps are involved: (1) Obtain the data sequence x(n) to be processed; (2) Parameters of the detected harmonic family, including the fundamental frequency extracted from the pth time frame to be analyzed And the time frame index m to be analyzed, initialize; (3) Obtain the normalized amplitude spectrum of the mth time frame to be analyzed of the data sequence to be processed x(n) The base frequency index k to be analyzed m and the extracted fundamental frequency index p; (4) Get the kth time frame of the mth time frame to be analyzed m The joint detection window center vector of the fundamental frequencies to be analyzed Joint detection of line spectrum frequency vectors Joint detection of line spectrum energy vector (5) Based on the normalized amplitude spectrum of the mth time frame to be analyzed and the kth time frame to be analyzed m The joint detection window center vector of the fundamental frequencies to be analyzed Update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector (6) Get the kth time frame of the mth time frame to be analyzed m The effective harmonic order vector of the fundamental frequencies to be analyzed is And update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector (7) Based on the kth time frame to be analyzed m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Calculate the greatest common divisor (8) Based on the kth time frame to be analyzed m The effective harmonic order vector of the fundamental frequencies to be analyzed is Joint detection of line spectrum energy vector and the greatest common divisor Calculate the overall evaluation index of the harmonic family And update the base frequency extracted from the pth time frame to be analyzed (9) Determine whether the baseband search end condition is met. If the condition is met, proceed to step (10); otherwise, return to step (4); (10) Determine whether the time frame traversal end condition is met. If the condition is met, go to step (11); otherwise, return to step (3); (11) Output the base frequency extracted from each time frame 2. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 1 is characterized in that: In step (1), the following method is used to obtain the data sequence x(n) to be processed, which specifically includes the following steps: Receive real-time data collected from N sampling points from the sensor as the data sequence to be processed x(n), n=0,1,…,N-1, or extract data from N sampling points from the memory as the data sequence to be processed x(n), n=0,1,…,N-1, where n is the index of the data sequence to be processed x(n), and N is the number of sampling points corresponding to the pulse width length of the data sequence to be processed x(n), and the value is N=2 κ , κ is a positive integer greater than or equal to 10.

3. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 2 is characterized in that: In step (2), the following method is used to initialize the parameters of the harmonic family detection, which specifically includes the following steps: (2-1) Initialize the parameters for obtaining the normalized amplitude spectrum of the time frame to be analyzed, including the initialization of the following parameters: ①Sliding time window length N w Initialized as: 4 < N w Positive even number less than N / 4; ② Sliding time window step N s Initialization: 1≤N s ≤round{N w / 4}, where round{·} is a rounding function; ③Total number of time frames M total Initialized to: M total =floor{(NN w ) / N s }-1, where floor{·} is the floor rounding function; (2-2) Initializing the parameters for obtaining the baseband search range, specifically including the initialization of the following parameters: ① Base frequency search lower limit Initialized to: A positive real number of ; ② Base frequency search upper limit Initialized to: A positive real number, where f s is the sampling frequency; ③ Base frequency search margin f liber Initialized to: 10≤f liber A positive real number ≤50; ④ Number of time frames of detected harmonic groups m detec Initialized to: m detec =0; ⑤ Detection of harmonic group time frame number threshold M detec Initialization: 1≤M detec ≤M total A positive integer of ; (2-3) Initializing the parameters for obtaining the joint detection window, specifically including the initialization of the following parameters: ① The joint detection harmonic order H is initialized to a positive integer of 5≤H≤10; ②Joint detection window length vector W d Initialized to: Where h is the harmonic index to be analyzed, is the detection window length of the hth harmonic to be analyzed, satisfying h=1,2,…,H, To satisfy A positive integer of ; (2-4) Initialize the parameters for calculating the local signal-to-noise ratio, including the initialization of the following parameters: ① Calculate the signal window length W of the local signal-to-noise ratio sig Initialized to: in, is the signal window length of the hth harmonic to be analyzed, satisfying h=1,2,…,H, To satisfy A positive integer of ; ② Calculate the background noise window length W of the local signal-to-noise ratio noi Initialized to: in, is the background noise window length of the hth harmonic to be analyzed, satisfying h=1,2,…,H, To satisfy A positive integer of ; ③ The parameter α1 for determining the relationship between the background noise and the maximum value on both sides is initialized to a positive real number of 0<α1≤2; ④ The parameter α2 for determining the relationship between the background noise and the minimum value on both sides is initialized to a positive real number of 1≤α2<5; (2-5) Initialize the parameters of the decision line spectrum, including the initialization of the following parameters: ① The line spectrum amplitude threshold A0 is initialized to a non-negative real number: 0≤A0<1; ② The line spectrum local signal-to-noise ratio threshold ξ0 is initialized to a non-negative real number of 0≤ξ0≤10; ③ The number of line spectrum search limits Y0 is initialized to a positive integer of 1≤Y0≤10; (2-6) Initialize the parameters for calculating the greatest common divisor, including the initialization of the following parameters: ① The frequency tolerance η1 for judging whether the frequencies are the same is initialized to a non-negative real number of 0≤η1≤20; ② The quality factor calculation threshold η2 is initialized to a positive real number of 10≤η2≤500; ③ The quality factor calculation parameter η3 when it is less than the threshold is initialized to a positive real number of 20≤η3≤100; ④ The quality factor calculation parameter η4 when it is greater than the threshold is initialized to a positive real number: η3≤η4≤200; (2-7) Initialize the parameters for evaluating the harmonic family, including the initialization of the following parameters: ① The harmonic order stability weight factor ω1 is initialized to a non-negative real number of 0≤ω1≤1; ② The harmonic greatest common divisor weight factor ω2 is initialized to a non-negative real number: 0≤ω2≤1-ω1; ③ The harmonic energy weight factor ω3 is initialized to a non-negative real number: ω3 = 1-ω1-ω2; ④ The total evaluation threshold of the harmonic family δ0 is initialized to a positive real number of 0<δ0<1; (2-8) Initializing the parameters for storing the baseband extracted from each time frame, specifically including initializing the following parameters: ① The maximum number of extractable fundamental frequencies in a single time frame, P, is initialized to a positive integer of 1≤P≤10; ② Base frequency extracted from each time frame Initialized to: m is the index of the time frame to be analyzed, m=1,2,…,M total , M total is the total number of time frames initialized in step (2-1), p is the extracted base frequency index, p = 1, 2, ..., P; (2-9) Initialize other parameters, including the initialization of the following parameters: ① The frame index m to be analyzed is initialized to: m=1.

4. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 3 is characterized in that: In step (3), the normalized amplitude spectrum of the mth time frame to be analyzed of the data sequence to be processed x(n) is obtained by the following method: The specific steps include: (3-1) Obtain the segmented data sequence s of the mth time frame to be analyzed m (n seg ): Where m is the index of the time frame to be analyzed initialized in step (2-9), n seg is the segmented data sequence s m (n seg ), N s and N w are the sliding time window step and sliding time window length initialized in step (2-1) respectively; (3-2) Obtain the mth time frame segment data sequence s to be analyzed m (n seg )'s amplitude spectrum S m (l): Among them, l is the frequency index and j is the imaginary unit, that is, |·| is the modulo function; (3-3) Obtain the mth time frame segment data sequence s to be analyzed m (n seg ) Among them, max{·} is the maximum value function; (3-4) Get the base frequency search lower limit and fundamental frequency search upper limit The corresponding frequency index and And calculate the upper limit of the fundamental frequency index of the mth time frame to be analyzed in, and are the base frequency search lower limit and base frequency search upper limit initialized in step (2-2), Δf = f s / N w is the resolution of the normalized amplitude spectrum, f s is the sampling frequency, round{·} is the rounding function; (3-5) Initialize the base frequency index k of the mth time frame to be analyzed m , and update the extracted fundamental frequency index p: k m ＝1 p＝1 Among them, k m Used to traverse each fundamental frequency between the lower limit and the upper limit of the fundamental frequency search. Corresponding to the kth m The frequency index of the fundamental frequency.

5. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 4 is characterized in that: In step (4), the following method is used to obtain the kth time frame to be analyzed: m The joint detection window center vector of the fundamental frequencies to be analyzed Joint detection of line spectrum frequency vectors Joint detection of line spectrum energy vector The specific steps include: (4-1) Get the kth time frame of the mth time frame to be analyzed m The joint detection window center vector of the fundamental frequencies to be analyzed in, is the kth time frame to be analyzed m The frequency index corresponding to the center of the detection window of the hth harmonic of the fundamental frequency to be analyzed satisfies h=1,2,…,H, where H is the joint detection harmonic order initialized in step (2-3), and The lower limit of the fundamental frequency search obtained in step (3-4) The corresponding frequency index; (4-2) Initialize the kth time frame to be analyzed m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed in, is the kth time frame to be analyzed m The line spectrum frequency of the hth harmonic of the fundamental frequency to be analyzed satisfies h=1,2,…,H; (4-3) Initialize the kth time frame to be analyzed m The energy vector of the joint detection line spectrum of the fundamental frequencies to be analyzed in, is the kth time frame to be analyzed m The line spectrum energy of the hth harmonic of the fundamental frequency to be analyzed satisfies h=1,2,…,H; (4-4) Update the harmonic index h to be analyzed: h＝1。 6. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 5, characterized in that: In step (5), the following method is used to update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector The specific steps include: (5-1) Determine the kth time frame of the mth time frame to be analyzed m The frequency index corresponding to the center of the detection window of the hth harmonic of the fundamental frequency to be analyzed Whether the following overflow conditions are met: in, is the kth time frame of the mth time frame to be analyzed obtained in step (4-1) m The frequency index corresponding to the center of the detection window of the hth harmonic of the fundamental frequency to be analyzed, N w / 2 is the normalized amplitude spectrum of the mth time frame to be analyzed obtained in step (3-3) Length; if the condition is met, go to step (6-1); otherwise, go to step (5-2); (5-2) Get the kth time frame of the mth time frame to be analyzed m The detection window search interval of the hth harmonic to be analyzed of the fundamental frequency to be analyzed in, is the detection window length of the hth harmonic to be analyzed initialized in step (2-3), max{·} and min{·} are the maximum value function and the minimum value function respectively, and round{·} is the rounding function; (5-3) Initialize the kth time frame to be analyzed m The detection window removal point set of the hth harmonic of the fundamental frequency to be analyzed Among them, Φ represents the empty set; (5-4) Initialize the kth time frame to be analyzed m The number of line spectrum searches of the hth harmonic of the fundamental frequency to be analyzed (5-5) Update the kth time frame of the mth time frame to be analyzed m The detection window search interval of the hth harmonic to be analyzed of the fundamental frequency to be analyzed Among them, - represents the complement operation, Indicated by All that do not belong to A set of elements of (5-6) Search for the kth time frame to be analyzed m Normalized amplitude spectrum within the detection window of the hth harmonic of the fundamental frequency to be analyzed The maximum value of and its corresponding frequency index in, and Respectively expressed in Search normalized amplitude spectrum within range The maximum value of and its corresponding frequency index; (5-7)Judgment Whether the following line spectrum judgment amplitude conditions are met: Wherein, A0 is the line spectrum amplitude threshold initialized in step (2-5); if the condition is met, proceed to step (5-8); otherwise, proceed to step (5-22); (5-8) Calculate the maximum value index Signal range Among them, l sig is the signal interval frequency index, is the signal window length of the hth harmonic to be analyzed initialized in step (2-4), max{·} and min{·} are the maximum value function and the minimum value function respectively, and round{·} is the rounding function; (5-9) Calculate the maximum value index The mean signal amplitude Among them, mean{·} is the mean function, Satisfaction with all The normalized amplitude spectrum of Take the mean; (5-10) Calculate the maximum value index Left background noise area Among them, l noi is the frequency index of the background noise interval, and are respectively the signal window length and background noise window length of the hth harmonic to be analyzed initialized in step (2-4); (5-11) Calculate the maximum value index Right background noise zone (5-12) Calculate the maximum value index Background noise range Among them, ∪ means taking the union of the left and right sets; (5-13) Initialize the background noise interval removal point set Among them, Φ represents the empty set; (5-14) Update the background noise interval frequency index l noi In the background noise range The corresponding index v in noi : v noi ＝1； (5-15) Update the background noise interval frequency index l noi : in, Indicates taking No.v noi elements; (5-16) Determine the background noise interval frequency index l noi The corresponding normalized amplitude spectrum Are the following local energy conditions met: Among them, α1 and α2 are the parameters initialized in step (2-4) for judging the relationship between a certain background noise and the maximum and minimum values ​​on both sides, N w / 2 is the normalized amplitude spectrum of the mth time frame to be analyzed obtained in step (3-3) Length; if the condition is met, then the background noise interval is removed from the point set [l noi ] indicates that a single element l noi If the set is a set of , then go to step (5-17); otherwise, go to step (5-17); (5-17) Determine the background noise interval frequency index l noi In the background noise range The corresponding index v in noi Whether the following traversal end conditions are met: Where length{·} represents the length of the vector. If the condition is met, proceed to step (5-18). Otherwise, let v noi =v noi +1, and return to step (5-15); (5-18) Update background noise range Among them, - represents the complement operation, Indicated by All that do not belong to The set of elements of ; (5-19) Calculate the maximum value index The mean background noise amplitude Among them, mean{·} is the mean function, Satisfaction with all The normalized amplitude spectrum of Take the mean; (5-20) Calculate the maximum value index The local signal-to-noise ratio in, and are the signal amplitude mean and background noise amplitude mean calculated in step (5-9) and step (5-19), respectively; (5-21) Determine the local signal-to-noise ratio Whether the following line spectrum judgment local signal-to-noise ratio conditions are met: Wherein, ξ0 is the line spectrum local signal-to-noise ratio threshold initialized in step (2-5); if the condition is met, proceed to step (5-24); otherwise, proceed to step (5-22); (5-22) Update the detection window excavation point set in, and Respectively represent The first and last element of ; (5-23) Determine the set of points to be removed from the detection window Number of line searches Whether the following line spectrum search cutoff conditions are met: or Where, Y0 is the number of line spectrum searches initialized in step (2-5), == means judging whether the elements contained in the left and right sets are exactly the same; if the condition is met, go to step (5-26); otherwise, let The number of detection window line spectrum searches initialized in step (5-4) is returned to step (5-5); (5-24) Update the kth time frame to be analyzed m The line spectrum frequency of the hth harmonic of the fundamental frequency to be analyzed and line spectrum energy Where Δf = f s / N w is the resolution of the normalized amplitude spectrum, f s is the sampling frequency, N w is the length of the sliding time window initialized in step (2-1), and are the normalized amplitude spectra within the detection window calculated in steps (5-6) The maximum value of and its corresponding frequency index; (5-25) Update the kth time frame to be analyzed m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed and line spectrum energy vector The hth element of and (5-26) Determine whether the harmonic index h to be analyzed satisfies the following traversal end condition: h≥H Wherein, H is the joint detection harmonic order initialized in step (2-3); if the condition is met, proceed to step (6-1); otherwise, set h=h+1 and return to step (5-1).

7. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 6 is characterized in that: In step (6), the following method is used to obtain the kth time frame to be analyzed: m The effective harmonic order vector of the fundamental frequencies to be analyzed is And update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector The specific steps include: (6-1) Get the kth time frame of the mth time frame to be analyzed m The effective harmonic order vector of the fundamental frequencies to be analyzed is in, is the kth time frame to be analyzed m The line spectrum frequency of the hth harmonic to be analyzed of the fundamental frequency to be analyzed, H is the joint detection harmonic order initialized in step (2-3), Indicates fixed m and k m When all satisfied The set composed of harmonic order h; (6-2) Update the kth time frame to be analyzed m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed in, Indicates that all satisfaction of The collection composed of (6-3) Update the kth time frame of the mth time frame to be analyzed m The energy vector of the joint detection line spectrum of the fundamental frequencies to be analyzed in, Indicates that all satisfaction of The collection composed of (6-4) Calculation Length Wherein, length{·} represents the length of the vector.

8. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 7, characterized in that: In step (7), the greatest common divisor is calculated using the following method: The specific steps include: (7-1)Judgment Length Whether the following conditions are met for extracting valid line spectrum: in, is calculated in step (6-4) The length of , == means to judge whether the left and right numbers are the same; if the condition is met, let the kth number of the mth time frame to be analyzed m The greatest common divisor extracted from the fundamental frequencies to be analyzed And go to step (8-1); otherwise, go to step (7-2); (7-2)Judgement Length Whether the following conditions are met to extract only one valid line spectrum: Among them, == means to judge whether the left and right numbers are the same; if the condition is met, let the kth number of the mth time frame to be analyzed be m The greatest common divisor extracted from the fundamental frequencies to be analyzed Indicates taking , and go to step (8-1); otherwise, go to step (7-3); (7-3) Ask The difference frequency matrix in, and Respectively represent The a-th and b-th elements of The horizontal and vertical indices of (7-4) Ask The difference frequency vector Among them, sort{·} represents the ascending sort operation; (7-5) Initialize the kth time frame to be analyzed m The set of de-redundant difference frequencies of the base frequencies to be analyzed Sum and difference frequency occurrence counts Among them, Φ is an empty set; (7-6) Get the difference frequency search end index and traverse the index Among them, length{·} represents the length of the vector; (7-7) Obtain the same difference frequency vector R same : Among them, r same is the difference frequency vector Elements of Indicates taking No. elements, η1 is the frequency tolerance for determining whether the frequencies are the same as initialized in step (2-6); (7-8) Update the kth time frame to be analyzed m The set of de-redundant difference frequencies of the base frequencies to be analyzed Sum and difference frequency occurrence counts Among them, ∪ means taking the union of the left and right sets, mean{·} is the mean function, Express With a single element mean {R same } composed of a set [mean{R same }] do the union, Express With a single element length {R same } composed of a set [length{R same }]Do the union; (7-9) Update the difference frequency search traversal index (7-10)Judgment Whether the following traversal end conditions are met: in, The end index of the difference frequency search. is the difference frequency search end index obtained in step (7-6); if the condition is met, proceed to step (7-11); otherwise, return to step (7-7); (7-11) Get the kth time frame of the mth time frame to be analyzed m The quality factor set of the fundamental frequency to be analyzed in, Indicates the length All zero vector of ; (7-12) Get the end index of the difference frequency search after redundancy removal and traverse the index (7-13) Calculate the The quality factor decision threshold d of the difference frequency after redundancy removal liber : in, express No. elements, η2, η3 and η4 are respectively the quality factor calculation threshold initialized in step (2-6), the quality factor calculation parameter when it is less than the threshold and the quality factor calculation parameter when it is greater than the threshold; (7-14) Get the end index of the line spectrum frequency search and traverse the index in, is the kth time frame of the mth time frame to be analyzed calculated in step (6-4) m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed Length; (7-15) Calculate the The line frequencies and The ratio of the difference frequencies after redundancy removal λ: in, and Respectively No. elements and No. elements; (7-16) Update the quality factor set No. Elements in, is the quality factor set obtained in step (7-11), d liber is the first The quality factor decision threshold of the difference frequency after redundancy removal, |·| is the modulo function, and round{·} is the rounding function; (7-17) Determine the line spectrum frequency and search the traversal index Whether the following traversal end conditions are met: in, is the line spectrum frequency search end index obtained in step (7-14); if the condition is met, go to step (7-18); otherwise, let And return to step (7-15); (7-18) Determine the difference frequency search traversal index after de-redundancy Whether the following traversal end conditions are met: in, is the end index of the difference frequency search after redundancy removal obtained in step (7-12); if the condition is met, go to step (7-19); otherwise, let And return to step (7-13); (7-19) Extract the greatest common divisor in, Indicates fixed m and k m When searching for the quality factor set The index corresponding to the maximum value.

9. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 8, characterized in that: In step (8), the overall evaluation index of the harmonic family is calculated using the following method: And update the base frequency extracted from the pth time frame to be analyzed The specific steps include: (8-1) Get the kth time frame of the mth time frame to be analyzed m The effective harmonic order vector of the fundamental frequencies to be analyzed is The difference sequence Among them, ∪ means taking the union of the left and right sets, dif{·} means calculating the vector difference sequence, It means that first the set [0] consisting of a single element 0 and Do the union and then do the difference of the union, that is The elements in are The difference between adjacent elements in ; (8-2) Calculation of harmonic order stability evaluation index Where H is the joint detection harmonic order initialized in step (2-3), length{·} represents the length of the vector, and sum{·} represents the sum of the elements of the vector. Express All elements in are reduced by 1; (8-3) Get the kth time frame of the mth time frame to be analyzed m The joint detection line spectrum frequency vector of the fundamental frequencies to be analyzed The difference sequence in, It means that first the set [0] consisting of a single element 0 and Do the union and then do the difference of the union, that is The elements in are The difference between adjacent elements in ; (8-4) Calculate the greatest common divisor of harmonics evaluation index in, is the greatest common divisor calculated in step (7), mean{·} is the mean function, floor{·} is the floor function, Express All elements in are divided by and round down; (8-5) Calculation of harmonic energy evaluation index in, is the kth time frame to be analyzed m The joint detection line spectrum energy vector of the fundamental frequencies to be analyzed; (8-6) Calculate the overall evaluation index of the harmonic family Wherein, ω1, ω2 and ω3 are the harmonic order stability weight factor, the harmonic greatest common divisor weight factor and the harmonic energy weight factor initialized in step (2-7), respectively; (8-7) Update the base frequency index p extracted from the mth frame to be analyzed and the pth extracted base frequency Wherein, δ0 is the total evaluation threshold of the harmonic family initialized in step (2-7); In step (9), the following method is used to determine whether the baseband search end condition is met, which specifically includes the following steps: Determine the base frequency index k to be analyzed m And whether the extracted fundamental frequency index p satisfies the following fundamental frequency search end condition: or p>P in, is the upper limit of the fundamental frequency index of the mth time frame to be analyzed calculated in step (3-4), and P is the maximum number of fundamental frequencies that can be extracted from a single time frame initialized in step (2-8); if the condition is met, proceed to step (10-1); otherwise, let k m =k m +1, and return to step (4-1).

10. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 9, characterized in that: In step (10), the following method is used to determine whether the time frame traversal end condition is met, which specifically includes the following steps: (10-1) Update the number of time frames m of the detected harmonic group detec : Among them, m detec is the number of time frames of the detected harmonic group initialized in step (2-2); (10-2) Determine whether the time frame index m to be analyzed satisfies the following time frame traversal end conditions: m≥M total Among them, M total is the total number of time frames initialized in step (2-1); if the condition is met, proceed to step (11); otherwise, proceed to step (10-3); (10-3) Update the base frequency search lower limit and fundamental frequency search upper limit in, f liber and M detec are the fundamental frequency search lower limit, fundamental frequency search upper limit, fundamental frequency search margin and harmonic family detection time frame number threshold initialized in step (2-2), m0 is the detected time frame index, satisfying m0=1,2,…,m, p is the extracted fundamental frequency index, satisfying p=1,2,…,P, P is the maximum number of fundamental frequencies that can be extracted from a single time frame initialized in step (2-8), and Indicated in Under the conditions, search separately The results corresponding to the minimum and maximum values ​​of ; (10-4) Update the time frame index m to be analyzed: m=m+1 And return to step (3-1); In step (11), the base frequency extracted from each time frame is outputted using the following method: Among them, m is the index of the time frame to be analyzed, M total is the total number of time frames initialized in step (2-1), p is the extracted fundamental frequency index, and P is the maximum number of fundamental frequencies that can be extracted from a single time frame initialized in step (2-8); The pth fundamental frequency extracted from the mth time frame to be analyzed is initialized in step (2-8) and updated in step (8-7). If its value is -1, it means that the value is invalid, that is, the fundamental frequency is not extracted; Indicates that m=1,2,…,M total ,p=1,2,…,P and Conditional

Images