Musical instrument classification method based on time-frequency fine analysis

Abstract

This invention discloses a musical instrument classification method based on time-frequency fine analysis. The method includes: fundamental frequency estimation, frequency domain harmonic labeling and time-domain single-cycle fundamental frequency sequence estimation, time-frequency fine feature extraction, and a random forest model. This invention designs a cyclic iterative fundamental frequency estimation algorithm based on autocorrelation and narrowband spectral energy frequency estimation for fundamental frequency estimation. Then, a harmonic labeling algorithm based on adaptive window search and a time-domain single-cycle fundamental frequency sequence estimation algorithm are used to accurately label harmonic positions and obtain the single-cycle fundamental frequency sequence. Time-frequency fine features are used to characterize and extract the fine spectral structure and single-cycle fundamental frequency changes. Finally, the time-frequency fine features are fused with 13-dimensional Mel-frequency cepstral coefficient features and input into the random forest model for training. The musical instrument classification error rate based on the random forest model is reduced by 6% compared to using only Mel-frequency cepstral coefficient features.

Description

Technical Field

[0001] This invention relates to the field of signal processing technology, and specifically to a method for classifying musical instruments based on fine time-frequency analysis. Background Technology

[0002] By parsing audio files and obtaining the characteristics of musical instruments, instrument recognition can quickly and accurately identify the types of instruments contained within them, which is beneficial for its application in various fields such as intelligent music classification, instrument teaching, and music creation.

[0003] In past studies, scholars have achieved some classification results using time-frequency features proposed for instrument recognition. However, these features are often large-scale, describing characteristics across the entire time or frequency domain. For example, amplitude envelope describes the overall change in audio amplitude over a period of time, but it cannot accurately depict the fine changes between individual fundamental frequency periods. Similarly, the spectral centroid feature only characterizes the center of gravity of frequency domain energy from a macroscopic perspective of the entire spectrogram, failing to finely characterize the magnitude relationships between different harmonics and between a particular harmonic and nearby non-harmonics. In other words, the spectral centroid cannot accurately convey information about the spectral structure. However, spectral structure often plays a decisive role in the timbre of an instrument. Therefore, this invention argues that many past feature studies have neglected the refined representation of the time-frequency characteristics of music, a problem that urgently needs to be addressed.

[0004] The foundation of refined time-frequency analysis lies in the accurate estimation of the fundamental frequency. However, accurate estimation of the fundamental frequency and accurate labeling of harmonics are often complex processes, requiring the integration of intricate theories and consideration of anomalies in real-world audio data, such as aliasing, noise superposition, and spectral shift. This places extremely high demands on the adaptability and robustness of the fundamental frequency estimation and harmonic labeling algorithms. Therefore, implementing accurate fundamental frequency estimation and harmonic labeling in code to adapt to various real-world anomalies, thus laying a foundation for subsequent research, presents a considerable challenge. Summary of the Invention

[0005] The purpose of this invention is to address the aforementioned deficiencies in the prior art and provide a musical instrument classification method based on fine time-frequency analysis. First, the harmonics of the musical instrument audio are accurately estimated and the precise positions of each harmonic in its spectrum are marked, and the time-domain single-cycle fundamental frequency change sequence is calculated. Then, the fine time-frequency features of the musical instrument audio are extracted and fused with traditional audio features to classify the musical instrument.

[0006] The objective of this invention can be achieved by adopting the following technical solutions:

[0007] A musical instrument classification method based on time-frequency fine analysis, the musical instrument classification method comprising the following steps:

[0008] S1. Input a single-tone audio note of a musical instrument and accurately estimate the fundamental frequency of the audio.

[0009] S2. Perform FFT frequency domain harmonic marking based on the fundamental frequency of the audio.

[0010] S3. Estimate the time-domain single-cycle fundamental frequency sequence based on the fundamental frequency of the audio.

[0011] S4. Extract time-frequency fine features and Mel-frequency cepstral coefficient features;

[0012] S5. Based on the extracted video fine features and Mel frequency cepstral coefficient features, feature correlation and coefficient of variation are evaluated, and a random forest model is established to classify musical instruments.

[0013] Furthermore, step S1 is performed as follows:

[0014] S101. The short-time energy method is used to detect the single-tone audio endpoints of musical instruments and extract the valid audio. The formula for calculating short-time energy is as follows:

[0015]

[0016] In the formula, s[n] is the input single-note audio signal of the musical instrument, w[n] is the window function, n is the time-domain index of s[n] and w[n], and N is the time-domain index of s[n] and w[n]. f It is the frame length of the signal frame, z is the shift amount of the window function, and E z It is the short-time energy corresponding to a translation amount of z; E corresponds to different translation amounts. z A short-time energy sequence is formed; a threshold is set for the short-time energy of the effective audio start and end points of the instrument's single-tone audio. When the short-time energy of the current frame exceeds the threshold of the start point, the midpoint of the current frame is recorded as the start point of the effective audio. When the short-time energy of the current frame is lower than the threshold of the end point, the midpoint of the current frame is recorded as the end point of the effective audio.

[0017] S102. The fundamental frequency f0 of the effective audio x[n] is roughly estimated by time-domain autocorrelation. The autocorrelation formula for discrete signals is as follows:

[0018]

[0019] In the formula, τ represents the lag value, i.e. the shift amount. After shifting x[n] to the left by τ units, it is multiplied by the original signal. Finally, the product results are summed to obtain the autocorrelation value under the current lag value. The autocorrelation value is the largest when the lag value is zero. Based on the relationship between the autocorrelation value and different lag values, the maximum point other than the lag value of zero is found. The difference between the lag value of the maximum point and the point where the lag value is zero is the approximate estimate of the fundamental frequency period, thereby estimating the fundamental frequency f0.

[0020] S103. Perform a Fast Fourier Transform (FFT) on the audio signal and record the frequency f corresponding to the current FFT maximum value. max ; Determine f based on time-domain autocorrelation results max Is it an interference component? If so, remove the component from the FFT spectrum and then re-obtain the frequency f corresponding to the maximum value of the FFT. max Repeat the above process until, based on the time-domain autocorrelation results, it is determined that f... max It is not an interfering component;

[0021] S104. Calculate the denoised frequency value f of the frequency corresponding to the maximum peak value of the Fast Fourier Transform using narrowband spectrum energy frequency estimation. c To eliminate noise interference, the denoised frequency value f of the FFT spectrum is calculated using narrowband spectral energy frequency estimation when the index of the FFT spectrum array starts from 1. c The formula is as follows:

[0022]

[0023] Where k0 is the FFT index corresponding to the maximum peak value of the Fast Fourier Transform, k is the index of the FFT spectrum, and |X k | is the FFT modulus corresponding to the k-th point in the FFT spectrum, F s N is the sampling rate of the audio signal, and N is the number of FFT points.

[0024] S105, Based on the noise reduction frequency value f c Based on the fundamental frequency f0 estimated by autocorrelation, the maximum peak value is calculated to be approximately located at n. m At the subharmonic level; where:

[0025]

[0026] The fundamental frequency estimate is obtained as f1:

[0027]

[0028] S106. Resample the audio time-domain signal based on the fundamental frequency estimate f1, and set the new sampling rate f. s The calculation formula is as follows:

[0029] f s =f1×2 m Formula (6)

[0030] The value of m must satisfy:

[0031]

[0032] To ensure that the new sampling rate is not too large, which would prevent the number of sampling points N from covering a sufficient fundamental frequency period, the resampled time-domain signal is executed again from step S102 until the fundamental frequency estimate f1 changes by less than 0.001 compared to the fundamental frequency estimate of the previous cycle or the number of cycles reaches the limit of 50, and the final fundamental frequency estimate f1 is output.

[0033] Furthermore, in step S2, the FFT frequency domain harmonic marking process based on the fundamental frequency of the audio is as follows:

[0034] S201. Use the short-time energy method to detect the audio endpoints of musical instruments and extract the valid audio portion;

[0035] S202. Based on the fundamental frequency estimate f1 from step S1, the audio time-domain signal is resampled to make the audio sampling rate a power of f1, thereby reducing spectral leakage and converting the FFT spectrum into a logarithmic spectrum.

[0036] S203. Determine the left endpoint l and right endpoint r of the next harmonic peak search range on the FFT spectrum. If it is the first 3 harmonics, the calculation formula is as follows:

[0037]

[0038]

[0039] Where P is the position index of the previous harmonic peak in the FFT spectrum. If the position of the first harmonic peak is marked, then P = 1.

[0040] If the 4th harmonic and subsequent harmonics are being marked, the formulas for calculating the left endpoint l and the right endpoint r are as follows:

[0041] l=P+bT-2 Formula (10)

[0042] r = P + 1.09T1 Formula (11)

[0043] Where bT is half the distance between the third harmonic peak and the first harmonic peak on the FFT spectrum; T1 is the FFT point interval between the two harmonics closest to the currently marked harmonic on the FFT spectrum;

[0044] S204. Find all maximum points within the search range, and take the maximum value of all maximum points as the desired peak point of the current harmonic. If there are no maximum points within the search range, expand the left and right boundaries of the search range when marking the first four harmonics, and only expand the right boundary when marking other harmonics. When expanding the search boundary of the first two harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions:

[0045]

[0046]

[0047] When expanding the search boundaries for the 3rd and 4th harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions:

[0048] lP≥0.7T1 formula (14)

[0049] rP≤1.3T1 formula (15)

[0050] When expanding the search boundary for the 5th and higher harmonics, the right endpoint r must satisfy the following condition:

[0051] rP≤1.3T1 formula (16)

[0052] Then, within the new search range, the maximum value of all maxima is found as the peak value of the current harmonic. All FFT logarithmic amplitude values ​​and the index value of the maximum value within 0.1 times the harmonic interval centered on the maximum value are recorded and stored in matrix HWB and sequence HWI, respectively. Each row of HWB stores the peak-to-peak value of each harmonic and the spectral logarithmic value of its vicinity. If a maximum value is still not found, the search range is expanded and the process is repeated.

[0053] Furthermore, in step S3, the process of estimating the time-domain single-cycle fundamental frequency sequence based on the fundamental frequency of the audio is as follows:

[0054] S301. Use the short-time energy method to detect the audio endpoints of musical instruments and extract the valid audio portion;

[0055] S302. Use a high-pass filter to filter out low-frequency interference;

[0056] S303. Use cross-correlation to mark the boundary point of the single fundamental frequency period. Taking a point o in the middle of the audio signal as the origin, extend equidistantly to both sides. The extension distance is the length of the single fundamental frequency period obtained in the previous calculation. The signal within the extension range is a shift vector. The length of the shift vector is twice the length of the single fundamental frequency period obtained in the previous calculation, thus having the characteristic of adaptive adjustment. If it is the first cross-correlation calculation, the single frequency period length refers to the fundamental frequency estimation result in step S2. Taking point o as the origin, take four signals of the single fundamental frequency period length to the left or right as a fixed direction. The quantity, where the single fundamental frequency period length is also the single fundamental frequency period length obtained in the previous calculation; the shift vector is used as a sliding window and cross-correlation operation is performed with the fixed vector to obtain the correspondence between the cross-correlation value and the relative displacement; the extreme point sequence of the cross-correlation value and the corresponding relative displacement sequence are recorded from the correspondence; reasonable values ​​are selected from the relative displacement sequence, and the next single fundamental frequency period dividing point is calculated with point o as the reference. The distance between the dividing point and o is the single fundamental frequency period length; this process is repeated to obtain the time-domain single-cycle fundamental frequency sequence. The calculation formula for the cross-correlation r[τ] is as follows:

[0057]

[0058] In the formula, τ represents the lag value, i.e. the shift amount. g[n] and v[n] are two different signals. After shifting v[n] to the left by a distance of τ units, it is multiplied by g[n]. The sum of the product results is the similarity under the current lag value.

[0059] Furthermore, in step S4, the process of extracting the time-frequency fine features and the Mel frequency cepstral coefficient features is as follows:

[0060] S401. Extract fine frequency domain features, including: centroid of harmonic energy with logarithmic values ​​greater than zero, average attenuation ratio of harmonic peaks and spikes, maximum harmonic order with logarithmic values ​​greater than zero, slope-bandwidth ratio between two harmonics, harmonic peak-to-valley ratio, harmonic energy decrease rate, energy ratio of odd to even harmonics in the first 6 harmonics, peak-to-valley ratio of the first harmonic, decrease rate of slope-to-bandwidth ratio between harmonics, ratio of harmonic energy with logarithmic values ​​greater than zero to total energy, harmonic peak-to-valley ratio to half-harmonic spacing, and peak-to-valley ratio to half-harmonic spacing of the first harmonic.

[0061] From the HWB obtained in step S2, the logarithmic amplitude sequence P[n] of the maximum value points of each harmonic peak can be obtained. h The numbers less than or equal to zero in the logarithmic amplitude array are set to zero, with indices starting from 1. The formula for calculating the centroid c of harmonic energy with a logarithmic value greater than zero is as follows:

[0062]

[0063] Where L P It is P[n hThe length of ] is also the maximum harmonic order, n h It is the current harmonic order, and also P[n h The index of ];

[0064] The maximum harmonic order with a logarithmic value greater than zero is the logarithmic amplitude array P[n]. h The largest index of the non-zero value corresponding to the less than or equal to zero numbers in the array after setting them to zero, with the index starting from 1;

[0065] The calculation of the average attenuation ratio of harmonic peak glitches requires that, in step S2, the logarithmic spectrum of the FFT spectrum is greater than zero, and that the maxima A[n] of the twenty non-harmonic peaks closest to the maximum harmonic peak within a harmonic spacing range centered on the maximum harmonic peak value be obtained. p ], n p It is A[n p If the number of maxima within the range is less than twenty, then all maxima are taken; if the number of maxima is less than or equal to 10, then A[n] is taken directly. p The mean of A[n] is used as the average amplitude of the spikes around the harmonic peak; if the number of maxima is greater than 10, then in A[n] p Find a continuous set of maxima of length 10 in the FFT index, where the FFT index span of this set of maxima either contains the index of the maximum harmonic peak or there are no other maxima between it and the index of the maximum harmonic peak, and the sum is the largest. Calculate the mean of this set of maxima as the average amplitude of the spikes around the harmonic peak. Let the average amplitude of the spikes around the harmonic peak be 'a'. The formula for the harmonic peak attenuation ratio (PDR) of a single harmonic is as follows:

[0066]

[0067] Where p is the maximum amplitude of the current harmonic peak, and the average of the harmonic peak attenuation ratios of multiple harmonics is the average harmonic peak attenuation ratio.

[0068] Let the logarithmic spectrum of the FFT spectrum in step S2 be mag[k], and the HWI obtained in step S2 be the index sequence of the maximum harmonic peaks in mag[k], abbreviated as hi[n]. h If all the indices of the sequence start from 1, then the formula for calculating the slope bandwidth ratio (SFWR) between the 3rd and 4th harmonics is:

[0069]

[0070] When calculating the interharmonic slope bandwidth ratio descent rate, first calculate the interharmonic slope bandwidth ratio between the DC component and the first harmonic in the logarithmic spectrum as r1. Then calculate the interharmonic slope bandwidth ratio for each harmonic interval between the 5th and 8th harmonics and take the average as r2. If there are less than 8 harmonics, then take the highest harmonic. The formula for calculating the interharmonic slope bandwidth ratio descent rate SFWRD is as follows:

[0071]

[0072] The harmonic peak-to-valley ratio is calculated based on the portion of the logarithmic spectrum from which the logarithm is greater than zero using the FFT. All elements less than or equal to zero in mag[k] are set to zero, resulting in a new sequence MAG[k]. For a harmonic peak, there are two troughs, left and right, with the left trough value being v. p1 Its value is the mean of MAG[k] within the trough range of the left harmonic spacing; the trough value on the right is v. p2 Its value is the average of MAG[k] within the trough range of the right harmonic spacing; in the calculation of the peak-to-trough ratio of the Mth harmonic, v p1 The calculation formula is as follows:

[0073]

[0074] Where L is the left boundary for calculating the trough value, and R is the right boundary for calculating the trough value;

[0075]

[0076]

[0077] The formula for calculating v2 in the peak-to-valley ratio of the Mth harmonic is as follows:

[0078]

[0079]

[0080]

[0081] Given that the logarithmic amplitude sequence of the peak maximum points of each harmonic is P[n] h The formula for calculating the peak-to-valley ratio (MR) of the Mth harmonic is as follows:

[0082]

[0083] Where P[M] represents the logarithmic amplitude of the peak maximum point of the Mth harmonic. The mean of the harmonic peak-to-valley ratio of multiple harmonics is calculated as the characteristic parameter of the harmonic peak-to-valley ratio;

[0084] The peak-to-valley ratio of the first harmonic is a special case of calculating the peak-to-valley ratio of the Mth harmonic, and the steps are similar; in the calculation of the peak-to-valley ratio of the first harmonic, v p1 The calculation formula is as follows:

[0085]

[0086]

[0087]

[0088] In the calculation of the peak-to-valley ratio of a first harmonic, v p2 The calculation formula is as follows:

[0089]

[0090]

[0091]

[0092] The formula for calculating the peak-to-valley ratio (FHPVR) of the first harmonic is as follows:

[0093]

[0094] The calculation of the harmonic energy decay rate (HDR) is divided into two cases: the maximum harmonic order num is greater than or equal to four, and the logarithmic value is greater than zero. In the first case, let head be the second largest logarithmic amplitude of the first three harmonics, with harmonic order s, and let tail be the mean of the harmonics with logarithmic amplitudes greater than zero among the harmonics of orders num-2, num-1, and num. In this case, the formula for calculating the harmonic energy decay rate HDR is:

[0095]

[0096] If the logarithmic value is greater than zero and the maximum harmonic order num is less than or equal to four, then let head be the maximum logarithmic amplitude value among the first num-1 harmonics, with the corresponding harmonic order being s, and tail be the logarithmic amplitude value of the num harmonic; in this case, the formula for calculating the harmonic energy decay rate HDR is:

[0097]

[0098] In the calculation of the energy ratio of odd harmonics to even harmonics in the first 6 harmonics, let matrix HWB be HB. ij The energy ratio R of the first 6th harmonic odd harmonics to even harmonics oe6 The calculation formula is as follows:

[0099]

[0100] The row and column indices of the matrix both start from 1, and wid is HB. ij The length of a line;

[0101] When calculating the ratio of harmonic energy with a logarithmic value greater than zero to the total energy in the frequency domain, set all elements less than or equal to zero in matrix HWB to zero, resulting in a new matrix UPHB. ij The formula for calculating the ratio of the frequency domain logarithmic value of the zero harmonic energy to the total energy, UHTR, is as follows:

[0102]

[0103] In the formula, row and col are matrices UPHB. ij The number of rows and columns;

[0104] The calculation of the harmonic peak-to-valley ratio frequency spacing requires calculating the value of the harmonic trough v and the absolute value of the difference between the maximum logarithmic value P[M] of the harmonic peak for the Mth harmonic. The absolute value of the difference is compared with the frequency spacing between the harmonic trough and peak. The average of this characteristic value for multiple harmonics is the final harmonic peak-to-valley ratio frequency spacing value. Generally, the calculation is performed on the fifth to eighth harmonics in the logarithmic spectrum, where the logarithmic value of the harmonic peak is greater than zero, and the average value is taken. For a harmonic peak, there are two troughs, left and right. The value of the left trough is v. p1 The value of the trough on the right is v p2 ;v p1 and v p2 The calculation method is the same as that used for the peak-to-valley ratio of the Mth harmonic; therefore, the formula for calculating the frequency spacing DFR (peak-to-valley ratio) of the Mth harmonic is as follows:

[0105]

[0106] The DFR of multiple harmonics is calculated and the average value is taken as the harmonic peak-to-valley ratio frequency spacing value.

[0107] The peak-to-valley ratio frequency spacing of the first harmonic is a special case of the peak-to-valley ratio frequency spacing of the Mth harmonic, where the left and right valley values ​​are v. p1 and v p2 The calculation is consistent with the peak-to-valence ratio of the first harmonic; its calculation formula is:

[0108]

[0109] The energy ratio of the first and second harmonics is the ratio of the maximum logarithmic amplitude of the first harmonic to the maximum logarithmic amplitude of the second harmonic.

[0110] S402. Extract the fine time-domain features, including the standard deviation of the single-cycle fundamental frequency sequence, the mean of the single-cycle fundamental frequency sequence, the median of the single-cycle fundamental frequency sequence, the distance between the mean and the median of the single-cycle fundamental frequency sequence, and the number of times the single-cycle fundamental frequency sequence crosses the mean per unit time; let the time-domain single-cycle fundamental frequency sequence obtained in step S3 be F[n]. T ], its length is Len,n T Indicates the index of a single period in the time domain. Mean f of the single-period fundamental frequency sequence. mean The calculation formula is:

[0111]

[0112] Standard deviation f of single-cycle fundamental frequency sequence std The calculation formula is:

[0113]

[0114] For F[n T If the fundamental frequency sequence is sorted, then the median f of the single-period fundamental frequency sequence is... median The calculation formula is:

[0115]

[0116] Distance f between the mean and median of a single-period fundamental frequency sequence mm The calculation formula is:

[0117] f mm =|f mean -f median | Formula (45)

[0118] The formula for calculating the number of times Nc, the fundamental frequency sequence of a single period passes the mean per unit time, is as follows:

[0119]

[0120] Where T1 is the duration of the time-domain audio signal, and the expression for u(t) is as follows:

[0121]

[0122] S403. Extracting the static features of 13-dimensional Mel frequency cepstral coefficients: First, preprocess the input audio data. Second, perform an FFT transform on the audio and square the FFT magnitude to obtain the audio energy spectrum. Then, use a Mel filter bank to filter and obtain the energy output of each filter. Finally, perform a logarithmic operation on the filter output and obtain the Mel frequency cepstral coefficient features through Discrete Cosine Transform (DCT).

[0123] Furthermore, in step S5, the process of evaluating feature correlation and coefficient of variation, and establishing a random forest model to classify musical instruments is as follows:

[0124] S501. Calculate the Pearson correlation coefficient between each element of the feature vector and evaluate the effectiveness of the feature;

[0125] S502. Calculate the coefficient of variation of each element of the feature vector and evaluate the effectiveness of the features;

[0126] S503. The time-frequency fine features and Mel frequency cepstral coefficient features are fused into a new feature vector. The fused features are then fed into a random forest model for training. Five-fold cross-validation is used during training.

[0127] The processing procedure of the random forest model is as follows: the training set is randomly sampled multiple times with replacement using the bootstrap method, and the sampled training samples are θ.

[0128] For each training sample θ, a corresponding decision tree is formed; in each constructed decision tree, q features are randomly selected from the Q features of each training sample θ for classification, where q ≤ Q. and Indicates to Round down; calculate the Gini impurity of the selected q features, and consider the feature with the smallest Gini impurity as the optimal classification feature; based on the selected optimal classification feature, split the node containing the optimal classification feature into two parts, and select the second-best classification feature from the remaining features; repeat the above steps until the number of remaining features θ is 0 or the current decision tree achieves the best classification effect.

[0129] All 30 decision trees obtained are combined to form a random forest model, and the final classification result is obtained through voting.

[0130] Random forest models can handle input samples with high-dimensional features without requiring dimensionality reduction. The basic unit of a random forest model is a decision tree, which combines multiple weak classifiers and uses voting or averaging to achieve high accuracy and generalization performance in the overall model.

[0131] The present invention has the following advantages and effects compared with the prior art:

[0132] (1) This invention proposes multiple fine-grained frequency domain features. Compared to the large-scale description of traditional frequency domain features, it focuses on describing the fine structure of the spectrum and proposes multiple indicators to characterize the fine features of harmonic peaks, spikes, and harmonic troughs, which can more accurately characterize the spectral structure that determines timbre. Fine-grained frequency domain features can better characterize the subtle differences in the spectrum caused by the instrument's playing characteristics and the characteristics of the instrument's resonating cavity. For example, the logarithmic spectrum diagrams of trumpets and violins have spikes around the maximum harmonic peak, which are proportionally larger than the maximum harmonic peak, while the harmonic peaks in the logarithmic spectrum diagram of a piano are very sharp and do not have this phenomenon. Fine-grained frequency domain features can describe these subtle differences in the spectrum more accurately than traditional frequency domain features.

[0133] (2) Compared with traditional time-domain features, the fine time-domain feature proposed in this invention accurately marks the interval of a single fundamental frequency period and records the fundamental frequency change sequence, focusing on describing the fine changes in the fundamental frequency. For different musical instruments, the fundamental frequency stability of their single-tone audio frequencies varies. For example, the fundamental frequency stability of a trumpet, which relies on lip vibration to control its fundamental frequency, is weaker than that of a piano string that vibrates freely after being struck. The fine time-domain feature proposed in this invention can accurately characterize this difference, while traditional time-domain features do not pay attention to this characteristic.

[0134] (3) The instrument classification error rate of the random forest model with input fusion features obtained in this invention is reduced by 5.7 percentage points compared with the random forest model with static features using only 13-dimensional Mel frequency cepstral coefficients.

[0135] (4) This invention designs a cyclic iterative fundamental frequency estimation algorithm based on autocorrelation algorithm and narrowband spectrum energy frequency estimation for fundamental frequency estimation, achieving accurate fundamental frequency estimation under low signal-to-noise ratio conditions. Based on the estimated fundamental frequency, this invention designs a harmonic marking algorithm based on adaptive window search and a time-domain single-period fundamental frequency sequence estimation algorithm to accurately mark the harmonic positions of the spectrum and obtain the single-period fundamental frequency sequence. On the one hand, this overcomes the problem of harmonic offset during harmonic marking, achieving accurate harmonic marking. On the other hand, the single-period fundamental frequency sequence is based on a continuous time-domain single period, which better characterizes the continuous change of the time-domain fundamental frequency. Attached Figure Description

[0136] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this application, illustrate exemplary embodiments of the invention and, together with their description, serve to explain the invention and do not constitute an undue limitation thereof. In the drawings:

[0137] Figure 1 This is a flowchart of a musical instrument classification method based on fine time-frequency analysis in an embodiment of the present invention;

[0138] Figure 2 This is a flowchart of the iterative fundamental frequency estimation algorithm based on autocorrelation algorithm and narrowband spectrum energy frequency estimation algorithm in an embodiment of the present invention;

[0139] Figure 3 This is a flowchart of the harmonic labeling algorithm based on adaptive window search in an embodiment of the present invention;

[0140] Figure 4 This is a flowchart of the time-domain single-period fundamental frequency sequence estimation algorithm in an embodiment of the present invention. Detailed Implementation

[0141] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0142] Example 1

[0143] Figure 1 This is a flowchart of a musical instrument classification method based on fine time-frequency analysis provided in an embodiment of the present invention.

[0144] S1. Input a single-tone audio note of a musical instrument and accurately estimate the fundamental frequency of the audio.

[0145] like Figure 2 This is a flowchart of the iterative fundamental frequency estimation algorithm based on autocorrelation algorithm and narrowband spectrum energy frequency estimation algorithm in an embodiment of the present invention;

[0146] S101. The short-time energy method is used to detect the single-tone audio endpoints of musical instruments and extract the valid audio. The formula for calculating short-time energy is as follows:

[0147]

[0148] In the formula, s[n] is the input single-note audio signal of the musical instrument, w[n] is the window function, n is the time-domain index of s[n] and w[n], and N is the time-domain index of s[n] and w[n]. f It is the frame length of the signal frame, z is the shift amount of the window function, and E z It is the short-time energy corresponding to a translation amount of z; E corresponds to different translation amounts. z A short-time energy sequence is formed; a threshold is set for the short-time energy of the effective audio start and end points of the instrument's single-tone audio. When the short-time energy of the current frame exceeds the threshold of the start point, the midpoint of the current frame is recorded as the start point of the effective audio. When the short-time energy of the current frame is lower than the threshold of the end point, the midpoint of the current frame is recorded as the end point of the effective audio.

[0149] In this example, the main focus is on detecting the start point of a single-tone audio signal. The short-time energy threshold used is 0.3 times the maximum value of the short-time energy of all audio frames.

[0150] S102. The fundamental frequency f0 of the effective audio x[n] is roughly estimated by time-domain autocorrelation. The autocorrelation formula for discrete signals is as follows:

[0151]

[0152] In the formula, τ represents the lag value, i.e., the shift. After shifting x[n] to the left by τ units, it is multiplied by the original signal. Finally, the product results are summed to obtain the autocorrelation value under the current lag value. The autocorrelation value is the largest when the lag value is zero. Based on the relationship between the autocorrelation value and different lag values, the maximum point other than the lag value of zero is found. The difference between the lag value of the maximum point and the lag value of the point with zero lag value is a rough estimate of the fundamental frequency period, thereby estimating the fundamental frequency f0.

[0153] S103. Perform a Fast Fourier Transform (FFT) on the audio signal and record the frequency f corresponding to the current FFT maximum value. max ; Determine f based on time-domain autocorrelation results maxIs it an interference component? If so, remove the component from the FFT spectrum and then re-obtain the frequency f corresponding to the maximum value of the FFT. max Repeat the above process until it is believed that f max It is not an interfering component;

[0154] In this example, the number of FFT points is 4096 when the fundamental frequency f0 of the autocorrelation estimation is less than 400 Hz, and 2048 when it is greater than or equal to 400 Hz.

[0155] S104. Calculate the denoised frequency value f of the frequency corresponding to the maximum peak value of the Fast Fourier Transform using narrowband spectrum energy frequency estimation. c To eliminate noise interference, the denoised frequency value f of the FFT spectrum is calculated using narrowband spectral energy frequency estimation when the index of the FFT spectrum array starts from 1. c The formula is as follows:

[0156]

[0157] Where k0 is the FFT index corresponding to the maximum peak value of the Fast Fourier Transform, k is the index of the FFT spectrum, and |X k | is the FFT modulus corresponding to the k-th point in the FFT spectrum, F s N is the sampling rate of the audio signal, and N is the number of FFT points.

[0158] S105, Based on the noise reduction frequency value f c Based on the fundamental frequency f0 estimated by autocorrelation, the maximum peak value is calculated to be approximately located at n. m At the subharmonic level; where:

[0159]

[0160] The fundamental frequency estimate is obtained as f1:

[0161]

[0162] S106. Resample the audio time-domain signal based on the fundamental frequency estimate f1, and set the new sampling rate f. s The calculation formula is as follows:

[0163] f s =f1×2 m Formula (6)

[0164] The value of m must satisfy:

[0165]

[0166] To ensure that the new sampling rate is not too large, which would prevent the number of sampling points N from covering a sufficient fundamental frequency period, the resampled time-domain signal is executed again from step S102 until the fundamental frequency estimate f1 changes by less than 0.001 compared to the fundamental frequency estimate of the previous cycle or the number of cycles reaches the limit of 50, and the final fundamental frequency estimate f1 is output.

[0167] S2. Perform FFT frequency domain harmonic marking based on the fundamental frequency of the audio.

[0168] Figure 3 This is the flowchart of the harmonic labeling algorithm based on adaptive window search in this embodiment of the invention;

[0169] S201. The short-time energy method is used to detect the audio endpoints of musical instruments and extract the effective audio portion. In this example, the main focus is on detecting the starting point of a single-tone audio signal. The short-time energy threshold used is 0.3 times the maximum value of the short-time energy of all audio frames.

[0170] S202. Based on the fundamental frequency estimate f1 from step S1, the audio time-domain signal is resampled to make the audio sampling rate a power of f1, thereby reducing spectral leakage, and the FFT spectrum is converted into a logarithmic spectrum. In this embodiment, due to the need for refined spectral analysis, the number of FFT points during harmonic labeling is 16384.

[0171] S203. Determine the left endpoint l and right endpoint r of the next harmonic peak search range on the FFT spectrum. If it is the first 3 harmonics, the calculation formula is as follows:

[0172]

[0173]

[0174] Where P is the position index of the previous harmonic peak in the FFT spectrum. If the position of the first harmonic peak is marked, then P = 1.

[0175] If the 4th harmonic and subsequent harmonics are being marked, the formulas for calculating the left endpoint l and the right endpoint r are as follows:

[0176] l=P+bT-2 Formula (10)

[0177] r = P + 1.09T1 Formula (11)

[0178] Where bT is half the distance between the third harmonic peak and the first harmonic peak on the FFT spectrum; T1 is the FFT point interval between the two harmonics closest to the currently marked harmonic on the FFT spectrum;

[0179] S204. Find all maximum points within the search range, and take the maximum value of all maximum points as the desired peak point of the current harmonic. If there are no maximum points within the search range, expand the left and right boundaries of the search range when marking the first four harmonics, and only expand the right boundary when marking other harmonics. When expanding the search boundary of the first two harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions:

[0180]

[0181]

[0182] When expanding the search boundaries for the 3rd and 4th harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions:

[0183] lP≥0.7T1 formula (14)

[0184] rP≤1.3T1 formula (15)

[0185] When expanding the search boundary for the 5th and higher harmonics, the right endpoint r must satisfy the following condition:

[0186] rP≤1.3T1 formula (16)

[0187] Then, within the new search range, the maximum value of all maxima is found as the peak value of the current harmonic. All FFT logarithmic amplitude values ​​and the index value of the maximum value within 0.1 times the harmonic interval centered on the maximum value are recorded and stored in matrix HWB and sequence HWI, respectively. Each row of HWB stores the peak-to-peak value of each harmonic and the spectral logarithmic value of its vicinity. If a maximum value is still not found, the search range is expanded and the process is repeated.

[0188] S3. The process of estimating the time-domain single-cycle fundamental frequency sequence based on the fundamental frequency of the audio is as follows:

[0189] Figure 4 This is a flowchart of the time-domain single-period fundamental frequency sequence estimation algorithm in an embodiment of the present invention;

[0190] S301. The short-time energy method is used to detect the audio endpoints of musical instruments and extract the effective audio portion. In this example, the main focus is on detecting the starting point of a single-tone audio signal. The short-time energy threshold used is 0.3 times the maximum value of the short-time energy of all audio frames.

[0191] S302. Use a high-pass filter to filter out low-frequency interference. In this example, the cutoff frequency is set to 0.3 times the fundamental frequency estimate f1 from step S1.

[0192] S303. Use cross-correlation to mark the boundary point of the single fundamental frequency period. Taking a point o in the middle of the audio signal as the origin, extend equidistantly to both sides. The extension distance is the length of the single fundamental frequency period obtained in the previous calculation. The signal within the extension range is a shift vector. The length of the shift vector is twice the length of the single fundamental frequency period obtained in the previous calculation, thus having the characteristic of adaptive adjustment. If it is the first cross-correlation calculation, the single frequency period length refers to the fundamental frequency estimation result in step S2. Taking point o as the origin, take four signals of the single fundamental frequency period length to the left or right as a fixed direction. The quantity, where the single fundamental frequency period length is also the single fundamental frequency period length obtained in the previous calculation; the shift vector is used as a sliding window to perform cross-correlation with the fixed vector to obtain the correspondence between the cross-correlation value and the relative displacement; the extreme point sequence of the cross-correlation value and the corresponding relative displacement sequence are recorded from the correspondence; reasonable values ​​are selected from the relative displacement sequence, and the next single fundamental frequency period dividing point is calculated with point o as the reference. The distance between the dividing point and o is the single fundamental frequency period length; this process is repeated to obtain the time-domain single-cycle fundamental frequency sequence. The calculation formula for the cross-correlation r[τ] is as follows:

[0193]

[0194] In the formula, τ represents the lag value, i.e. the shift amount. g[n] and v[n] are two different signals. After shifting v[n] to the left by a distance of τ units, it is multiplied by g[n]. The sum of the product results is the similarity under the current lag value.

[0195] S4. The process of extracting time-frequency fine features and Mel-frequency cepstral coefficient features is as follows:

[0196] S401. Extract fine frequency domain features, including: centroid of harmonic energy with logarithmic values ​​greater than zero, average attenuation ratio of harmonic peaks and spikes, maximum harmonic order with logarithmic values ​​greater than zero, slope-bandwidth ratio between two harmonics, harmonic peak-to-valley ratio, harmonic energy decrease rate, energy ratio of odd to even harmonics in the first 6 harmonics, peak-to-valley ratio of the first harmonic, decrease rate of slope-to-bandwidth ratio between harmonics, ratio of harmonic energy with logarithmic values ​​greater than zero to total energy, harmonic peak-to-valley ratio to half-harmonic spacing, and peak-to-valley ratio to half-harmonic spacing of the first harmonic.

[0197] From the HWB obtained in step S2, the logarithmic amplitude sequence P[n] of the maximum value points of each harmonic peak can be obtained. h The numbers less than or equal to zero in the logarithmic amplitude array are set to zero, with indices starting from 1. The formula for calculating the centroid c of harmonic energy with a logarithmic value greater than zero is as follows:

[0198]

[0199] Where L P It is P[n h The length of ] is also the maximum harmonic order, nh It is the current harmonic order, and also P[n h The index of ];

[0200] The maximum harmonic order with a logarithmic value greater than zero is the logarithmic amplitude array P[n]. h The largest index of the non-zero value corresponding to the less than or equal to zero numbers in the array after setting them to zero, with the index starting from 1;

[0201] The calculation of the average attenuation ratio of harmonic peak glitches requires that, in step S2, the logarithmic spectrum of the FFT spectrum is greater than zero, and that the maxima A[n] of the twenty non-harmonic peaks closest to the maximum harmonic peak within a harmonic spacing range centered on the maximum harmonic peak value be obtained. p ], n p It is A[n p If the number of maxima within the range is less than twenty, then all maxima are taken; if the number of maxima is less than or equal to 10, then A[n] is taken directly. p The mean of A[n] is used as the average amplitude of the spikes around the harmonic peak; if the number of maxima is greater than 10, then in A[n] p Find a continuous set of maxima of length 10 in the FFT index, where the FFT index span of this set of maxima either contains the index of the maximum harmonic peak or there are no other maxima between it and the index of the maximum harmonic peak, and the sum is the largest. Calculate the mean of this set of maxima as the average amplitude of the spikes around the harmonic peak. Let the average amplitude of the spikes around the harmonic peak be 'a'. The formula for the harmonic peak attenuation ratio (PDR) of a single harmonic is as follows:

[0202]

[0203] Where p is the maximum amplitude of the current harmonic peak, and the average of the harmonic peak attenuation ratios of multiple harmonics is the average harmonic peak attenuation ratio.

[0204] Let the logarithmic spectrum of the FFT spectrum in step S2 be mag[k], and the HWI obtained in step S2 be the index sequence of the maximum harmonic peaks in mag[k], abbreviated as hi[n]. h If all the indices of the sequence start from 1, then the formula for calculating the slope bandwidth ratio (SFWR) between the 3rd and 4th harmonics is:

[0205]

[0206] When calculating the interharmonic slope bandwidth ratio descent rate, first calculate the interharmonic slope bandwidth ratio between the DC component and the first harmonic in the logarithmic spectrum as r1. Then calculate the interharmonic slope bandwidth ratio for each harmonic interval between the 5th and 8th harmonics and take the average as r2. If there are less than 8 harmonics, then take the highest harmonic. The formula for calculating the interharmonic slope bandwidth ratio descent rate SFWRD is as follows:

[0207]

[0208] The harmonic peak-to-valley ratio is calculated based on the portion of the logarithmic spectrum from which the logarithm is greater than zero using the FFT. All elements less than or equal to zero in mag[k] are set to zero, resulting in a new sequence MAG[k]. For a harmonic peak, there are two troughs, left and right, with the left trough value being v. p1 Its value is the mean of MAG[k] within the trough range of the left harmonic spacing; the trough value on the right is v. p2 Its value is the average of MAG[k] within the trough range of the right harmonic spacing; in the calculation of the peak-to-trough ratio of the Mth harmonic, v p1 The calculation formula is as follows:

[0209]

[0210] Where L is the left boundary for calculating the trough value, and R is the right boundary for calculating the trough value;

[0211]

[0212]

[0213] The formula for calculating v2 in the peak-to-valley ratio of the Mth harmonic is as follows:

[0214]

[0215]

[0216]

[0217] Given that the logarithmic amplitude sequence of the peak maximum points of each harmonic is P[n] h The formula for calculating the peak-to-valley ratio (MR) of the Mth harmonic is as follows:

[0218]

[0219] Where P[M] represents the logarithmic amplitude of the peak maximum point of the Mth harmonic. The mean of the harmonic peak-to-valley ratio of multiple harmonics is calculated as the characteristic parameter of the harmonic peak-to-valley ratio;

[0220] The peak-to-valley ratio of the first harmonic is a special case of calculating the peak-to-valley ratio of the Mth harmonic, and the steps are similar; in the calculation of the peak-to-valley ratio of the first harmonic, v p1 The calculation formula is as follows:

[0221]

[0222]

[0223]

[0224] In the calculation of the peak-to-valley ratio of a first harmonic, v p2 The calculation formula is as follows:

[0225]

[0226]

[0227]

[0228] The formula for calculating the peak-to-valley ratio (FHPVR) of the first harmonic is as follows:

[0229]

[0230] The calculation of the harmonic energy decay rate (HDR) is divided into two cases: the maximum harmonic order num is greater than or equal to four, and the logarithmic value is greater than zero. In the first case, let head be the second largest logarithmic amplitude of the first three harmonics, with harmonic order s, and let tail be the mean of the harmonics with logarithmic amplitudes greater than zero among the harmonics of orders num-2, num-1, and num. In this case, the formula for calculating the harmonic energy decay rate HDR is:

[0231]

[0232] If the logarithmic value is greater than zero and the maximum harmonic order num is less than or equal to four, then let head be the maximum logarithmic amplitude value among the first num-1 harmonics, with the corresponding harmonic order being s, and tail be the logarithmic amplitude value of the num harmonic; in this case, the formula for calculating the harmonic energy decay rate HDR is:

[0233]

[0234] In the calculation of the energy ratio of odd harmonics to even harmonics in the first 6 harmonics, let matrix HWB be HB. ij The energy ratio R of the first 6th harmonic odd harmonics to even harmonics oe6 The calculation formula is as follows:

[0235]

[0236] The row and column indices of the matrix both start from 1, and wid is HB. ij The length of a line;

[0237] When calculating the ratio of harmonic energy with a logarithmic value greater than zero to the total energy in the frequency domain, set all elements less than or equal to zero in matrix HWB to zero, resulting in a new matrix UPHB. ij The formula for calculating the ratio of the frequency domain logarithmic value of the zero harmonic energy to the total energy, UHTR, is as follows:

[0238]

[0239] In the formula, row and col are matrices UPHB. ij The number of rows and columns;

[0240] The calculation of the harmonic peak-to-valley ratio frequency spacing requires calculating the value of the harmonic valley v for the Mth harmonic and the absolute value of the difference between the maximum logarithmic value P[M] of the harmonic peak. The absolute value of the difference is compared with the frequency spacing between the harmonic valley and peak to obtain the Mth harmonic peak-to-valley ratio frequency spacing. The final harmonic peak-to-valley ratio frequency spacing value is obtained by averaging this characteristic value for multiple harmonics. Generally, the calculation is performed on the harmonics with a logarithmic peak value greater than zero in the 5th to 8th harmonics in the logarithmic spectrum, and the average value is taken. For a harmonic peak, there are two valleys, left and right. The value of the left valley is v. p1 The value of the trough on the right is v p2 ;v p1 and v p2 The calculation method is the same as that used for the peak-to-valley ratio of the Mth harmonic; therefore, the formula for calculating the frequency spacing DFR (peak-to-valley ratio) of the Mth harmonic is as follows:

[0241]

[0242] The DFR of multiple harmonics is calculated and the average value is taken as the harmonic peak-to-valley ratio frequency spacing value.

[0243] The peak-to-valley ratio frequency spacing of the first harmonic is a special case of the peak-to-valley ratio frequency spacing of the Mth harmonic, where the left and right valley values ​​are v. p1 and v p2 The calculation is consistent with the peak-to-valence ratio of the first harmonic; its calculation formula is:

[0244]

[0245] The energy ratio of the first and second harmonics is the ratio of the maximum logarithmic amplitude of the first harmonic to the maximum logarithmic amplitude of the second harmonic.

[0246] S402. Extract the fine time-domain features, including the standard deviation of the single-cycle fundamental frequency sequence, the mean of the single-cycle fundamental frequency sequence, the median of the single-cycle fundamental frequency sequence, the distance between the mean and the median of the single-cycle fundamental frequency sequence, and the number of times the single-cycle fundamental frequency sequence crosses the mean per unit time; let the time-domain single-cycle fundamental frequency sequence obtained in step S3 be F[n]. T ], its length is Len,n T Indicates the index of a single period in the time domain. Mean f of the single-period fundamental frequency sequence. mean The calculation formula is:

[0247]

[0248] Standard deviation f of single-cycle fundamental frequency sequence std The calculation formula is:

[0249]

[0250] For F[n T If the fundamental frequency sequence is sorted, then the median f of the single-period fundamental frequency sequence is... median The calculation formula is:

[0251]

[0252] Distance f between the mean and median of a single-period fundamental frequency sequence mm The calculation formula is:

[0253] f mm =|f mean -f median | Formula (45)

[0254] The number of times N of a single-period fundamental frequency sequence passes the mean per unit time c The calculation formula is:

[0255]

[0256] Where T1 is the duration of the time-domain audio signal, and the expression for u(t) is as follows:

[0257]

[0258] S403. Extracting the static features of 13-dimensional Mel frequency cepstral coefficients: First, preprocess the input audio data. Second, perform an FFT transform on the audio and square the FFT magnitude to obtain the audio energy spectrum. Then, use a Mel filter bank to filter and obtain the energy output of each filter. Finally, perform a logarithmic operation on the filter output and obtain the Mel frequency cepstral coefficient features through Discrete Cosine Transform (DCT).

[0259] S5. Based on the extracted video fine features and Mel frequency cepstral coefficient features, feature correlation and coefficient of variation are evaluated, and a random forest model is established to classify musical instruments as follows:

[0260] S501. Calculate the Pearson correlation coefficient between each element of the feature vector and evaluate the effectiveness of the feature;

[0261] S502. Calculate the coefficient of variation of each element of the feature vector and evaluate the effectiveness of the features;

[0262] S503. The time-frequency fine features and Mel frequency cepstral coefficient features are fused into a new feature vector. In this embodiment, the fused features of 822 audio samples from 5 musical instruments are input into the random forest model for training, and five-fold cross-validation is used during training. Then, the 13-dimensional Mel frequency cepstral coefficient static features are input into the random forest model separately for training, and five-fold cross-validation is used during training. Random forests can handle input samples with high-dimensional features without dimensionality reduction. The basic unit of random forest is the decision tree, which combines multiple weak classifiers, and the final result is obtained by voting or averaging, so that the overall model has high accuracy and generalization performance. The musical instrument classification error rate of the random forest model with fused input features obtained in this example is 6% lower than that of the random forest model using only the 13-dimensional Mel frequency cepstral coefficient static features, as shown in Table 1.

[0263] Table 1. Feature Types and Classification Accuracy

[0264]

[0265]

[0266] Example 2

[0267] Based on the musical instrument classification method based on time-frequency fine analysis disclosed in Embodiment 1 above, this embodiment further provides the classification implementation process of the musical instrument classification method based on time-frequency fine analysis, as follows:

[0268] S1. Referring to the corresponding steps in Example 1, obtain the fundamental frequency f1 of the input violin single-note audio.

[0269] S2. Record all FFT logarithmic amplitude values ​​and index values ​​of the maximum values ​​within 0.1 times the harmonic interval centered on the maximum value of each harmonic, and store them in matrix HWB and sequence HWI respectively. Each row of HWB stores the peak logarithmic value of each harmonic and the spectral logarithmic value of its vicinity.

[0270] S3. Refer to the corresponding steps in Example 1, which will not be repeated here;

[0271] S4. Extract the time-frequency fine features and the 13-dimensional Mel frequency cepstral coefficient static features. Use MATLAB code to automatically calculate each feature. The calculated feature values ​​of the audio are shown in the table below.

[0272] Table 2. Time-frequency fine features and static features of 13-dimensional Vimel frequency cepstral coefficients of the corresponding audio.

[0273]

[0274]

[0275] S5. Based on the input fusion features obtained in Example 1, the random forest model, with the above features input into the model, yields a classification result of "violin", which is correct.

[0276] The above embodiments are preferred embodiments of the present invention, but the embodiments of the present invention are not limited to the above embodiments. Any changes, modifications, substitutions, combinations, or simplifications made without departing from the spirit and principle of the present invention shall be considered equivalent substitutions and shall be included within the protection scope of the present invention.

Claims

1. A musical instrument classification method based on fine time-frequency analysis, characterized in that, The musical instrument classification method includes the following steps: Step S1: Input the single-note audio of the instrument and accurately estimate the fundamental frequency of the audio. Step S2: Perform FFT frequency domain harmonic labeling based on the fundamental frequency of the audio; the logarithmic spectrum of the FFT spectrum obtained in step S2 is mag[k], and the index sequence of the maximum harmonic peaks in mag[k] is simply called hi[n]. h ], n h It is the current harmonic order; Step S3: Estimate the time-domain single-cycle fundamental frequency sequence based on the fundamental frequency of the audio. Step S4: Extract time-frequency fine features and Mel-frequency cepstral coefficient features, wherein, The time-frequency fine features include frequency domain fine features and time domain fine features. The frequency domain fine features include: centroid of harmonic energy with logarithmic values ​​greater than zero, average attenuation ratio of harmonic peaks and spikes, maximum harmonic order with logarithmic values ​​greater than zero, slope-bandwidth ratio between two harmonics, harmonic peak-to-valley ratio, harmonic energy decay rate, energy ratio of odd to even harmonics in the first 6 harmonics, peak-to-valley ratio of the first harmonic, decay rate of the slope-to-bandwidth ratio between harmonics, ratio of harmonic energy with logarithmic values ​​greater than zero to total energy, harmonic peak-to-valley ratio frequency spacing, and first harmonic peak-to-valley ratio frequency spacing. The time-domain fine features include: standard deviation of a single-cycle fundamental frequency sequence, mean of a single-cycle fundamental frequency sequence, median of a single-cycle fundamental frequency sequence, distance between the mean and median of a single-cycle fundamental frequency sequence, and number of times the single-cycle fundamental frequency sequence crosses the mean per unit time. The calculation of the average attenuation ratio of the harmonic peak glitch requires that, in step S2, the logarithmic spectrum of the FFT spectrum is greater than zero, and that the 20 non-harmonic peak maxima closest to the harmonic peak maximum value are obtained within a harmonic spacing range centered on the harmonic peak maximum value. , yes If the number of maxima within the range is less than 20, then all maxima are taken; if the number of maxima is less than or equal to 10, then A is taken directly. The mean of A[] is used as the average amplitude of the spikes around the harmonic peak; if the number of maxima is greater than 10, then in A[ Find a continuous set of maxima of length 10 in the FFT index, where the FFT index span of this set of maxima either contains the index of the maximum harmonic peak or there are no other maxima between it and the index of the maximum harmonic peak, and the sum is the largest. Calculate the mean of this set of maxima as the average amplitude of the spikes around the harmonic peak. Let the average amplitude of the spikes around the harmonic peak be 'a'. The formula for the harmonic peak attenuation ratio (PDR) of a single harmonic is as follows: , formula (19), Where p is the maximum amplitude of the current harmonic peak, and the average of the harmonic peak attenuation ratios of multiple harmonics is the average harmonic peak attenuation ratio. The calculation of the harmonic peak-to-valley ratio frequency spacing requires calculating the value of the harmonic valley for the Mth harmonic. The absolute value of the difference between the maximum logarithmic values ​​of the harmonic peaks P[M] and the frequency spacing between the harmonic troughs and peaks is the Mth harmonic peak-trough difference frequency spacing. The final harmonic peak-trough difference frequency spacing value is obtained by averaging this characteristic value for multiple harmonics. The average value is calculated for the fifth to eighth harmonics in the logarithmic spectrum where the logarithmic values ​​of the harmonic peaks are greater than zero. For a harmonic peak, there are two troughs, left and right. The value of the left trough is... The trough value on the right is ; and The calculation method is consistent with the algorithm used in the peak-to-valley ratio of the Mth harmonic; therefore, the peak-to-valley difference ratio of the Mth harmonic is the frequency spacing. The calculation formula is as follows: Formula (40), where Fs is the sampling rate of the audio signal and N is the number of FFT points; Calculation of multiple harmonics The average value is the harmonic peak-to-valley ratio frequency spacing value; The peak-to-valley ratio frequency spacing of the first harmonic is a special case of the peak-to-valley ratio frequency spacing of the Mth harmonic, with the left and right valley values... and The calculation is consistent with the peak-to-valley ratio of the first harmonic, and the calculation formula is: , formula (41), The energy ratio of the first and second harmonics is the ratio of the maximum logarithmic amplitude of the first harmonic to the maximum logarithmic amplitude of the second harmonic. The average of the peak-to-valley ratios of multiple harmonics is used as the harmonic peak-to-valley ratio. Step S5: Based on the extracted time-frequency fine features and Mel frequency cepstral coefficient features, perform feature correlation and coefficient of variation evaluation, and establish a random forest model to classify musical instruments.

2. The musical instrument classification method based on time-frequency fine analysis according to claim 1, characterized in that, The process of step S1 is as follows: S101. The short-time energy method is used to detect the single-tone audio endpoints of musical instruments and extract the valid audio. The formula for calculating short-time energy is as follows: , formula (1), In the formula It is the input single-note audio signal of the musical instrument. It is a window function. yes and Time-domain index, It is the frame length of the signal frame. It is the translation amount of the window function. The translation amount is The corresponding short-time energy; the corresponding short-time energy for different translation amounts. A short-time energy sequence is formed; a threshold is set for the short-time energy of the effective audio start and end points of the instrument's single-tone audio. When the short-time energy of the current frame exceeds the threshold of the start point, the midpoint of the current frame is recorded as the start point of the effective audio. When the short-time energy of the current frame is lower than the threshold of the end point, the midpoint of the current frame is recorded as the end point of the effective audio. S102. Roughly estimate the effective audio through time-domain autocorrelation. base frequency The formula for the autocorrelation of a discrete signal is as follows: , Official (2), In the formula This represents the lag value, i.e., the shift amount, which shifts x[n] to the left. The distance of a given unit is multiplied by the original signal, and the product is summed to obtain the autocorrelation value at the current lag value. The autocorrelation value is maximum when the lag value is zero. Based on the relationship between the autocorrelation value and different lag values, the maximum locality (MLT) is found except when the lag value is zero. The difference between the MLT and the point where the lag value is zero is a rough estimate of the fundamental frequency period, thus estimating the fundamental frequency. ; S103. Perform a Fast Fourier Transform (FFT) on the audio signal and record the frequency corresponding to the current FFT maximum value. ; Based on the time-domain autocorrelation results Is it an interference component? If so, remove the component from the FFT spectrum and then re-obtain the frequency corresponding to the maximum value of the FFT. Repeat the above process until it is concluded, based on the time-domain autocorrelation results, that... It is not an interfering component; S104. Calculate the denoised frequency value corresponding to the maximum peak value of the Fast Fourier Transform using narrowband spectrum energy frequency estimation. To eliminate noise interference, the denoised frequency values ​​of the FFT spectrum are calculated using narrowband spectral energy-frequency estimation when the index of the FFT spectrum array starts from 1. The formula is as follows: , formula (3), in, It is the FFT index corresponding to the maximum peak value of the Fast Fourier Transform. It is an index of the FFT spectrum. It is the first in the FFT spectrum The FFT modulus corresponding to each point; S105, Based on the noise reduction frequency value The fundamental frequency estimated by autocorrelation The maximum peak value was calculated to be approximately located at At the subharmonic level; where: , formula (4), The fundamental frequency estimate is obtained as follows: : , formula (5); S106, Based on the estimated fundamental frequency value Resampling of the audio time-domain signal at a new sampling rate The calculation formula is as follows: , Official (6), The value of m must satisfy: , formula (7), The resampled time-domain signal is executed again starting from step S102 until the fundamental frequency estimate is reached. If the change in fundamental frequency estimate compared to the previous iteration is less than 0.001, or if the iteration count reaches the limit of 50, output the final fundamental frequency estimate. .

3. The musical instrument classification method based on fine time-frequency analysis according to claim 2, characterized in that, In step S2, the FFT frequency domain harmonic marking process based on the fundamental frequency of the audio is as follows: S201. Use the short-time energy method to detect the audio endpoints of musical instruments and extract the valid audio portion; S202, Based on the fundamental frequency estimate from step S1 Resampling the audio time-domain signal increases the audio sampling rate. The power of , and convert the FFT spectrum to a logarithmic spectrum; S203. Determine the left endpoint l and right endpoint r of the next harmonic peak search range on the FFT spectrum. If it is the first 3 harmonics, the calculation formula is as follows: , formula (8), , Official (9), in This is the index of the previous harmonic peak's position in the FFT spectrum. If the position of the first harmonic peak is marked, then... ; If the 4th harmonic and subsequent harmonics are being marked, the formulas for calculating the left endpoint l and the right endpoint r are as follows: Official (10), Official (11), in, It is half the distance between the third harmonic peak and the first harmonic peak on the FFT spectrum; It is the interval between the two harmonic FFT points closest to the currently labeled harmonic in the FFT spectrum; S204. Find all maximum points within the search range, and take the maximum value of all maximum points as the desired peak point of the current harmonic. If there are no maximum points within the search range, expand the left and right boundaries of the search range when marking the first four harmonics, and only expand the right boundary when marking other harmonics. When expanding the search boundary of the first two harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions: Official (12), Official (13), When expanding the search boundaries for the 3rd and 4th harmonics, the left endpoint l and the right endpoint r must satisfy the following conditions: Official (14), Official (15), When expanding the search boundary for the 5th and higher harmonics, the right endpoint r must satisfy the following condition: , Official (16) Then, within the new search range, the maximum value of all maxima is calculated as the peak value of the current harmonic. All FFT logarithmic amplitude values ​​and the index value of the maximum value within 0.1 times the harmonic interval centered on the maximum value are recorded and stored in matrix HWB and sequence HWI, respectively. Each row of HWB stores the peak-to-peak value of each harmonic and the spectral logarithmic value of its vicinity. If the maximum value is still not found, the search range is expanded and the process is repeated until the maximum value is found.

4. The musical instrument classification method based on fine time-frequency analysis according to claim 3, characterized in that, In step S3, the process of estimating the time-domain single-cycle fundamental frequency sequence based on the fundamental frequency of the audio is as follows: S301. Use the short-time energy method to detect the audio endpoints of musical instruments and extract the valid audio portion; S302. Use a high-pass filter to filter out low-frequency interference; S303. Use cross-correlation to mark the boundary points of the single fundamental frequency period. Taking a point o in the middle of the audio signal as the origin, extend equidistantly to both sides. The extension distance is the length of the single fundamental frequency period obtained in the previous calculation. The signal within the extension range is a shift vector. The length of the shift vector is twice the length of the single fundamental frequency period obtained in the previous calculation, thus having the characteristic of adaptive adjustment. If it is the first cross-correlation operation, the single-cycle length is taken from the fundamental frequency estimation result in step S2. Taking point o as the origin, take four signals of single fundamental frequency period length to the left or right as fixed vectors. The single fundamental frequency period length is also the length of the single fundamental frequency period obtained in the previous calculation. Perform cross-correlation operation on the shift vector as a sliding window and the fixed vector to obtain the correspondence between the cross-correlation value and the relative displacement. Record the extreme point sequence of the cross-correlation value and the corresponding relative displacement sequence from the correspondence. Select reasonable values ​​from the relative displacement sequence, and calculate the next single fundamental frequency period dividing point with point o as the reference. The distance between the dividing point and o is the length of the single fundamental frequency period. Repeat this process to obtain the time-domain single-cycle fundamental frequency sequence and cross-correlation. The calculation formula is as follows: , formula (17), In the formula This represents the lag value, i.e., the translation amount. v[n] and v[n] are two different signals. Shift v[n] to the left. After a distance of one unit, and Multiply the products and sum them up to obtain the similarity at the current lag value.

5. The musical instrument classification method based on time-frequency fine analysis according to claim 4, characterized in that, The process of step S4 is as follows: S401. Extract fine features in the frequency domain, as follows: Obtain the logarithmic amplitude sequence of the maximum value points of each harmonic peak from the HWB obtained in step S2. The numbers less than or equal to zero in the logarithmic amplitude array are set to zero, with the index starting from 1. The formula for calculating the centroid c of harmonic energy with a logarithmic value greater than zero is as follows: , formula (18), in It is P[ The length of ] is also the maximum harmonic order. It is the current harmonic order, and also P[ The index of ]; The maximum harmonic order with a logarithmic value greater than zero is the logarithmic amplitude array P[ The largest index of the non-zero value corresponding to the less than or equal to zero numbers in the array after setting them to zero, with the index starting from 1; Let the logarithmic spectrum of the FFT spectrum obtained in step S2 be... The HWI obtained in step S2 is The index sequence containing the maximum value of the middle harmonic peak is simply referred to as If all the indices of the sequence start from 1, then the formula for calculating the slope bandwidth ratio (SFWR) between the 3rd and 4th harmonics is: , formula (20), When calculating the rate of decrease in the slope-bandwidth ratio between harmonics, first calculate the slope-bandwidth ratio between the DC component and the first harmonic in the logarithmic spectrum. Then calculate the inter-harmonic slope bandwidth ratio of each harmonic spacing between the 5th and 8th harmonics and take the average value. If there are fewer than 8 harmonics, then the highest harmonic is used; the formula for calculating the harmonic slope bandwidth ratio descent rate (SFWRD) is as follows: Official (21), The harmonic peak-to-valley ratio is calculated based on the portion of the logarithmic spectrum from which the logarithm is greater than zero, let If all elements less than or equal to zero are zero, a new sequence is obtained. For a harmonic peak, there are two troughs, one on the left and one on the right. The value of the left trough is... Its value is The mean value within the trough range of the left harmonic spacing; the trough value on the right is... Its value is The mean value within the trough range of the right harmonic spacing; In the calculation of the peak-to-valley ratio of the Mth harmonic The calculation formula is as follows: Official (22), in It is the left boundary for calculating the trough value. It is the right boundary for calculating the trough value; Official (23), Official (24), In the calculation of the peak-to-valley ratio of the Mth harmonic, v p2 The calculation formula is as follows: Official (25), Official (26), Official (27), Given that the logarithmic amplitude sequence of the peak maximum points of each harmonic is P[ The formula for calculating the peak-to-valley ratio (MR) of the Mth harmonic is as follows: , formula (28), in This represents the logarithmic amplitude of the peak maximum point of the Mth harmonic; In the calculation of the peak-to-valley ratio of the first harmonic The calculation formula is as follows: Official (29), Official (30), Official (31), In the calculation of the peak-to-valley ratio of a first harmonic wave The calculation formula is as follows: Official (32), Official (33), Official (34), The formula for calculating the peak-to-valley ratio (FHPVR) of the first harmonic is as follows: Official (35), The calculation of the harmonic energy decay rate is divided into the maximum harmonic order with a logarithmic value greater than zero. There are two cases: greater than four and less than or equal to four. In the first case, let... It is the second largest logarithmic amplitude of the first three harmonics, and its harmonic order is ,make for , , The mean of the harmonics with a logarithmic amplitude greater than zero in the second harmonic series. In this case, the formula for calculating the harmonic energy decay rate (HDR) is: , formula (36), If the logarithm is greater than zero, the maximum harmonic number If the number is less than or equal to four, then let For the front The maximum logarithmic amplitude value in the second harmonic, and the corresponding harmonic order is , for The logarithmic amplitude value of the second harmonic; in this case, the formula for calculating the harmonic energy decay rate (HDR) is: , formula (37), In the calculation of the energy ratio of odd harmonics to even harmonics in the first 6 harmonics, let matrix HWB be... The energy ratio of odd harmonics to even harmonics in the first 6 harmonics The calculation formula is as follows: , formula (38), The row and column indices of the matrix both start from 1. yes The length of a line; When calculating the ratio of harmonic energy with a logarithmic value greater than zero to the total energy in the frequency domain, set all elements less than or equal to zero in matrix HWB to zero to obtain a new matrix. The formula for calculating the ratio of the frequency domain logarithmic value of the zero harmonic energy to the total energy, UHTR, is as follows: Official (39), In the formula and They are matrices The number of rows and columns; S402. Extract fine features in the temporal domain, as follows: Let the time-domain single-period fundamental frequency sequence obtained in step S3 be F[n] T ], with length Len,n T Indicates the sequence number of a single period in the time domain; mean of the single-period fundamental frequency sequence. The calculation formula is: , formula (42), Standard deviation of single-period fundamental frequency sequence The calculation formula is: Official (43), For F[n T If sorted, the median of the single-period fundamental frequency sequence is... The calculation formula is: Official (44), Distance between mean and median of a single-period fundamental frequency sequence The calculation formula is: Official (45), Number of times a single-period fundamental frequency sequence passes the mean per unit time The calculation formula is: Official (46), in It is the duration of the time-domain audio signal. The expression is as follows: Official (47), S403. Extracting the static features of 13-dimensional Mel frequency cepstral coefficients: First, preprocess the input audio data. Second, perform an FFT transform on the audio and square the FFT magnitude to obtain the audio energy spectrum. Then, use a Mel filter bank to filter and obtain the energy output of each filter. Finally, perform a logarithmic operation on the filter output and obtain the Mel frequency cepstral coefficient features through Discrete Cosine Transform (DCT).

6. The musical instrument classification method based on time-frequency fine analysis according to claim 1, characterized in that, In step S5, the process of evaluating feature correlation and coefficient of variation, and establishing a random forest model to classify musical instruments is as follows: S501. Calculate the Pearson correlation coefficient between each element of the feature vector and evaluate the effectiveness of the feature; S502. Calculate the coefficient of variation of each element of the feature vector and evaluate the effectiveness of the features; S503. The time-frequency fine features and Mel frequency cepstral coefficient features are fused into a new feature vector. The fused features are then fed into a random forest model for training. Five-fold cross-validation is used during training. The random forest model's processing procedure is as follows: The training set is randomly sampled multiple times with replacement using a bootstrap method; the resulting training samples are then... ; For each training sample This forms the corresponding decision tree; in each constructed decision tree, from each training sample... of From features, q features are randomly selected for classification, where ... , ,and Indicates to " "Round down operation; calculate the selected..." The Gini impurity of each feature is calculated, and the feature with the smallest Gini impurity is considered the optimal classification feature. Based on the selected optimal classification feature, the node containing the optimal classification feature is bisected, and the second-best classification feature is selected from the remaining features. This process is repeated until... The number of remaining features is 0, or the current decision tree achieves the best classification performance; All 30 decision trees obtained are combined to form a random forest model, and the final classification result is obtained through voting.

Images