A method of harmonic structure analysis

Abstract

The application relates to the field of bionic hearing signal processing, in particular to a harmonic structure analysis method. The method comprises the following steps: obtaining multiple frames of short-time signals by performing frame division, windowing and pre-emphasis preprocessing on an original time-domain acoustic signal, then inputting the multiple frames of short-time signals into a hearing filter bank for frequency band decomposition and performing phase shift compensation by using a continuously adjustable fractional order parameter, outputting phase-aligned multichannel sub-band signals, then sequentially performing half-wave rectification, nonlinear amplitude compression and low-pass filtering on the signals of each channel to extract an electro-physiological envelope, locking a fundamental frequency by using an autocorrelation function and extracting multiple initial harmonic peak values, finally constructing a two-dimensional space-time random matrix and combining Karhunen-Loeve expansion to extract orthogonal basis functions, and reconstructing concurrent harmonic time-frequency trajectories. The application improves time distortion caused by filter group delay, effectively suppresses dynamic floor noise under non-stationary environment, and improves noise stability and anti-interference ability when analyzing harmonic structure.

Description

Technical Field

[0001] This invention relates to the field of bionic auditory signal processing, and more particularly to a method for harmonic structure analysis. Background Technology

[0002] Harmonic structure analysis is a core component of acoustic signal processing and target state recognition. Existing bionic auditory technologies typically utilize auditory filter banks to simulate the frequency response characteristics of the basilar membrane, decomposing broadband time-domain acoustic signals into multi-channel sub-band sequences. To further analyze complex sound source characteristics, the industry often extracts features from each channel and constructs a time-series matrix, thereby reconstructing the harmonic evolution trajectory of the target signal in the spatiotemporal domain to meet the urgent need for high-resolution audio analysis in modern industrial applications.

[0003] However, traditional auditory filters are prone to introducing nonlinear group delays when decomposing frequency bands, leading to severe time distortion and phase shifts in multi-channel signals. Meanwhile, conventional feature extraction methods struggle to suppress complex dynamic background noise in non-stationary and harsh environments, resulting in weak overall system noise resistance and interference immunity. Summary of the Invention

[0004] To overcome the above shortcomings, this invention provides a harmonic structure analysis method, which aims to improve the time distortion caused by the group delay of traditional filters and the dynamic noise interference in complex environments.

[0005] This invention provides the following technical solution: a harmonic structure analysis method, comprising: S1. Obtain the original time-domain acoustic signal to be tested, and sequentially perform framing, windowing and pre-emphasis operations on the original time-domain acoustic signal to obtain multiple frames of short-time preprocessed signals; S2. Input the short-time preprocessed signal of each frame into the auditory filter bank for frequency band decomposition to obtain multiple sub-frequency band signals. Perform phase shift compensation on each sub-frequency band signal using continuously adjustable fractional order parameters to output phase-aligned multi-channel sub-frequency band signals. S3. Perform half-wave rectification, nonlinear amplitude compression and low-pass filtering operations on each channel of the multi-channel sub-band signal in sequence to extract the electrophysiological envelope signal corresponding to each channel. S4. Calculate the autocorrelation function of each of the electrophysiological envelope signals, extract the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency, and extract the initial harmonic peak value that is an integer multiple of the fundamental frequency along the frequency axis. S5. The initial harmonic peak values ​​corresponding to each frame are spliced ​​together in time order to construct a two-dimensional spatiotemporal random matrix. The covariance function of the two-dimensional spatiotemporal random matrix is ​​calculated and eigenvalue decomposition is performed. Orthogonal basis functions are extracted through Karhunen-Loève expansion to reconstruct the concurrent harmonic time-frequency trajectory.

[0006] Preferably, in step S1, the steps of sequentially performing framing, windowing, and pre-emphasis operations include: The original time-domain acoustic signal to be tested is divided into multiple short-time data frames that overlap continuously according to preset frame length and frame shift parameters. The Hamming window function is used to perform time-domain multiplication on each of the short-time data frames; A first-order high-pass differential filter is constructed, and a pre-emphasis filtering operation is performed on each of the windowed short-time data frames to generate the multi-frame short-time preprocessed signal.

[0007] Preferably, in step S2, the step of performing frequency band decomposition on the input auditory filter bank includes: Based on the equivalent rectangular bandwidth scaling formula, multiple center frequency nodes within the target analysis frequency band are calculated. Based on each of the aforementioned center frequency nodes, a Gammatone auditory filter bank is constructed, and a bandwidth control parameter positively correlated with the aforementioned center frequency nodes is set. The short-time preprocessed signals of each frame are subjected to discrete convolution operations with the channel impulse responses of the Gammatone auditory filter bank to extract the multiple sub-band signals.

[0008] Preferably, in step S2, the step of performing phase shift compensation on each of the sub-band signals includes: By introducing the Grünwald-Letnikov operator, a fractional phase shift transfer function is constructed in the complex frequency domain; Calculate the nonlinear group delay distribution of the auditory filter bank for the multiple sub-band signals at each center frequency, and determine the phase compensation difference for each channel; The order parameter of the fractional phase shift transfer function is adjusted according to the phase compensation difference, and the phase spectrum of the multiple sub-band signals is reverse-calibrated to output the phase-aligned multi-channel sub-band signal.

[0009] Preferably, in step S3, the steps of sequentially performing half-wave rectification, nonlinear amplitude compression, and low-pass filtering include: A half-wave rectification and truncation operation is performed on each channel of the phase-aligned multi-channel sub-band signal to remove negative amplitude components and extract the corresponding positive vibration signal of each channel. The positive vibration signals are nonlinearly mapped using a logarithmic compression function to generate amplitude compression signals for each channel. By using a low-pass filter with a cutoff frequency lower than the corresponding center frequency, high-frequency carrier rejection is performed on each of the amplitude compression signals to extract the electrophysiological envelope signal.

[0010] Preferably, in step S4, the step of extracting the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency includes: Calculate the autocorrelation function sequence of the electrophysiological envelope signal of each channel, and sum the autocorrelation function sequences of each channel along the time delay axis to generate a global autocorrelation graph; Within a preset fundamental frequency range, the global autocorrelation graph is searched for extreme points to locate the correlation peak with the largest amplitude; The time delay corresponding to the largest correlation peak is extracted as the basic period, and the reciprocal of the basic period is calculated as the base frequency of the current frame.

[0011] Preferably, in step S4, the step of extracting the initial harmonic peak values ​​along the frequency axis that are integer multiples of the fundamental frequency includes: Based on the determined fundamental frequency, multiple ideal harmonic frequency nodes that are integer multiples of the fundamental frequency are calculated and generated. A tolerance search window of a preset width is constructed around each of the ideal harmonic frequency nodes; Within the auditory spectrum band, the frequency and amplitude of the local energy maximum point within each tolerance search window are extracted and spliced ​​together to form the initial harmonic peak value corresponding to the current frame.

[0012] Preferably, in step S5, the step of calculating the covariance function of the two-dimensional spatiotemporal random matrix includes: The initial harmonic peak values ​​corresponding to multiple consecutive frames are obtained and concatenated as column vectors according to the time sequence to construct the two-dimensional spatiotemporal random matrix. Calculate the mean of each row vector of the two-dimensional spatiotemporal random matrix and perform centering processing; The covariance matrix is ​​obtained by multiplying the centered two-dimensional spatiotemporal random matrix with its transpose and dividing by the total number of frames.

[0013] Preferably, in step S5, the step of reconstructing the concurrent harmonic time-frequency trajectory includes: Singular value decomposition is performed on the covariance function of the two-dimensional spatiotemporal random matrix to extract the eigenvalues ​​and their corresponding orthogonal eigenvectors arranged in descending order of energy contribution rate. Based on a preset cumulative energy contribution rate threshold, the principal component feature vectors whose cumulative energy contribution rate reaches the threshold are extracted, and the Karhunen-Loève orthogonal basis functions are constructed. The two-dimensional spatiotemporal random matrix is ​​projected onto a low-dimensional subspace spanned by the Karhunen-Loève orthogonal basis functions for inverse reconstruction, and the concurrent harmonic time-frequency trajectory is output.

[0014] The present invention has the following beneficial effects: 1. In this invention, by introducing fractional-order calculus operators to compensate for the phase shift of sub-band signals, the time distortion caused by group delay is improved. Furthermore, by combining the Karhunen-Loève expansion to perform orthogonal projection and inverse reconstruction on the two-dimensional spatiotemporal matrix, complex background noise is suppressed, thereby improving the system's noise resistance and anti-interference capability under non-stationary and harsh environments.

[0015] 2. In this invention, the Gammatone auditory filter bank, combined with half-wave rectification and nonlinear logarithmic amplitude compression, is used to simulate the nonlinear dynamic response of the bionic cochlear mechanism to broadband acoustic signals. This breaks through the inherent limitations of traditional linear transformation in time-frequency resolution and achieves high-resolution extraction of weak fundamental frequencies and high-order overtone structures in complex mixed sound sources.

[0016] 3. In this invention, a rigorous computational link including frequency domain tolerance search and multidimensional covariance analysis is constructed. By accumulating global autocorrelation graphs and stitching together the dimensions of continuous multi-frame sequences, discrete time slices are transformed into a spatiotemporal continuous model containing the evolution law of sound waves, filling the analysis blind spot caused by signal mutations and ensuring the continuity and reliability of dynamic tracking of transient targets. Attached Figure Description

[0017] Figure 1 This is a flowchart of a harmonic structure analysis method proposed in this invention; Figure 2 This is a flowchart of the phase shift compensation based on fractional calculus proposed in this invention; Figure 3 This is a flowchart of the concurrent harmonic trajectory reconstruction based on KL expansion proposed in this invention. Detailed Implementation

[0018] The technical solutions in 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, and 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.

[0019] In embodiments of the present invention, the present invention provides a harmonic structure analysis method, such as... Figure 1 As shown, it includes: S1. Obtain the original time-domain acoustic signal to be tested, and perform framing, windowing and pre-emphasis operations on the original time-domain acoustic signal in sequence to obtain multiple frames of short-time preprocessed signals; Further, in step S1, the steps of framing, windowing, and pre-emphasis operations are performed sequentially, including: The original time-domain acoustic signal to be tested is divided into multiple short-time data frames that overlap continuously according to preset frame length and frame shift parameters. The Hamming window function is used to perform time-domain multiplication on each short-time data frame; A first-order high-pass differential filter is constructed, and pre-emphasis filtering is performed on each short-time data frame after windowing to generate multi-frame short-time preprocessed signals.

[0020] Specifically, the discrete sound pressure level sequence acquired by the acoustic sensor is used as the raw time-domain acoustic signal to be measured. The system sampling frequency is set, and preset frame length and frame shift parameters are determined based on the physical steady-state duration of the target sound source. In this scheme, the preferred system sampling frequency is... The original time-domain acoustic signal is divided into multiple consecutive, overlapping short-time data frames. Let the original time-domain acoustic signal be a sequence. ,in This is the index for discrete-time sampling points. The discrete extraction formula for a short-time data frame is expressed as: ; in This indicates the preset number of frame shift samples. This represents the index of a local sampling point within a single frame, and its value range is limited to 1. , This is the preset number of samples per frame length. In this scheme, Preferred , Preferred .

[0021] After extracting each short-time data frame, a time-domain multiplication process is performed on each short-time data frame using a Hamming window function to reduce the energy amplitude at the signal truncation edges. The discrete calculation formula for the Hamming window sequence is as follows: ; The generated Hamming window sequence is multiplied point-by-point in the time domain with each of the defined short-time data frames to generate the windowed short-time data frame signal sequence. The specific formula for the multiplication operation is as follows: ; In the formula Representing the The acoustic amplitude sequence of frames after windowing and smoothing.

[0022] For the windowed data, a first-order high-pass differential filter is constructed to perform pre-emphasis filtering. This differential filter amplifies the high-frequency acoustic components in the signal proportionally by extracting the weighted difference between adjacent acoustic sampling sequence points. The time-domain difference equation for pre-emphasis filtering is: ; In the formula The pre-set pre-weighting coefficient, This is the preprocessed acoustic sequence output after the current single frame has been processed. In this scheme, Parameters are preferred The system iterates through all the divided data frames and performs the above-mentioned time-domain truncation, windowing smoothing, and differential filtering operations, finally outputting a multi-frame short-time preprocessed signal arranged in the original timeline order.

[0023] The above processing effectively reduces spectral leakage caused by acoustic signal segmentation, compensates for the energy loss of high-frequency sound wave components during propagation in physical space, and provides a standard and stable input data source for subsequent frequency band decomposition.

[0024] S2. Input the short-time preprocessed signal of each frame into the auditory filter bank for frequency band decomposition to obtain multiple sub-frequency band signals. Use continuously adjustable fractional order parameters to perform phase shift compensation on each sub-frequency band signal and output phase-aligned multi-channel sub-frequency band signals. Further, in step S2, the step of performing frequency band decomposition on the input auditory filter bank includes: Based on the equivalent rectangular bandwidth scaling formula, multiple center frequency nodes within the target analysis frequency band are calculated. Based on each center frequency node, a Gammatone auditory filter bank is constructed, and a bandwidth control parameter that is positively correlated with the center frequency node is set. Each frame of short-time preprocessed signal is subjected to discrete convolution operation with the channel impulse response of the Gammatone auditory filter bank to extract multiple sub-band signals.

[0025] Furthermore, step S2, the step of performing phase shift compensation on each sub-band signal, includes: By introducing the Grünwald-Letnikov operator, a fractional phase shift transfer function is constructed in the complex frequency domain; Calculate the nonlinear group delay distribution of the auditory filter bank for multiple sub-band signals at each center frequency, and determine the phase compensation difference for each channel; The order parameters of the fractional-order phase shift transfer function are adjusted based on the phase compensation difference to perform inverse calibration of the phase spectrum of multiple sub-band signals, and output phase-aligned multi-channel sub-band signals.

[0026] Specifically, the lowest and highest acoustic signal frequencies of the analysis band are defined, and multiple center frequency nodes within the target analysis band are calculated based on the equivalent rectangular bandwidth scaling formula. In this scheme, the target analysis band is preferably... to The specific nonlinear mapping relationship between the equivalent rectangular bandwidth and the center frequency is expressed as follows: ; In the formula The center frequency node of the filter channel, in Hertz. This represents the bandwidth limit at the corresponding center frequency node. A Gammatone auditory filter bank simulating the frequency response characteristics of the basilar membrane is constructed based on each center frequency node, and a bandwidth control parameter positively correlated with the center frequency node is set. The bandwidth control parameter is calculated as the equivalent rectangular bandwidth of the corresponding center frequency node multiplied by a constant 1.019. The discrete-time domain... The impulse response sequence of each channel is represented as follows: ; In the formula For the first Discrete impulse response sequence of each channel For time sampling point index, To control the gain coefficient of the acoustic output amplitude, The preset fourth-order filter constant, For the first Bandwidth control parameters for each channel, The system sampling frequency of the acoustic sensor acquisition device. For the first The center frequency node corresponding to each channel This represents the initial phase of the signal. In this scheme, Parameters are preferred , Parameters are preferred , Parameters are preferred , Parameters are preferred The input multi-frame short-time preprocessed signal sequence is subjected to one-dimensional discrete convolution operation with the discrete impulse response sequence of each channel of the auditory filter bank to extract multiple sub-band signals covering the entire analysis frequency band.

[0027] The nonlinear group delay distribution of the auditory filter bank at each center frequency for multiple sub-band signals is calculated. The actual acoustic group delay of each channel is obtained by taking the negative derivative of the phase spectrum of the filter bank's transfer function in the complex frequency domain with respect to the angular frequency. The channel with the largest absolute delay value is selected as the time synchronization reference. This reference delay time is subtracted from the actual group delay of each of the other channels to determine the phase compensation difference required for synchronization of each channel.

[0028] A fractional-order phase-shift transfer function with continuous phase adjustment capability in the complex frequency domain is constructed by introducing the Grünwald-Letnikov operator. To ensure that the acoustic signal does not experience amplitude attenuation when adjusting the phase, a frequency domain transfer function in pure phase form is constructed: ; In the formula For the first Fractional phase shift transfer function for each channel For discrete acoustic subband signals, the independent angular frequency variables are... It is a continuously adjustable fractional order parameter. Let be the sign function defined in the frequency domain. Based on the calculated phase compensation difference, the order parameters of the fractional-order phase shift transfer function are adjusted inversely to ensure that the local phase deflection generated at the corresponding center frequency after adjusting the order parameters accurately cancels the time asynchrony caused by the group delay. The formula for adjusting the fractional-order order parameters is: ; In the formula For the first The acoustic center angular frequency corresponding to each channel The phase compensation difference for the corresponding channel obtained above is used. The extracted sub-band signals are transformed to the complex frequency domain via Fourier transform and multiplied by the fractional-order phase shift transfer function with adjusted parameters. This performs nonlinear inverse calibration of the phase spectra of the multiple sub-band signals. Finally, the frequency domain sequence is mapped back to the time domain via inverse Fourier transform, outputting a multi-channel sub-band signal sequence with strictly aligned phases across all channels.

[0029] This step completes the frequency band separation of broadband acoustic features and eliminates the group delay bias introduced by the nonlinear response of the filter through fractional calculus, thus ensuring the time alignment accuracy of the parallel acoustic channel data.

[0030] S3. Perform half-wave rectification, nonlinear amplitude compression and low-pass filtering operations on each channel of the multi-channel sub-band signal in sequence to extract the electrophysiological envelope signal corresponding to each channel. Further, in step S3, the steps of performing half-wave rectification, nonlinear amplitude compression, and low-pass filtering operations in sequence include: Half-wave rectification and truncation operation is performed on each channel of the phase-aligned multi-channel sub-band signal to remove negative amplitude components and extract the corresponding positive vibration signal of each channel; The nonlinear amplitude mapping of each positive vibration signal is performed using a logarithmic compression function to generate the amplitude compression signal corresponding to each channel; By using a low-pass filter with a cutoff frequency lower than the corresponding center frequency, high-frequency carrier rejection is performed on each amplitude compressed signal to extract the electrophysiological envelope signal.

[0031] Specifically, let the acoustic amplitude sequence value of the k-th frequency band channel in the input multi-channel sub-band signal at the m-th discrete-time sampling point be... The half-wave rectification truncation operation is performed independently on each channel of the multi-channel sub-band signal. The truncation operation directly removes the negative numerical physical components from the discrete acoustic sequence, extracting the corresponding positive vibration signal sequence. The specific half-wave rectification mapping operator is expressed as follows: ; In the formula It is the positive acoustic vibration amplitude sequence output by the k-th channel after rectification and truncation.

[0032] After extracting the positive vibration signal sequence, a logarithmic amplitude compression function is constructed to continuously perform nonlinear numerical mapping on each positive vibration signal, thereby compressing the physical dynamic range of high-energy surge noise and amplifying weak acoustic feature components. The nonlinear logarithmic amplitude mapping formula is constructed as follows: ; In the formula The amplitude-compressed signal sequence of each channel generated by the mapping. To limit the global scaling factor of the physical output range of signal energy, To control the nonlinear adjustment parameters of the curvature amplified by the weak acoustic envelope feature, this scheme... Parameters are preferred , Parameters are preferred .

[0033] After generating the amplitude-compressed signals for each channel, a low-pass smoothing filter is connected in series with each signal link to perform high-frequency acoustic carrier rejection. It is strictly stipulated that the physical cutoff frequency of this low-pass filter must be lower than the reference value of the corresponding channel's center frequency node. A first-order differential iterative architecture is used to extract the electrophysiological envelope signal reflecting the slowly varying energy characteristics of the sound source. The filtering iterative operator is configured as follows: ; In the formula To extract the discrete values ​​of the final electrophysiological envelope signal generated at the current sampling time, Let be the envelope physical memory constant of the previous discrete sampling time. This is a smoothing forgetting factor determined by both the preset filter cutoff frequency and the system audio sampling rate. In this scheme, Parameters are preferred .

[0034] This process removes the original high-frequency oscillation components from the acoustics and outputs a baseband envelope feature sequence that can directly characterize the slowly varying energy intensity of the sound wave.

[0035] S4. Calculate the autocorrelation function of each electrophysiological envelope signal, extract the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency, and extract the initial harmonic peak value that is an integer multiple of the fundamental frequency along the frequency axis. Further, in step S4, the step of extracting the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency includes: Calculate the autocorrelation function sequence of the electrophysiological envelope signal of each channel, and sum the autocorrelation function sequences of each channel along the time delay axis to generate a global autocorrelation plot; Within a preset fundamental frequency range, the extreme point search of the global autocorrelation plot is performed to locate the correlation peak with the largest amplitude; Extract the time delay corresponding to the largest correlation peak as the basic period, and calculate the reciprocal of the basic period as the base frequency of the current frame.

[0036] Further, step S4, the step of extracting the initial harmonic peak values ​​that are integer multiples of the fundamental frequency along the frequency axis, includes: Based on a defined fundamental frequency, multiple ideal harmonic frequency nodes that are integer multiples of the fundamental frequency are calculated and generated. A tolerance search window of preset width is constructed around each ideal harmonic frequency node; Within the auditory spectrum band, the frequency and amplitude of the local energy maximum point within each tolerance search window are extracted and spliced ​​together to form the initial harmonic peak value corresponding to the current frame.

[0037] Specifically, the discrete sequences of electrophysiological envelope signals for each channel extracted in the aforementioned steps are obtained. Let the electrophysiological envelope signal of the k-th frequency channel be... , where m is the index of the discrete-time sampling point in the current frame. The discrete autocorrelation function sequence is calculated independently for each channel, and the formula for calculating the autocorrelation function sequence is expressed as: ; In the formula The delay time corresponding to the k-th channel The autocorrelation discrete value, where M is the total number of sampling points in the current acoustic data frame. Let M be the discrete-time delay variable. In this scheme, the preferred M parameter is... The autocorrelation function sequences of all channels are summed at corresponding points along the time delay axis to generate a global autocorrelation map sequence representing the global acoustic periodicity. The formula for synthesizing the global autocorrelation map is: ; In the formula Here, K represents the amplitude sequence of the global autocorrelation plot, and K is the total number of channels in the auditory filter bank. In this scheme, the parameter K is preferably... .

[0038] Based on the physical properties of the target sound source, a lower and upper frequency limit for the fundamental frequency search range are preset. This frequency range is then mapped to a discrete time delay range by combining the system's audio sampling frequency. Within this time delay search interval, a maximum extremum point search is performed on the global autocorrelation graph sequence to locate the significant correlation peak with the highest energy aggregation. The formula for locating the maximum correlation peak is: ; In the formula This represents the time delay corresponding to the maximum correlation peak. and These represent the number of discrete delay sample boundaries corresponding to the upper and lower limits of the fundamental frequency, respectively. In this scheme, the preset fundamental frequency search range is preferably... to Corresponding Parameters are preferred , Parameters are preferred The time delay corresponding to the maximum correlation peak is extracted as the fundamental period. The reciprocal of the fundamental period is calculated by dividing the system sampling frequency by the number of discrete samples in the fundamental period, yielding the exact fundamental frequency of the current acoustic data frame. The fundamental frequency calculation formula is: ; In the formula To determine the target acoustic fundamental frequency in the current frame, The system sampling frequency is the original acoustic signal.

[0039] Using a defined fundamental frequency as the benchmark, multiple ideal harmonic frequency nodes are calculated and generated along the frequency axis. The frequency value of the h-th ideal harmonic frequency node is the product of the fundamental frequency and the corresponding integer order h. A tolerance search window with a set bandwidth is constructed around each ideal harmonic frequency node. Using the auditory spectrum amplitude data obtained by performing a Fourier transform on the current short-time preprocessed signal, local amplitude extrema are retrieved within the real frequency sub-interval covered by each tolerance search window. In this scheme, the bandwidth of the tolerance search window is preferably... Locate the true frequency coordinates and acoustic amplitude of the local energy maximum points within each search window, extract and splice them to form an initial harmonic peak sequence characterizing the acoustic multi-order overtone structure of the current single frame.

[0040] This step enables accurate quantification of the periodic characteristics of complex acoustic signals and effectively extracts the core frequency domain structure information of the target sound source.

[0041] S5. The initial harmonic peak values ​​corresponding to each frame are spliced ​​together in time order to construct a two-dimensional spatiotemporal random matrix. The covariance function of the two-dimensional spatiotemporal random matrix is ​​calculated and eigenvalue decomposition is performed. Orthogonal basis functions are extracted through Karhunen-Loève expansion to reconstruct the concurrent harmonic time-frequency trajectory.

[0042] Further, step S5, the step of calculating the covariance function of the two-dimensional spatiotemporal random matrix, includes: The initial harmonic peak values ​​corresponding to multiple consecutive frames are obtained and concatenated as column vectors according to the time sequence to construct a two-dimensional spatiotemporal random matrix. Calculate the mean of each row vector of a two-dimensional spatiotemporal random matrix and perform centering processing; The covariance matrix is ​​obtained by multiplying the centered two-dimensional spatiotemporal random matrix with its transpose and dividing by the total number of frames.

[0043] Furthermore, step S5, the step of reconstructing the time-frequency trajectory of the concurrent harmonics, includes: Singular value decomposition is performed on the covariance function of the two-dimensional spatiotemporal random matrix to extract the eigenvalues ​​and their corresponding orthogonal eigenvectors arranged in descending order of energy contribution rate. Based on a preset cumulative energy contribution rate threshold, the principal component feature vectors that reach the cumulative energy contribution rate threshold are extracted, and the Karhunen-Loève orthogonal basis functions are constructed. The two-dimensional spatiotemporal random matrix is ​​projected onto a low-dimensional subspace spanned by Karhunen-Loève orthogonal basis functions for inverse reconstruction, and the concurrent harmonic time-frequency trajectory is output.

[0044] Specifically, the initial harmonic peak sequences corresponding to multiple consecutive acoustic data frames are obtained. The total number of consecutively extracted acoustic data frames is set to N, and the total number of extracted harmonic orders in a single frame is H. The initial harmonic peak sequences of each data frame are sequentially used as column vectors and horizontally concatenated to construct a two-dimensional spatiotemporal random matrix P with dimensions H rows and N columns. In this scheme, the parameter N is preferably... The H parameter is preferably Calculate the statistical mean of acoustic amplitude for each row vector in the two-dimensional spatiotemporal random matrix P across different time frames, generating a feature mean column vector. Subtract this feature mean column vector from each column of the original matrix P to perform centering, eliminating the bias effect of static acoustic background energy in subsequent covariance calculations, and outputting the centered two-dimensional spatiotemporal random matrix X.

[0045] The centered two-dimensional spatiotemporal random matrix X is multiplied by its transpose, and each element of the product matrix is ​​divided by the total number of extracted frames N to calculate the covariance matrix, which characterizes the acoustic energy linkage correlation between different harmonic orders. The formula for calculating the covariance matrix is ​​as follows: ; In the formula, C is the acoustic feature covariance matrix of dimension H x H, and X is the centered two-dimensional spatiotemporal random matrix. Let X be the transpose of matrix X.

[0046] Singular value decomposition (SVD) is performed on the calculated acoustic feature covariance matrix C. Since the covariance matrix is ​​a real symmetric matrix, SVD extracts a set of non-negative eigenvalues ​​and a corresponding set of orthogonal eigenvectors. Let the decomposed eigenvalues ​​be arranged in descending order of energy contribution rate. to The feature space matrix formed by the corresponding orthogonal eigenvectors is U. The ratio of the sum of the eigenvalues ​​of the first K principal components to the sum of all eigenvalues ​​is calculated and used as the current cumulative energy contribution rate. Based on a preset cumulative energy contribution rate threshold, the minimum number of principal components K such that the above ratio is greater than or equal to the preset threshold is found. In this scheme, the preferred cumulative energy contribution rate threshold is... Extract the top K column vectors from the feature space matrix U to construct a Karhunen-Loève orthogonal basis function matrix for filtering occasional noise interference. .

[0047] The centered two-dimensional spatiotemporal random matrix X is projected onto a low-dimensional subspace using the constructed orthogonal basis function matrix to extract core features. Then, the transpose of the orthogonal basis function matrix is ​​used to perform an inverse spatial reconstruction operation. Finally, the previously subtracted acoustic feature mean is added back to the reconstructed matrix data, outputting the final concurrent harmonic time-frequency trajectory. The reconstruction operation formula is expressed as: ; In the formula To reconstruct the output concurrent harmonic time-frequency trajectory matrix, The truncated Karhunen-Loève orthogonal basis function matrix is... Let M be the transpose of the orthogonal basis function matrix, and let M be the static background mean matrix formed by copying and expanding the mean vectors of each row of a two-dimensional spatiotemporal random matrix column by column.

[0048] This step removes random noise interference and information redundancy in the original harmonic sequence through orthogonal expansion technology, reconstructing a pure core harmonic evolution trajectory, which effectively improves the anti-interference capability and noise resistance stability of the complex acoustic signal analysis process.

[0049] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method of harmonic structure analysis, characterized by, include: S1. Obtain the original time-domain acoustic signal to be tested, and sequentially perform framing, windowing and pre-emphasis operations on the original time-domain acoustic signal to obtain multiple frames of short-time preprocessed signals; S2. Input the short-time preprocessed signal of each frame into the auditory filter bank for frequency band decomposition to obtain multiple sub-frequency band signals. Perform phase shift compensation on each sub-frequency band signal using continuously adjustable fractional order parameters to output phase-aligned multi-channel sub-frequency band signals. S3. Perform half-wave rectification, nonlinear amplitude compression and low-pass filtering operations on each channel of the multi-channel sub-band signal in sequence to extract the electrophysiological envelope signal corresponding to each channel. S4. Calculate the autocorrelation function of each of the electrophysiological envelope signals, extract the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency, and extract the initial harmonic peak value that is an integer multiple of the fundamental frequency along the frequency axis. S5. The initial harmonic peak values ​​corresponding to each frame are spliced ​​together in time order to construct a two-dimensional spatiotemporal random matrix. The covariance function of the two-dimensional spatiotemporal random matrix is ​​calculated and eigenvalue decomposition is performed. Orthogonal basis functions are extracted through Karhunen-Loève expansion to reconstruct the concurrent harmonic time-frequency trajectory.

2. The method of claim 1, wherein, In step S1, the steps of sequentially performing framing, windowing, and pre-emphasis operations include: The original time-domain acoustic signal to be tested is divided into multiple short-time data frames that overlap continuously according to preset frame length and frame shift parameters. The Hamming window function is used to perform time-domain multiplication on each of the short-time data frames; A first-order high-pass differential filter is constructed, and a pre-emphasis filtering operation is performed on each of the windowed short-time data frames to generate the multi-frame short-time preprocessed signal.

3. The method of claim 1, wherein Step S2, the step of performing frequency band decomposition on the input auditory filter bank, includes: Based on the equivalent rectangular bandwidth scaling formula, multiple center frequency nodes within the target analysis frequency band are calculated. Based on each of the aforementioned center frequency nodes, a Gammatone auditory filter bank is constructed, and a bandwidth control parameter positively correlated with the aforementioned center frequency nodes is set. The short-time preprocessed signals of each frame are subjected to discrete convolution operations with the channel impulse responses of the Gammatone auditory filter bank to extract the multiple sub-band signals.

4. The method of claim 1, wherein Step S2, the step of performing phase shift compensation on each of the sub-band signals, includes: By introducing the Grünwald-Letnikov operator, a fractional phase shift transfer function is constructed in the complex frequency domain; Calculate the nonlinear group delay distribution of the auditory filter bank for the multiple sub-band signals at each center frequency, and determine the phase compensation difference for each channel; The order parameter of the fractional phase shift transfer function is adjusted according to the phase compensation difference, and the phase spectrum of the multiple sub-band signals is reverse-calibrated to output the phase-aligned multi-channel sub-band signal.

5. The method of claim 1, wherein In step S3, the steps of sequentially performing half-wave rectification, nonlinear amplitude compression, and low-pass filtering include: A half-wave rectification and truncation operation is performed on each channel of the phase-aligned multi-channel sub-band signal to remove negative amplitude components and extract the corresponding positive vibration signal of each channel. The positive vibration signals are nonlinearly mapped using a logarithmic compression function to generate amplitude compression signals for each channel. By using a low-pass filter with a cutoff frequency lower than the corresponding center frequency, high-frequency carrier rejection is performed on each of the amplitude compression signals to extract the electrophysiological envelope signal.

6. The harmonic structure analysis method according to claim 1, characterized in that, In step S4, the step of extracting the time delay corresponding to the peak value in the autocorrelation function as the fundamental period to determine the fundamental frequency includes: Calculate the autocorrelation function sequence of the electrophysiological envelope signal of each channel, and sum the autocorrelation function sequences of each channel along the time delay axis to generate a global autocorrelation graph; Within a preset fundamental frequency range, the global autocorrelation graph is searched for extreme points to locate the correlation peak with the largest amplitude; The time delay corresponding to the largest correlation peak is extracted as the basic period, and the reciprocal of the basic period is calculated as the base frequency of the current frame.

7. The harmonic structure analysis method according to claim 1, characterized in that, Step S4, the step of extracting the initial harmonic peak values ​​along the frequency axis that are integer multiples of the fundamental frequency, includes: Based on the determined fundamental frequency, multiple ideal harmonic frequency nodes that are integer multiples of the fundamental frequency are calculated and generated. A tolerance search window of a preset width is constructed around each of the ideal harmonic frequency nodes; Within the auditory spectrum band, the frequency and amplitude of the local energy maximum point within each tolerance search window are extracted and spliced ​​together to form the initial harmonic peak value corresponding to the current frame.

8. The method of claim 1, wherein, Step S5, the step of calculating the covariance function of the two-dimensional spatiotemporal random matrix, includes: The initial harmonic peak values ​​corresponding to multiple consecutive frames are obtained and concatenated as column vectors according to the time sequence to construct the two-dimensional spatiotemporal random matrix. Calculate the mean of each row vector of the two-dimensional spatiotemporal random matrix and perform centering processing; The covariance matrix is ​​obtained by multiplying the centered two-dimensional spatiotemporal random matrix with its transpose and dividing by the total number of frames.

9. The method of claim 1, wherein, Step S5, the step of reconstructing the time-frequency trajectory of the concurrent harmonics, includes: Singular value decomposition is performed on the covariance function of the two-dimensional spatiotemporal random matrix to extract the eigenvalues ​​and their corresponding orthogonal eigenvectors arranged in descending order of energy contribution rate. Based on a preset cumulative energy contribution rate threshold, the principal component feature vectors whose cumulative energy contribution rate reaches the threshold are extracted, and the Karhunen-Loève orthogonal basis functions are constructed. The two-dimensional spatiotemporal random matrix is ​​projected onto a low-dimensional subspace spanned by the Karhunen-Loève orthogonal basis functions for inverse reconstruction, and the concurrent harmonic time-frequency trajectory is output.

Images