A harmonic family detection method based on multi-window joint detection and harmonic feature evaluation

Abstract

The present invention discloses a harmonic family detection method based on multi-window joint detection and harmonic feature evaluation, which includes the following steps: first, obtaining a data sequence; second, initializing parameters; third, obtaining a normalized amplitude spectrum; fourth, obtaining a joint detection window center vector, a line spectrum frequency vector, and a line spectrum energy vector; fifth, updating the line spectrum frequency vector and the line spectrum energy vector; sixth, obtaining an effective harmonic order vector; seventh, calculating a greatest common divisor; eighth, updating the extracted fundamental frequency; ninth, judging whether a fundamental frequency search end condition is met, if so, proceeding to step ten, otherwise, returning to step four; tenth, judging whether a time frame traversal end condition is met, if so, proceeding to step eleven, otherwise, returning to step three; and eleventh, outputting the fundamental frequency extracted from each time frame. The method can effectively detect harmonic families under the condition that there are single or multiple harmonic families.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of signal processing, and particularly relates to a harmonic family detection method based on multi-window joint detection and harmonic feature evaluation. BACKGROUND

[0002] Signal features are important basis for signal detection and analysis, including time domain, frequency domain, time-frequency domain and fuzzy domain features. The propeller aerodynamic noise of unmanned aerial vehicle, the propeller cavitation noise of ship and the vital sign signal all have harmonic features in the frequency domain, that is, there are multiple line spectra with multiple relationships in the frequency spectrum, which are called harmonic families, and the fundamental frequency of the harmonic family has actual physical meaning. By detecting the harmonic family, the unmanned aerial vehicle can be detected; by using the extracted fundamental frequency of the harmonic family, the influence of the propeller cavitation noise of the ship on the received signal in wartime can be suppressed, and it also has important significance in the monitoring of vital sign signals.

[0003] The current methods for detecting the harmonic family mainly include: (1) traditional signal processing-based method. The traditional signal processing-based method mainly obtains the frequency domain information of the signal based on the fast Fourier transform method, and judges whether the line spectrum extracted from the frequency domain has a multiple relationship to realize the judgment of the existence of the harmonic family. The correct rate is high only under the condition of high signal-to-noise ratio, and it is not applicable in the presence of line spectrum interference and multiple harmonic families; (2) machine learning and deep learning-based method. The machine learning and deep learning-based method usually uses mel feature coefficients and other feature vectors to represent the signal, and trains the feature vectors through the network to realize the binary classification problem of harmonic family detection, but a large amount of data is needed as the training set of the model, and the actual physical meaning is difficult to extract. SUMMARY

[0004] Technical problem: The present application aims to provide a harmonic family detection method which can detect the existence and number of harmonic families and extract the fundamental frequency information of each harmonic family, and can also realize effective detection of the harmonic family under the condition of single line spectrum interference, multiple harmonic families or absence of fundamental frequency.

[0005] Technical scheme: In order to solve the above technical problems, as shown in the accompanying drawings, the present application provides a harmonic family detection method based on multi-window joint detection and harmonic feature evaluation, which comprises the following steps: Figure 1

[0006] (1) obtaining a data sequence x(n) to be processed;

[0007] (2) initializing each parameter for detecting the harmonic family, including the mth extracted fundamental frequency of the pth analyzed time frame and the analyzed time frame index m;

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

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

[0010] (5) Based on the normalized amplitude spectrum of the mth time frame to be analyzed and the joint detection window center vector of the kth m fundamental frequency to be analyzed of the mth time frame to be analyzed update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector

[0011] (6) Obtain the joint detection effective harmonic number vector of the kth m fundamental frequency to be analyzed of the mth time frame to be analyzed and update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector

[0012] (7) Based on the joint detection line spectrum frequency vector of the kth m fundamental frequency to be analyzed of the mth time frame to be analyzed calculate the greatest common divisor

[0013] (8) Based on the joint detection effective harmonic number vector of the kth m fundamental frequency to be analyzed of the mth time frame to be analyzed the joint detection line spectrum energy vector and the greatest common divisor calculate the total evaluation index of the harmonic family and update the pth extracted fundamental frequency of the mth time frame to be analyzed

[0014] (9) Determine whether the fundamental frequency search end condition is met, if the condition is met, go 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, go to step (11); otherwise, return to step (3);

[0016] (11) Output the extracted fundamental frequencies of each time frame

[0017] Further, in step (1), the data sequence x(n) to be processed is obtained by the following method, specifically comprising the following steps:

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

[0019] Further, in step (2), the parameters of the detected harmonic family are initialized by the following method, specifically comprising the following steps:

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

[0021] ① the length N of the sliding time window w initialized as: 4<N w ; and

[0022] ② the step N of the sliding time window s initialized as: 1≤N s ≤round{N w / 4}, wherein round{·} is a rounding function;

[0023] ③ the total number M of time frames total initialized as: M total =floor{(N-N w ) / N s}-1, wherein floor{·} is a downward rounding function;

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

[0025] ① the lower limit of the fundamental frequency search initialized as a positive real number: ; and

[0026] ② the upper limit of the fundamental frequency search initialized as a positive real number: , wherein f s is the sampling frequency;

[0027] ③ the search margin f liberinitialized as: 10 ≤ f < 50, a positive real number; liber initialized as: 10 ≤ f < 50, a positive real number;

[0028] ④ The frame number m of the detected harmonic family detec initialized as: m = 0; detec detec initialized as: m = 0; detec

[0029] ⑤ The frame number threshold M of the detected harmonic family detec initialized as: 1 ≤ M < 50, a positive integer; detec initialized as: 1 ≤ M < 50, a positive integer; total initialized as: 1 ≤ M < 50, a positive integer;

[0030] (2-3) Initialize the parameters for obtaining the joint detection window, including the following parameter initialization:

[0031] ① The joint detection harmonic number H is initialized as: 5 ≤ H ≤ 10, a positive integer;

[0032] ② The joint detection window length vector W d initialized as: where h is the index of the harmonic to be analyzed, is the detection window length of the hth harmonic to be analyzed, satisfying is a positive integer satisfying

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

[0034] ① The signal window length W sig initialized as: where h is the index of the harmonic to be analyzed, is the signal window length of the hth harmonic to be analyzed, satisfying is a positive integer satisfying

[0035] ② The background noise window length W noi initialized as: where h is the index of the harmonic to be analyzed, is the background noise window length of the hth harmonic to be analyzed, satisfying is a positive integer satisfying

[0036] ③ The parameter α1 for judging the relationship between background noise and maximum value on both sides is initialized as: 0 < α1 ≤ 2, a positive real number;

[0037] ④ The parameter α2 for judging the relationship between background noise and minimum value on both sides is initialized as: 1 ≤ α2 < 5, a positive real number;

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

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

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

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

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

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

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

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

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

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

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

[0049] ② The harmonic greatest common divisor weight factor ω2 is initialized as a non-negative real number satisfying 0≤ω2≤1-ω1;

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

[0051] ④ The harmonic family total evaluation threshold δ0 is initialized as a positive real number satisfying 0<δ0<1;

[0052] (2-8) The parameters for saving the extracted fundamental frequencies of each time frame are initialized, specifically including the initialization of the following parameters:

[0053] ① The maximum number of extractable fundamental frequencies P for a single time frame is initialized as a positive integer satisfying 1≤P≤10;

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

[0055] (2-9) Initialize other parameters, 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 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:

[0058] (3-1) Obtain the segmented data sequence s of the mth 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 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;

[0061] (3-2) Get the mth time frame segmented data sequence s to be analyzed m (n seg )'s amplitude spectrum S m (l):

[0062]

[0063] Among them, l is the frequency index, j is the imaginary unit, that is, |·| is the modulo function;

[0064] (3-3) Get the mth time frame segmented data sequence s to be analyzed m (n seg )

[0065]

[0066] where max{•} is a maximum function;

[0067] (3-4) Obtain the lower limit of the fundamental frequency search and the upper limit of the fundamental frequency search corresponding to the frequency index and and calculate the upper limit of the fundamental frequency index of the mth time frame to be analyzed

[0068]

[0069] wherein, and are the lower limit and the upper limit of the fundamental frequency search initialized in step (2-2), respectively, and Δf = f s / N w is the resolution of the normalized amplitude spectrum, f s is the sampling frequency, and round{•} is a rounding function.

[0070] (3-5) Initialize the analysis fundamental frequency index k of the mth time frame to be analyzed m and update the extracted fundamental frequency index p:

[0071] k m = 1

[0072] p = 1

[0073] wherein k m is used to traverse each fundamental frequency between the lower limit and the upper limit of the fundamental frequency search, corresponding to the frequency index of the k m th fundamental frequency.

[0074] Further, in step (4), the following method is used to obtain the joint detection window center vector of the k m th analysis fundamental frequency of the mth time frame to be analyzed joint detection line spectrum frequency vector joint detection line spectrum energy vector Specifically, the following steps are included:

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

[0076]

[0077] wherein, is the frequency index corresponding to the detection window center of the hth analysis harmonic of the k m th analysis fundamental frequency of the mth time frame to be analyzed, satisfying H is the joint detection harmonic number initialized in step (2-3), and is the base frequency search lower limit obtained in step (3-4) corresponding to the frequency index;

[0078] (4-2) initializing the joint detection line spectrum frequency vector of the kth m analysis base frequency of the mth analysis time frame

[0079]

[0080] wherein, is the line spectrum frequency of the hth analysis harmonic of the kth m analysis base frequency of the mth analysis time frame, satisfying

[0081] (4-3) initializing the joint detection line spectrum energy vector of the kth m analysis base frequency of the mth analysis time frame

[0082]

[0083] wherein, is the line spectrum energy of the hth analysis harmonic of the kth m analysis base frequency of the mth analysis time frame, satisfying

[0084] (4-4) updating the analysis harmonic index h:

[0085] h = 1.

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

[0087] (5-1) judging whether the frequency index corresponding to the center of the detection window of the hth analysis harmonic of the kth m analysis base frequency of the mth analysis time frame satisfies the following overflow condition:

[0088]

[0089] wherein, is the frequency index corresponding to the center of the detection window of the hth analysis harmonic of the kth m analysis base frequency of the mth analysis time frame obtained in step (4-1), N​w / 2 is the normalized amplitude spectrum of the mth time frame to be analyzed obtained in step (3-3) ; if the condition is satisfied, go to step (6-1); otherwise, go to step (5-2);

[0090] (5-2) obtain the detection window search interval of the hth analyzed harmonic of the kth analyzed fundamental frequency of the mth time frame to be analyzed m

[0091]

[0092] wherein, is the detection window length of the hth analyzed harmonic 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;

[0093] (5-3) initialize the detection window hole set of the hth analyzed harmonic of the kth analyzed fundamental frequency of the mth time frame to be analyzed m

[0094]

[0095] wherein, Φ represents an empty set;

[0096] (5-4) initialize the number of spectrum search of the detection window of the hth analyzed harmonic of the kth analyzed fundamental frequency of the mth time frame to be analyzed m

[0097]

[0098] (5-5) update the detection window search interval of the hth analyzed harmonic of the kth analyzed fundamental frequency of the mth time frame to be analyzed m

[0099]

[0100] wherein, - represents the complement set operation, represents a set composed of all elements in that do not belong to ;

[0101] (5-6) search for the maximum value m of the normalized amplitude spectrum of the hth analyzed harmonic of the kth analyzed fundamental frequency of the mth time frame to be analyzed and the corresponding frequency index ​​​​​

[0102]

[0103] wherein, and respectively represent searching the maximum value of the normalized amplitude spectrum and its corresponding frequency index within the range of ;

[0104] (5-7) judging whether the following line spectrum decision amplitude condition is satisfied:

[0105]

[0106] wherein, A0 is the line spectrum amplitude threshold initialized in step (2-5); if the condition is satisfied, then entering step (5-8); otherwise, entering step (5-22);

[0107] (5-8) calculating the signal interval of the maximum value index

[0108]

[0109] wherein, 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 integer function;

[0110] (5-9) calculating the signal amplitude mean value of the maximum value index

[0111]

[0112] wherein, mean{·} is the mean value function, represents taking the mean value of all normalized amplitude spectra satisfying ;

[0113] (5-10) calculating the left background noise interval of the maximum value index

[0114]

[0115] wherein, l noi is the background noise interval frequency index, and ​The signal window length and the background noise window length of the hth harmonic to be analyzed initialized in step (2-4), respectively;

[0116] (5-11) Calculate the maximum value index Right background noise interval

[0117]

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

[0119]

[0120] Wherein, ∪ represents taking the union of the left and right two sets;

[0121] (5-13) Initialize the background noise interval dig-out point set

[0122]

[0123] Wherein, Φ represents the empty set;

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

[0125] v noi =1;

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

[0127]

[0128] Wherein, The v noi th element of ;

[0129] (5-16) Judge whether the normalized amplitude spectrum corresponding to the background noise interval frequency index l noi Satisfies the following local energy condition:

[0130]

[0131] Wherein, α1 and α2 are respectively the parameters for judging the relationship between the certain background noise and the maximum value and the minimum value on both sides initialized in step (2-4), and N w / 2 is the normalized amplitude spectrum of the mth time frame to be analyzed obtained in step (3-3)​ length{v } < length{v } then let the set of background noise interval carve-out points

[0132] noi} denotes the set consisting of a single element l noi , and go to step (5-17); otherwise, go to step (5-17);

[0132] (5-17) judge whether the background noise interval frequency index l noi in the background noise interval corresponding index v noi satisfies the following traversal end condition:

[0133]

[0134] wherein length{·} denotes the length of a vector; if the condition is satisfied, go to step (5-18); otherwise, let v noi = v noi + 1, and return to step (5-15);

[0135] (5-18) update the background noise interval

[0136]

[0137] wherein - denotes the complement set operation, denotes the set consisting of all elements in that do not belong to ;

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

[0139]

[0140] wherein mean{·} is the mean value function, denotes the mean value of all normalized amplitude spectrum satisfying ;

[0141] (5-20) calculate the local signal-to-noise ratio of the maximum value index

[0142]

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

[0144] (5-21) judging local signal-to-noise ratio whether the following line spectrum decision local signal-to-noise ratio condition is satisfied:

[0145]

[0146] wherein, ξ0 is the line spectrum local signal-to-noise ratio threshold initialized in step (2-5); if the condition is satisfied, go to step (5-24); otherwise, go to step (5-22);

[0147] (5-22) updating the detection window hole point set

[0148]

[0149] wherein, and respectively represent the 1st and last elements of ;

[0150] (5-23) judging the detection window hole point set and the line spectrum search number whether the following line spectrum search cut-off condition is satisfied:

[0151] or

[0152] wherein, Y0 is the line spectrum search limit number initialized in step (2-5), == = indicates whether the elements contained in the left and right sets are completely same; if the condition is satisfied, go to step (5-26); otherwise, let be the detection window line spectrum search number initialized in step (5-4), and return to step (5-5);

[0153] (5-24) updating the line spectrum frequency of the hth to-be-analyzed harmonic of the kth to-be-analyzed fundamental frequency of the mth to-be-analyzed time frame m and the line spectrum energy

[0154]

[0155] wherein, Δ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 spectrum in the detection window and the normalized amplitude spectrum in the detection window calculated in step (5-6) respectively ​The maximum value of and its corresponding frequency index;

[0156] (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

[0157]

[0158]

[0159] (5-26) Determine whether the harmonic index h to be analyzed meets the following traversal end conditions:

[0160] h≥H

[0161] 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).

[0162] Furthermore, in step (6), the following method is used to obtain the kth time frame to be analyzed: m The joint detection effective harmonic order vector of the fundamental frequencies to be analyzed And update the joint detection line spectrum frequency vector and the joint detection line spectrum energy vector The specific steps include:

[0163] (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

[0164]

[0165] in, is the kth time frame of the mth 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 number of joint detection harmonics initialized in step (2-3), Indicates fixed m and k m When all satisfied The set of harmonic orders h;

[0166] (6-2) Update 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

[0167]

[0168] wherein, denotes a set consisting of all satisfying ;

[0169] (6-3) updating the joint detection line spectrum energy vector of the kth m analysis fundamental frequency of the mth analysis time frame

[0170]

[0171] wherein, denotes a set consisting of all satisfying ;

[0172] (6-4) calculating the length of

[0173]

[0174] wherein, length{·} denotes the length of a vector.

[0175] Further, in step (7), the greatest common divisor is calculated by using the following method and specifically comprising the following steps:

[0176] (7-1) judging whether the length of satisfies the following condition that no valid line spectrum is extracted:

[0177]

[0178] wherein, is the length of calculated in step (6-4), == denotes judging whether the left and right two numbers are same; if the condition is satisfied, the greatest common divisor extracted by the kth m analysis fundamental frequency of the mth analysis time frame is set as and step (8-1) is entered; otherwise, step (7-2) is entered;

[0179] (7-2) judging whether the length of satisfies the following condition that only one valid line spectrum is extracted:

[0180]

[0181] wherein, == denotes judging whether the left and right two numbers are same; if the condition is satisfied, the greatest common divisor extracted by the kth m analysis fundamental frequency of the mth analysis time frame is set as and step (8-1) is entered; otherwise, step (7-2) is entered;the greatest common divisor extracted for the base frequency to be analyzed denotes the first element of , and step (8-1) is entered; otherwise, step (7-3) is entered;

[0182] (7-3) calculates the difference frequency matrix

[0183]

[0184] wherein, and denote the a-th and b-th elements of , and a and b are the horizontal coordinate index and the vertical coordinate index of the difference frequency matrix , respectively,

[0185]

[0186] (7-4) calculates the difference frequency vector

[0187]

[0188] wherein, sort{·} denotes the ascending order sorting operation;

[0189] (7-5) initializes the k-th m redundancy-removed difference frequency set of the m-th time frame to be analyzed for the base frequency to be analyzed and the difference frequency occurrence number set

[0190]

[0191] wherein, Φ is an empty set;

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

[0193]

[0194]

[0195] wherein, length{·} denotes the length of the vector;

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

[0197] ​​

[0198] wherein r same is the element of the difference frequency vector , denotes the k th element of the set , and η1 is the frequency tolerance for judging whether the frequencies are the same, which is initialized in step (2-6);

[0199] (7-8) updating the de-redundant difference frequency set and the difference frequency occurrence number set of the k m th to-be-analyzed fundamental frequency of the mth to-be-analyzed time frame

[0200]

[0201] wherein ∪ denotes taking the union set of the left and right sets, mean{·} is a mean value function, denotes taking the union set of and the set [mean{R same}] consisting of a single element mean{R same}, denotes taking the union set of and the set [length{R same}] consisting of a single element length{R same};

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

[0203]

[0204] (7-10) judging whether the following traversal end condition is met:

[0205]

[0206] wherein is the difference frequency search end index, and is the difference frequency search end index obtained in step (7-6); if the condition is met, step (7-11) is entered; otherwise, step (7-7) is returned;

[0207] (7-11) obtaining the quality factor set of the k m th to-be-analyzed fundamental frequency of the mth to-be-analyzed time frame

[0208]

[0209] wherein denotes the set of lengths ​​All zero vectors of ;

[0210] (7-12) Get the end index of the difference frequency search after de-redundancy and traverse the index

[0211]

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

[0213]

[0214] in, express No. elements, η2, η3 and η4 are 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 respectively;

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

[0216]

[0217]

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

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

[0220]

[0221] in, and Respectively No. elements and No. elements;

[0222] (7-16) Update the quality factor set No. Elements

[0223]

[0224] in, is the quality factor set obtained in steps (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;

[0225] (7-17) Determine the line spectrum frequency and search the traversal index Whether the following traversal end conditions are met:

[0226]

[0227] 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);

[0228] (7-18) Determine the difference frequency search traversal index after de-redundancy Whether the following traversal end conditions are met:

[0229]

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

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

[0232]

[0233] in, Indicates fixed m and k m When searching for the quality factor set The index corresponding to the maximum value.

[0234] Furthermore, in step (8), the following method is used to calculate the overall evaluation index of the harmonic family: And update the base frequency extracted from the pth time frame of the mth time frame to be analyzed The specific steps include:

[0235] (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

[0236]

[0237] where ∪ denotes taking the union of the left and right sets, and dif{·} denotes calculating the difference sequence of vectors, denotes first taking the union of the set [0] consisting of a single element 0 and , and then taking the difference of the union, i.e. each element in is the difference between adjacent elements in

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

[0239]

[0240] where H is the joint detection harmonic number initialized in step (2-3), length{·} denotes taking the length of a vector, and sum{·} denotes taking the sum of elements of a vector, denotes subtracting 1 from all elements in ;

[0241] (8-3) Obtain the difference sequence of the joint detection line spectrum frequency vector of the kth analysis fundamental frequency of the mth analysis time frame m

[0242]

[0243] where denotes first taking the union of the set [0] consisting of a single element 0 and , and then taking the difference of the union, i.e. each element in is the difference between adjacent elements in

[0244] (8-4) Calculate the harmonic greatest common divisor evaluation index

[0245]

[0246] where is the greatest common divisor calculated in step (7), mean{·} is the mean function, and floor{·} is the floor function, denotes dividing all elements in by and taking the floor;

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

[0248]

[0249] wherein, is the joint detection line spectrum energy vector of the k m th analysis fundamental frequency of the m th analysis time frame;

[0250] (8-6) calculating the total evaluation index of the harmonic family (8-7) updating the extracted fundamental frequency index p of the m th analysis time frame and the p

[0251] th extracted fundamental frequency

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

[0253] (8-7) updating the extracted fundamental frequency index p of the m th analysis time frame and the p th extracted fundamental frequency

[0254]

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

[0256] Further, in step (9), the following method is used to determine whether the fundamental frequency search end condition is met, which specifically includes the following steps:

[0257] determining whether the analysis fundamental frequency index k m and the extracted fundamental frequency index p meet the following fundamental frequency search end condition:

[0258] or p > P

[0259] wherein, is the upper limit of the fundamental frequency index of the m m th analysis time frame calculated in step (3-4), and P is the maximum extractable fundamental frequency number of a single time frame initialized in step (2-8); if the condition is met, step (10-1) is entered; otherwise, k m = k detec + 1, and step (4-1) is returned.

[0260] Further, 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:

[0261] (10-1) updating the number m detec of the detected harmonic family time frame:

[0262]

[0263] Among them, m detec is the number of time frames of the detected harmonic family initialized in step (2-2);

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

[0265] m≥M total

[0266] 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);

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

[0268]

[0269] 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 Indicates Under the conditions, search The results corresponding to the minimum and maximum values ​​of ;

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

[0271] m=m+1

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

[0273] Furthermore, in step (11), the following method is used to output the base frequency extracted from each time frame:

[0274] and

[0275] 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 for the mth time frame to be analyzed is initialized in step (2-8) and updated in step (8-7), and if the value is -1, it is an invalid value, i.e. the fundamental frequency is not extracted; and represents a value satisfying m=1, 2,..., M total , p=1, 2,..., P and

[0276] Advantages: Compared with the prior art, the technical scheme of the present application has the following advantages:

[0277] 1. The present application makes full use of the harmonic characteristics of the harmonic family, and designs a joint detection window to realize joint detection of each harmonic, as shown in steps 2 and 4. The designed joint detection window includes window length and center varying with the harmonic number. On the basis of ensuring that the subsequently extracted line spectrum has harmonic characteristics, it can also adapt to the case that the harmonic frequency and the integer multiple of the fundamental frequency exist certain deviations.

[0278] 2. The present application makes full use of the effective harmonic number, the greatest common divisor and the energy of the joint detection line spectrum, and comprehensively evaluates the extracted harmonic family characteristics, as shown in step 8. The designed three evaluation indexes can reduce the influence of noise and the like on line spectrum detection, and improve the correctness of harmonic family detection and fundamental frequency extraction.

[0279] 3. The present application makes full use of the extracted fundamental frequency of the analyzed time frame, and realizes self-constriction of the fundamental frequency search range, as shown in step 10. The self-constriction of the fundamental frequency search range reduces the calculation amount and improves the correctness of harmonic family detection and fundamental frequency extraction. BRIEF DESCRIPTION OF DRAWINGS

[0280] Figure 1 It is a flowchart of the method of the present application;

[0281] Figure 2 It is a normalized amplitude spectrum diagram of each time frame of Example 1;

[0282] Figure 3 It is a diagram of the extracted fundamental frequency of each time frame of Example 1;

[0283] Figure 4 It is a normalized amplitude spectrum diagram of each time frame of Example 2;

[0284] Figure 5 It is a diagram of the extracted fundamental frequency of each time frame of Example 2. DETAILED DESCRIPTION

[0285] In order to better understand the purpose, structure and function of the present application, the present application will be further described below in combination with the drawings.​

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

[0287]

[0288] Among them, U count is the number of harmonic groups, u count is the harmonic family index, H real is the number of harmonics in a single harmonic family, h real is the signal harmonic index; and u count The amplitude vector, frequency vector and initial phase vector of the harmonic family, with a length of H real ; and u count Harmonic group h real amplitude, frequency and initial phase; τ0 is the signal start time, τ is the signal pulse width; ο(t) is the mean value of 0 and the variance of σ 2 Gaussian white noise with variance σ 2 The size depends on the signal-to-noise ratio SNR: SNR = 10log 10 {A max / (2σ 2 )},A max is the maximum value of all harmonic amplitudes in all harmonic groups.

[0289] With sampling frequency f s The signal x(t) received in the above simulation is discretely sampled to obtain the signal sampling data sequence x(n):

[0290]

[0291] Among them, 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: number of harmonic groups U count =1, harmonic family index u count =1, the number of harmonics in a single harmonic family is H 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 initial phase vector Signal start time τ0= 0 s, signal pulse width τ = 5 s, signal-to-noise ratio SNR = -5 dB, sampling frequency f s = 4 kHz.

[0294] The simulation signal is then subjected to parameter estimation as follows:

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

[0296] According to step (2), the sliding time window length N w = 512, the sliding time window step N s = 128, the total number of time frames M total = 151; the lower limit of the fundamental frequency search the upper limit of the fundamental frequency search the fundamental frequency search width f liber = 30, the number of detected harmonic time frames m detec = 0, the threshold of the number of detected harmonic time frames M detec = 10; the number of harmonics to be jointly detected H = 8, the window length vector of the joint detection window the signal window length for calculating the local signal-to-noise ratio the background noise window length for calculating the local signal-to-noise ratio 40], the parameter α1 = 1 for judging the relationship between certain background noise and the maximum values on both sides, the parameter α2 = 2.2 for judging the relationship between certain background noise and the minimum values on both sides; the line spectrum amplitude threshold A0 = 0.1, the line spectrum local signal-to-noise ratio threshold ξ0 = 5, the line spectrum search limit number Y0 = 2; the frequency tolerance η1 = 10 for judging whether the frequencies are the same, the quality factor calculation threshold η2 = 100, the quality factor calculation parameter η3 = 50 when the threshold is less than the threshold, the quality factor calculation parameter η4 = 100 when the threshold is greater than the threshold; the harmonic number stability weight factor ω1 = 0.7, the harmonic greatest common divisor weight factor ω2 = 0.2, the harmonic energy weight factor ω3 = 0.1, the total evaluation threshold δ0 = 0.65 of the harmonic family; the maximum number of extractable fundamental frequencies P = 5 in a single time frame; the time frame index m = 1 to be analyzed;

[0297] According to step (3), for the data sequence x(n), the normalized amplitude spectrum of the mth time frame to be analyzed is obtained the lower limit frequency index of the fundamental frequency search is obtained and the upper limit frequency index is obtained the upper limit of the fundamental frequency index of the mth time frame to be analyzed is calculated the analysis fundamental frequency index k of the mth time frame to be analyzed is initialized 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 kth time frame to be analyzed m The fundamental frequency index to be analyzed is used to obtain the center vector of the joint detection window. Initialize the joint detection line spectrum frequency vector Joint detection of line spectrum energy vectors And update the harmonic index to be analyzed h=1;

[0299] According to step (5), for the kth time frame to be analyzed m The frequency index corresponding to the detection window center of the hth harmonic of the fundamental frequency to be analyzed and detection window length Get the detection window search range 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 kth time frame to be analyzed m The line spectrum frequency of the hth harmonic of the fundamental frequency to be analyzed Get 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 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

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

[0303] According to step (9), for the base frequency index k of the mth time frame to be analyzed, m and the extracted base frequency index p to determine whether the base frequency search end condition is met or p > P, if the condition is satisfied, go to step (10); otherwise, let k m = k m + 1, and return to step (4);

[0304] According to step (10), update the frame number m detec for the detected harmonic family, and judge whether the frame traversal end condition m ≥ M total is satisfied, if the condition is satisfied, go to step (11); otherwise, update the lower limit of the fundamental frequency search and the upper limit of the fundamental frequency search , and let m = m + 1, and return to step (3);

[0305] According to step (11), output the extracted fundamental frequency of each frame and as shown in Figure 3 .

[0306] Example 2:

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

[0308] The following is the parameter estimation of the simulation signal:

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

[0310] According to step (2), set the sliding time window length N w = 1280, and the sliding time window step N s=320, total number of time frames M total =88; fundamental frequency search lower limit Fundamental frequency search upper limit Fundamental frequency search margin f liber =30, number of detected harmonic family time frames m detec =0, the threshold of the number of frames for detecting harmonic groups is M detec =10; joint detection harmonic order H=8, joint detection window length vector Signal window length for calculating local signal-to-noise ratio Background noise window length for calculating local signal-to-noise ratio Parameter α1 for judging the relationship between a certain background noise and the maximum value on both sides = 1, parameter α2 for judging the relationship between a certain background noise and the minimum value on both sides = 2.2; 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 for judging whether the frequencies are the same = 10, 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.7; maximum number of extractable fundamental frequencies in a single time frame P = 5; time frame index m = 1;

[0311] According to step (3), for the data sequence x(n), obtain the normalized amplitude spectrum of the mth time frame segmented data sequence to be analyzed Get the fundamental frequency search lower limit frequency index and upper frequency index And calculate the upper limit of the fundamental frequency index of the mth time frame to be analyzed Initialize the base frequency index k of the mth 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 kth time frame to be analyzed m The fundamental frequency index to be analyzed is used to obtain the center vector of the joint detection window. Initialize the joint detection line spectrum frequency vector Joint detection of line spectrum energy vectors And update the harmonic index to be analyzed h=1;

[0313] According to step (5), for the kth time frame to be analyzed m The frequency index corresponding to the detection window center of the hth harmonic of the fundamental frequency to be analyzed and detection window length acquiring a detection window search interval line spectrum frequency and line spectrum energy and updating a joint detection line spectrum frequency vector and line spectrum energy vector

[0314] According to step (6), for the line spectrum frequency of the hth harmonic of the kth fundamental frequency to be analyzed of the mth time frame to be analyzed m , a joint detection effective harmonic number vector is acquired and a joint detection line spectrum frequency vector is updated line spectrum energy vector and its length

[0315] According to step (7), for the joint detection line spectrum frequency vector of the kth fundamental frequency to be analyzed of the mth time frame to be analyzed m , a greatest common divisor is calculated

[0316] According to step (8), for the joint detection effective harmonic number vector of the kth fundamental frequency to be analyzed of the mth time frame to be analyzed m , a greatest common divisor is calculated greatest common divisor and line spectrum energy vector a total evaluation index of the harmonic family is calculated and the extracted fundamental frequency is updated

[0317] According to step (9), for the fundamental frequency index k to be analyzed of the mth time frame to be analyzed m and the extracted fundamental frequency index p, it is determined whether a fundamental frequency search end condition is met or p > P, if the condition is met, step (10) is entered; otherwise, k m = k m + 1, and step (4) is returned to;

[0318] According to step (10), the number of time frames of the detected harmonic family m is updated detec , for the time frame index m to be analyzed, it is determined whether a time frame traversal end condition m ≥ M total is met, if the condition is met, step (11) is entered; otherwise, the lower limit of the fundamental frequency search and the upper limit of the fundamental frequency search are updated, and m = m + 1, and step (3) is returned to;

[0319] According to step (11), the extracted fundamental frequency of each time frame is outputted and​​ As Figure 5 shown.

[0320] It is to be understood that the present application is described by way of example only, and that modifications or alterations can be made to the features and embodiments described without departing from the spirit and scope of the application as set out in the claims. In addition, modifications can be made to the features and embodiments described to adapt them to particular situations and materials without departing from the spirit and scope of the application. Accordingly, the application is not limited to the specific embodiments disclosed herein, but rather, all such modifications as would be apparent to one skilled in the art are intended to be included within the scope of the application as defined in the claims.

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) The 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, initialized; (3) Obtain the normalized amplitude spectrum of the mth time frame to be analyzed of the data sequence to be processed x(n) The fundamental 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 vectors (5) Based on the normalized amplitude spectrum of the mth time frame to be analyzed and 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 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 joint detection effective harmonic order vector of the fundamental frequencies to be analyzed 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 joint detection effective harmonic order vector of the fundamental frequencies to be analyzed Joint detection of line spectrum energy vectors 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 of the mth time frame to be analyzed (9) Determine whether the base frequency 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, proceed to step (11); otherwise, return to step (3); (11) Output the base frequency extracted from each time frame In step (1), the data sequence x(n) to be processed is obtained by the following method, 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), which is N = 2 κ , κ is a positive integer greater than or equal to 10.

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 (2), the parameters for detecting the harmonic family are initialized using the following method, which specifically includes the following steps: (2-1) Initialize the parameters for obtaining the normalized amplitude spectrum of the time frame to be analyzed, specifically 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 function; (2-2) Initialize the parameters for obtaining the baseband search range, specifically including the initialization of the following parameters: ① Base frequency search lower limit Initialized as: A positive real number of ② Base frequency search upper limit Initialized as: A positive real number, where f s is the sampling frequency; ③ Base frequency search margin f liber Initialized to: 10≤f liber 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; (2-3) Initialize 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 as: Where h is the harmonic index to be analyzed, is the detection window length of the hth harmonic to be analyzed, satisfying 1,2,…,H, To satisfy A positive integer; (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 as: in, is the signal window length of the hth harmonic to be analyzed, satisfying To satisfy A positive integer; ② Calculate the background noise window length W of the local signal-to-noise ratio noi Initialized as: in, is the background noise window length of the hth harmonic to be analyzed, satisfying To satisfy A positive integer; ③ The parameter α1 that determines the relationship between the background noise and the maximum value on both sides is initialized to a positive real number: 0<α1≤2; ④ The parameter α2 used to determine the relationship between the background noise and the minimum value on both sides is initialized to a positive real number with a value 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: 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 determining whether the frequencies are the same is initialized to a non-negative real number: 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 the value 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: 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 δ0 of the harmonic family is initialized to a positive real number: 0<δ0<1; (2-8) Initialize the parameters for storing the fundamental frequencies extracted from each time frame, specifically including the initialization of 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; ②The fundamental frequency extracted from each time frame Initialized as: 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 fundamental frequency index, p = 1, 2, ..., P; (2-9) Initialize other parameters, including the following parameters: ① The frame index m to be analyzed is initialized to: m=1.

3. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 2, 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 ): s m (n seg )=x((m-1)N s +n seg ),n seg =0,1,…,N w -1 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) Get the mth time frame segmented data sequence s to be analyzed m (n seg )'s amplitude spectrum S m (l): Among them, l is the frequency index, j is the imaginary unit, that is, |·| is the modulo function; (3-3) Get the mth time frame segmented data sequence s to be analyzed m (n seg ) Among them, max{·} is the maximum value function; (3-4) Get the lower limit of the fundamental frequency search 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 They are the lower limit and upper limit of the fundamental frequency search 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.

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 (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 vectors 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 of the mth time frame to be analyzed m The frequency index corresponding to the detection window center of the hth harmonic of the fundamental frequency to be analyzed satisfies H is the number of joint detection harmonics 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 of the mth time frame to be analyzed m The line spectrum frequency of the hth harmonic of the fundamental frequency to be analyzed satisfies (4-3) Initialize the kth time frame to be analyzed m The joint detection line spectrum energy vector of the fundamental frequencies to be analyzed in, is the kth time frame of the mth time frame to be analyzed m The line spectrum energy of the hth harmonic of the fundamental frequency to be analyzed satisfies (4-4) Update the harmonic index h to be analyzed: h＝1。 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 (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 detection window center 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) If the condition is met, go to step 6; 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 of the fundamental frequency to be analyzed is 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 set of detection window removal points 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 Number of line spectrum searches in the detection window for the hth harmonic of the fundamental frequency to be analyzed (5-5) Update the kth time frame to be analyzed m The detection window search interval of the hth harmonic of the fundamental frequency to be analyzed is Among them, - represents the complement operation, Indicated by All that do not belong to A set of elements; (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 and its corresponding frequency index in, and Respectively expressed in Search normalized amplitude spectrum within the 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, Express 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 the signal window length and background noise window length of the hth harmonic to be analyzed initialized in steps (2-4) respectively; (5-11) Calculate the maximum value index Right background noise area (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 background noise interval frequency index l noi In the background noise range The corresponding index v in noi : v noi ＝1； (5-15) Update 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 steps (2-4) to judge 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) The length of; 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 A set of elements; (5-19) Calculate the maximum value index The mean background noise amplitude Among them, mean{·} is the mean function, Express 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 mean signal amplitude and the mean background noise amplitude calculated in step (5-9) and step (5-19), respectively; (5-21) Determining 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 elements of ; (5-23) Determine the detection window excavation point set 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 sliding time window length 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 meets the following traversal end conditions: h≥H Wherein, H is the joint detection harmonic order initialized in step (2-3); if the condition is met, proceed to step 6; otherwise, set h=h+1 and return to step (5-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 (6), the following method is used to obtain the kth time frame to be analyzed m The joint detection effective harmonic order vector of the fundamental frequencies to be analyzed 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 joint detection effective harmonic order vector of the fundamental frequencies to be analyzed in, is the kth time frame of the mth 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 number of joint detection harmonics initialized in step (2-3), Indicates fixed m and k m When all satisfied The set of harmonic orders h; (6-2) Update 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 in, Indicates all satisfaction of the collection it consists of; (6-3) Update the kth time frame of the mth time frame to be analyzed m The joint detection line spectrum energy vector of the fundamental frequencies to be analyzed in, Indicates all satisfaction of the collection it consists of; (6-4) Calculation Length Wherein, length{·} represents the length of the vector.

7. The harmonic family detection method based on multi-window joint detection and harmonic feature evaluation according to claim 6, 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 no valid line spectrum: in, is calculated in step (6-4) The length of , == means to judge whether the two numbers on the left and right are the same; if the condition is met, let the kth number of the mth time frame to be analyzed be the same. m The greatest common divisor extracted from the fundamental frequencies to be analyzed And go to step 8; 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 two numbers on the left and right are the same; if the condition is met, let the kth number of the mth time frame to be analyzed be the same. m The greatest common divisor extracted from the fundamental frequencies to be analyzed Indicates taking , and go to step 8; otherwise, go to step (7-3); (7-3) Ask The difference frequency matrix in, and Respectively represent The ath and bth elements of , a and b are the difference frequency matrices The horizontal and vertical indexes 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 A set of de-redundant difference frequencies of the base frequencies to be analyzed Sum and difference frequency occurrence counts Among them, Φ is the 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, 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 A set of de-redundant difference frequencies of the base frequencies to be analyzed Sum and difference frequency occurrence counts Among them, ∪ represents 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, is the end index of the difference frequency search, 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 A set of quality factors of the fundamental frequency to be analyzed in, Indicates the length is All zero vectors of ; (7-12) Get the end index of the difference frequency search after de-redundancy 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 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 respectively; (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 spectrum frequency 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 steps (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 de-redundancy 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.

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 (8), the following method is used to calculate the overall evaluation index of the harmonic family: And update the base frequency extracted from the pth time frame of the mth 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 joint detection effective harmonic order vector of the fundamental frequencies to be analyzed 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 is compared with Do the union and then do the difference of the union, that is The elements in 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 is compared with Do the union and then do the difference of the union, that is The elements in 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 of the mth time frame to be analyzed m The joint detection line spectrum energy vector of the fundamental frequency 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, harmonic greatest common divisor weight factor and harmonic energy weight factor initialized in steps (2-7), respectively; (8-7) Update the base frequency index p extracted from the mth time 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 base frequency search end condition is met, specifically including the following steps: Determine the fundamental frequency index k to be analyzed m And whether the extracted fundamental frequency index p meets the following fundamental frequency search end conditions: 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).

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 (10), the following method is used to determine whether the time frame traversal end condition is met, specifically including the following steps: (10-1) Update the number of time frames m of the detected harmonic family detec : Among them, m detec is the number of time frames of the detected harmonic family 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 lower limit of the fundamental frequency search 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 Indicates Under the conditions, search 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 following method is used to output the base frequency extracted from each time frame: 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