A heart sound segmentation method and storage medium based on instantaneous frequency and instantaneous damping

Abstract

The present invention provides a heart sound segmentation method based on instantaneous damping and instantaneous frequency and a storage medium; the method is: collecting the heart sound signal to be segmented, preprocessing the heart sound signal, dividing the signal into frames and adding windows after preprocessing, calculating the power spectrum density envelope vector, instantaneous damping and instantaneous frequency of each frame of heart sound, constructing an observation vector as the input of the LR classifier, and finally using the Viterbi algorithm to decode the LR-HSMM model to obtain the best segmentation sequence. The method adds instantaneous damping, instantaneous frequency, and power spectrum density envelope values ​​of multiple consecutive frames as segmentation features, thereby improving the anti-noise performance of the segmentation method and more accurately describing the energy characteristics of heart sound.

Description

Technical Field

[0001] The present invention relates to the field of medical signal processing, and in particular to a heart sound segmentation method based on instantaneous frequency and instantaneous damping and a storage medium. Background Art

[0002] Heart sound segmentation divides a single heart sound cycle into four stages: the first heart sound S1, systole, second heart sound S2, and diastole; it is one of the important steps to achieve automatic auscultation of heart sounds. Most pathological murmurs occur at a specific time point in the heart sound cycle, so accurately segmenting heart sounds and then determining the time point of murmur occurrence is a prerequisite for correctly identifying pathological murmurs.

[0003] The heart sound segmentation method based on the Logistic Regression Hidden Semi-Markov Model (LR-HSMM model) is one of the algorithms with the best segmentation effect. Existing LR-HSMM models often extract features such as the heart sound homomorphic envelope, Hilbert envelope, power spectral density envelope and wavelet envelope to construct the observation sequence of the model. When the parameters are set reasonably, these features can effectively reflect the changes in the intensity of heart sounds over time, which helps the LR-HSMM model to determine the positions of the first and second heart sounds. However, when the time-frequency distribution of the ambient noise overlaps with the first or second heart sound, the calculation results of these features will be interfered by the noise and obtain abnormal envelope amplitudes, which will eventually lead to segmentation errors.

[0004] Therefore, considering that the heart sound signal acquisition process is easily interfered by various noises, how to extract appropriate heart sound features to construct the observation sequence so that the segmentation method can still have a high accuracy even under conditions of various environmental noise interferences is a problem that needs to be solved. Summary of the invention

[0005] In order to improve the anti-noise ability of existing heart sound segmentation methods, the present invention provides a heart sound segmentation method and storage medium based on instantaneous frequency and instantaneous damping; the segmentation method adds instantaneous damping, instantaneous frequency, and power spectral density envelope values ​​of multiple consecutive frames as segmentation features, thereby improving the anti-noise performance of the segmentation method, and can more accurately characterize the energy characteristics of heart sounds and improve the segmentation accuracy.

[0006] In order to achieve the above object, the present invention provides the following technical solution: a heart sound segmentation method based on instantaneous frequency and instantaneous damping, characterized in that it includes the following steps:

[0007] Y1: For a single subject, collect the heart sound signal to be segmented;

[0008] Y2: Preprocess the segmented heart sound signal to obtain the heart sound signal s(l); after the heart sound signal s(l) is framed and windowed, the power spectrum density envelope vector E(m), instantaneous damping η(m) and instantaneous frequency f(m) of the m-th frame of heart sound are calculated; construct the observation vector F(m) of the m-th frame of heart sound = [E(M), E(m-1), E(m-2), η(m), f(m)] T ;

[0009] Y3: Arrange the observation vectors F(1),…,F(m) of the first m frames of heart sounds in time sequence to construct the observation sequence O m =[F(1),F(2),…,F(m)];

[0010] Y4: The observation sequence O m As the input of the LR classifier in the LR-HSMM model, the Viterbi algorithm is used to decode the LR-HSMM model to obtain the optimal segmentation sequence S of the heart sound signal. m .

[0011] Preferably, said Y2 comprises the following sub-steps:

[0012] Y2.1: Downsample the segmented heart sound signal to f s ;

[0013] Y2.2: The downsampled heart sound signal is bandpass filtered to obtain the heart sound signal s(l); where l = 1, 2, ..., L s , L s is the number of sampling points of the heart sound signal after downsampling;

[0014] Y2.3: Frame and window the heart sound signal s(l) to obtain the mth frame of windowed heart sound y w (m) and windowed heart sound z w (m):

[0015] y w (m) = [s((m-1)Δ+1)·w y (1),…,s((m-1)Δ+L y )·w y (L y )];

[0016] z w (m) = [s((m-1)Δ+1)·w z (1),…,s((m-1)Δ+L z )·w z (L z )];

[0017] Among them, y w (m) is used to extract the power spectral density envelope, zw (m) is used to extract the instantaneous damping ratio and instantaneous natural oscillation frequency; L y is the frame length; L z is the frame length; w y (·) and w z (·) is the window function; Δ is the frame shift representing the time difference between the starting positions of the previous frame and the next frame; y w (m) and z w The frame shifts Δ of (m) are equal;

[0018] Y2.4: Calculate the windowed heart sound y of the mth frame w The power spectral density envelope E of (m) E (m), and normalized to the power spectral density envelope vector E(m);

[0019] Y2.5: Calculate the windowed heart sound z of the mth frame w (m) instantaneous damping ratio η η (m) and the instantaneous natural oscillation frequency f f (m) and normalized to the instantaneous damping η(m) and the instantaneous frequency f(m).

[0020] Y2.6: Construct the observation vector of the mth frame heart sound F(m) = [E(m), E(m-1), E(m-2), η(m), f(m)] T .

[0021] Preferably, in Y2.4, the method for obtaining the power spectrum density envelope vector E(m) is: calculating the windowed heart sound y of the mth frame w The power spectrum matrix Y(m,k) of (m), where k = 1, 2, ..., L y is the frequency number; calculate the power spectrum density envelope E E (m):

[0022]

[0023] in, f L and f H They are the set lower and upper frequency limits respectively; Indicates rounding to the nearest integer;

[0024] Calculate the windowed heart sound y of the mth frame w The power spectral density envelope vector E(m) of (m):

[0025]

[0026] Among them, μ E is the power spectral density envelope E of all previous m frames E The mean of (m); Eis the power spectral density envelope E of all previous m frames E The standard deviation of (m).

[0027] Preferably, in Y2.5, the instantaneous damping η(m) and the instantaneous frequency f(m) are obtained by: w (m) is subjected to a second-order linear prediction analysis to obtain a second-order full-pole model, and the two poles are z1(m) and z2(m); the instantaneous damping ratio η is calculated through the amplitude-angle relationship η (m) and the instantaneous natural oscillation frequency f f (m):

[0028]

[0029] Among them, |z1(μ)| represents the modulus of the pole z1(m); ∠z1(m) represents the magnitude of the argument of the pole z1(m) (in radians), f s is the sampling frequency of the heart sound signal s(l);

[0030] The solution of the linear prediction analysis model is implemented using the Levinson-Durbin algorithm.

[0031] Calculate the windowed heart sound z of the mth frame w The instantaneous damping η(m) and instantaneous frequency f(n) of (m):

[0032]

[0033] Among them, μ η and μ f are all instantaneous damping η of the first m frames respectively η (m) and instantaneous frequency f f The mean of (m); ∑ and ∑ f are all instantaneous damping η of the first M frames respectively η (m) and instantaneous frequency f f The standard deviation of (m).

[0034] Preferably, in Y4, the LR classifier refers to a trained LR classifier;

[0035] The LR classifier training process refers to:

[0036] The four-class LR classifier is established using the Logistic function to fit the conditional probability P(q(m)=ε j |F(m)); The input of LR classifier is the observation vector of heart sound signal F(m)=[E(m),E(m-1),E(m-2),η(m),f(m)] T ,ω=[ω1,ω2,ω3,ω4,ω5] Tand ω0 are the unknown parameters of the Logistic function;

[0037] Collect heart sound signals from P subjects, and manually cut out heart sound signals of N cardiac cycles for each subject. One cardiac cycle represents one sample, and construct a training set with P·N samples; for each sample, construct the observation vector of the sample according to Y2 Where 1≤h≤P·N; using the sample observation vector F h (m) Training is performed based on the gradient descent method to train the unknown parameters ω and ω0 of the LR classifier to obtain the optimal weight parameter ω * and

[0038] Preferably, the parameters set and calculated for the LR-HSMM model include the state transition probability matrix A, the observation probability matrix B, the initial state probability vector π, the probability density function p of the duration j (d) Observation sequence O m The probability of generation P(O m ).

[0039] Preferably, the LR-HSMM model, setting and calculating parameters comprises the following sub-steps:

[0040] Z1: Calculate cardiac cycle T a The length of time T that the heart sounds are in S1 and systolic states b ;

[0041] Z2: Set the initial state probability vector π, the state transition probability matrix A = {a ij};

[0042] Z3: Statistics of F of all training samples h The mean vector μ of (m) F and the covariance matrix ∑ F , construct observation sequence O m , calculate and generate the observation sequence O m The probability P(O m ) and the observation probability matrix B = {b j (O m )}; where b j (O m ) is the observation probability;

[0043] Z4: Calculate the mean μ of the duration of the jth state of the heart sound j and standard deviation∑ j , calculate the probability density function p of the duration j (d); where j = 1, 2,…, 4.

[0044] Preferably, in Z2, the initial state probability vector π is set to The state transition probability matrix A is set as:

[0045]

[0046] With q(m-1)=ε i As input, estimate the m-th frame heart sound state q(m) = ε j The conditional probability a ij :

[0047] a ij =P(q(m)=ε j |q(m-1)=ε i )(1≤i,j≤4)

[0048] Among them, q(m)∈{ε1,ε2,ε3,ε4} and q(m-1)∈{ε1,ε2,ε3,ε4} are the states of heart sounds at time m and time m-1 respectively; ε1,ε2,ε3,ε4 represent the four states of heart sounds respectively; ε1 indicates that the heart sound is in the S1 state, ε2 indicates that the heart sound is in the systolic state, ε3 indicates that the heart sound is in the S2 state, and ε4 indicates that the heart sound is in the diastolic state.

[0049] Preferably, in Z3, the observation vectors F(1), ..., F(m) of the first m frames of heart sounds are arranged in time to construct an observation sequence O m =[F(1),F(2),…,F(m)];

[0050] Observation probability b j (O m ) is calculated as:

[0051]

[0052] in, is the estimated heart sound state of the mth frame q(m) = ε j The conditional probability of generating the observation sequence O m The probability P(O m ) is estimated using a multivariate Gaussian model, and the corresponding mean vector and covariance matrix are μ F and∑ F ; The state is ε j The probability P(ε j ) is equal to the initial state probability of Z2

[0053] In Z4, the probability density function p of the duration j (d) is:

[0054]

[0055] Among them, j=1 represents that the heart sound is in the S1 state; j=2 represents that the heart sound is in the systolic state; j=3 represents that the heart sound is in the S2 state; j=4 represents that the heart sound is in the diastolic state.

[0056] The principle of the segmentation method of the present invention is:

[0057] Since the energy of the first heart sound in the S1 state and the second heart sound in the S2 state is mainly distributed in the frequency band below 150 Hz, and there is a peak near 50 Hz, the existing envelope calculation method usually uses f L to f H Frequency range (usually f L =40Hz, f H =60Hz) calculate the mean of the power spectrum Y(m,k) as the power spectrum density envelope E(m).

[0058] The present invention takes into account the characteristics of the intensity or energy of heart sounds changing over time in this specific frequency band (40-60Hz): the power spectrum density envelope values ​​of heart sound segments in the systolic and diastolic stages are small, while the power spectrum density envelope values ​​of heart sound segments in the S1 and S2 states are large. Therefore, using the power spectrum density envelope values ​​of adjacent frames can more accurately describe the characteristics of the heart sound energy changing over time than using only the single frame envelope value E(m). Therefore, the present invention uses the power spectrum density envelope vector [E(m), E(m-1), E(m-2)] T As part of the input to the LR classifier.

[0059] The first heart sound, the second heart sound, various abnormal heart sounds and various cardiovascular murmur signals are all generated by the blood flow impacting the valve, so there is a specific correlation between the heart sound values ​​at different times, and they can be approximately described as attenuated oscillation signals; and the damping ratio and natural oscillation frequency of the first and second heart sounds are quite different from the damping ratio and natural oscillation frequency of abnormal heart sounds and cardiovascular murmurs, so they can also be used as features for heart sound segmentation. The present invention calculates the damping ratio and natural oscillation frequency in a less computationally intensive manner, performs a second-order linear prediction analysis on a single-frame heart sound signal to obtain a second-order full-pole model, and then calculates the instantaneous damping η(m) and instantaneous frequency f(m) based on the poles of the second-order full-pole model. Therefore, Y2.5 of the present invention uses the instantaneous damping η(m) and instantaneous frequency f(m) based on linear prediction coding as another part of the input of the LR classifier.

[0060] Since the characteristics of the heart sound signal remain basically unchanged in a short period of time and can be regarded as a short-term steady-state process, the present invention Y2.3 frames and windows the heart sound signal s(l) to ensure its short-term stable characteristics; the framing adopts the method of overlapping two adjacent frames to make a smooth transition between the two frames and prevent discontinuous changes in the signal characteristics.

[0061] A storage medium, characterized in that: the storage medium stores a computer program, and the computer program can be used by a processor to implement the above-mentioned heart sound segmentation method based on instantaneous damping and instantaneous frequency.

[0062] Compared with the prior art, the present invention has the following beneficial results:

[0063] 1. Based on the classic heart sound segmentation method based on the LR-HSMM model, the present invention adds instantaneous damping, instantaneous frequency, and power spectrum density envelope values ​​of multiple consecutive frames as segmentation features according to the mechanism and characteristics of heart sound signals, thereby improving the anti-noise performance of the segmentation method and more accurately describing the energy characteristics of heart sounds. The segmentation method has a small amount of calculation and is suitable for heart sound segmentation in environments with low signal-to-noise ratio and complex noise types;

[0064] 2. The present invention uses the power spectrum density envelope value of multiple consecutive frames as the segmentation feature, which can more accurately describe the characteristics of the change of heart sound energy over time;

[0065] 3. Aiming at the correlation characteristics of heart sounds, the present invention extracts instantaneous damping and instantaneous frequency from the heart sound signal as new segmentation features, which is conducive to suppressing noises such as music, metal impact and nail friction, and still has a good segmentation accuracy in an environment with a low signal-to-noise ratio;

[0066] 4. The present invention calculates instantaneous frequency and instantaneous damping based on a second-order linear prediction model, and the amount of calculation is relatively small. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 This is a flow chart of a heart sound segmentation method according to Embodiment 1 of the present invention;

[0068] Figure 2(a) to Figure 2(f) They are respectively a time domain graph of noisy heart sounds of a healthy subject according to the first embodiment of the present invention, a spectrum graph, a time domain graph of power spectrum density envelope, a time domain graph of instantaneous damping ratio, and a time domain graph of instantaneous natural frequency;

[0069] Figure 3(a) to Figure 3(c) They are respectively a time domain curve diagram of noisy heart sounds of a healthy subject in Example 1 of the present invention, a state sequence comparison diagram of an algorithm segmentation result, and a manual labeling result. DETAILED DESCRIPTION

[0070] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods, but the present invention is not limited to these embodiments.

[0071] Embodiment 1

[0072] This embodiment provides a heart sound segmentation method based on instantaneous damping and instantaneous frequency. The process is as follows: Figure 1As shown, the following steps are included:

[0073] Y1: For a single subject, collect the heart sound signal to be segmented;

[0074] Y2: Downsample the heart sound signal to be segmented to f s =1000Hz; after bandpass filtering the downsampled heart sound, the heart sound signal s(l) is obtained, and the passband is [25Hz,400Hz]; where l=1,2,…,L s , L s is the number of sampling points of heart sounds after downsampling;

[0075] The heart sound signal s(l) is framed and windowed to obtain the mth frame of windowed heart sound y w (m) and windowed heart sound z w (m):

[0076] y w (m) = [s((m-1)Δ+1)·w y (1),…,s((m-1)Δ+L y )·w y (L y )];

[0077] z w (m) = [s((m-1)Δ+1)·w z (1),…,s((m-1)Δ+L z )·w z (L z )];

[0078] Among them, y w (m) is used to extract the power spectral density envelope, and the frame length is set to 25ms, that is, L y =25; z w (m) is used to extract the instantaneous damping ratio and instantaneous natural oscillation frequency. The frame length is set to 100ms, i.e., L z =100; frame shift (the time difference between the starting position of the previous frame and the next frame) is unified to 10ms, that is, Δ=10; w y (·) and w z (·) is the Hamming window function;

[0079] Calculate the windowed heart sound y of the mth frame w The power spectrum matrix Y(m,k) of (m), where k = 1, 2, ..., L y is the frequency number; calculate the power spectrum density envelope E E (m):

[0080]

[0081] in, f L =40Hz and f H =60Hz are the lower and upper frequency settings respectively, Indicates rounding to the nearest integer;

[0082] First calculate the power spectral density envelope E of a single subject E The mean value μ of (m) E and standard deviation∑ E , and then calculate the windowed heart sound y of the mth frame w The normalized power spectral density envelope vector E(m) of (m):

[0083]

[0084] Add window heart sound z to the mth frame w (m) Perform second-order linear prediction coding to obtain a second-order full-pole model H(z):

[0085]

[0086] Among them, the coefficient a x and gain G are the parameters of the obtained model, and p = 2 is the model order;

[0087] Through the second-order full-pole model H(z), we get two poles, namely:

[0088]

[0089] Among them, z1 and z2 are a pair of poles with a conjugate relationship, that is, |z1(m)|=|z2(m)|, ∠z1(m)=-∠z2(m);

[0090] Instantaneous damping ratio η η (m) and instantaneous natural frequency f f (m) is calculated as:

[0091]

[0092] Among them, |z1(m)| represents the modulus of the pole z1(m), ∠z1(m) represents the angle of the pole z1(m) (in radians), f s is the sampling frequency of the heart sound signal s(l), f s =1000Hz;

[0093] Calculate the instantaneous damping η for a single subject η (m) and instantaneous frequency f f The mean value μ of (m) η , μ f and standard deviation∑ η ,∑f , and then calculate the windowed heart sound z of the mth frame w The instantaneous damping η(m) and instantaneous frequency f(m) of (m):

[0094]

[0095] Construct the observation vector of the mth frame heart sound F(m) = [E(m), E(m-1), E(m-2), η(m), f(m)] T .

[0096] Taking the heart sound samples of healthy subjects with music noise as an example, Figure 2(a) is the time domain curve of the noisy heart sound, Figure 2(b) is the spectrum diagram of the noisy heart sound, Figure 2(c) is the artificially marked state sequence diagram of the noisy heart sound, Figure 2(d) is the time domain curve of the power spectral density envelope of the noisy heart sound, Figure 2(e) is the time domain curve of the damping coefficient of the noisy heart sound, and Figure 2(f) is the time domain curve of the natural frequency of the noisy heart sound.

[0097] Y3: Arrange the observation vectors F(1),…,F(m) of the first m frames of heart sounds in time sequence to construct the observation sequence O m =[F(1),F(2),…,F(m)];

[0098] Y4: The observation sequence O m As the input of the LR classifier in the LR-HSMM model, the Viterbi algorithm is used to decode the LR-HSMM model to obtain the optimal segmentation sequence S of the heart sound signal. m .

[0099] FIG3(a) is a time domain graph of noisy heart sounds, FIG3(b) is a state sequence diagram of the algorithm segmentation results of noisy heart sounds, and FIG3(c) is a state sequence diagram of the manual labeling results of noisy heart sounds.

[0100] LR-HSMM model, the parameters set and calculated include state transition probability matrix A, observation probability matrix B, initial state probability vector π, and duration probability density function p j (d) Observation sequence O m The probability of generation P(O m ). Setting and calculating parameters includes the following sub-steps:

[0101] Z1: Calculate cardiac cycle T a The length of time T that the heart sounds are in S1 and systolic states b :First calculate the energy spectrum density envelope E(m) of the heart sound signal s(l), and then calculate the autocorrelation function r of E(m) E , search for r in the range [500ms,1200ms] E The time length corresponding to the maximum value, that is, the cardiac cycle Ta ; Search [200ms, 0.5·T a ] within the range of r E The time length corresponding to the maximum value, that is, the time length T b .

[0102] Z2: The initial state probability vector π is set to The state transition probability matrix A is set as:

[0103]

[0104] With q(m-1)=ε i As input, estimate the m-th frame heart sound state q(m) = ε j The conditional probability a ij :

[0105] a ij =P(q(m)=ε j |q(m-1)=ε i )(1≤i,j≤4)

[0106] Among them, q(m)∈{ε1,ε2,ε3,ε4} and q(m-1)∈{ε1,ε2,ε3,ε4} are the states of heart sounds at time m and time m-1 respectively; ε1,ε2,ε3,ε4 represent the four states of heart sounds respectively; ε1 indicates that the heart sound is in the S1 state, ε2 indicates that the heart sound is in the systolic state, ε3 indicates that the heart sound is in the S2 state, and ε4 indicates that the heart sound is in the diastolic state.

[0107] Z3: Calculate the observation vectors F(1),…,F(m) of the previous m frames of heart sounds based on Y2, arrange them in time sequence, and construct the observation sequence O m =[F(1),F(2),…,F(m)];

[0108] According to the F of all the above training samples h (m) The statistical mean vector μ F and the covariance matrix ∑ F , calculate and generate the observation sequence O m The probability P(O m ) and the observation probability matrix B = {b j (O m )}; Observation probability b j (O m ) is calculated as:

[0109]

[0110] in, The m-th frame heart sound state q(m) is estimated by the LR classifier jThe conditional probability of generating the observation sequence O m The probability P(O m ) is estimated using a multivariate Gaussian model, and the corresponding mean vector and covariance matrix are μ F and∑ F ; The state is ε j The probability P(ε j ) is equal to the initial state probability of Z2

[0111] Z4: Calculate the mean μ of the duration of the jth state of the heart sound j and standard deviation∑ j , calculate the probability density function p of the duration j (d); mean μ j and standard deviation Σ j The calculation formula is:

[0112] μ1=122ms;∑1=22ma

[0113] μ2=T b -μ1;∑2=25ms

[0114] μ3=92ms;Σ3=22ms

[0115] μ4=T a -T b -μ3;∑4=0.07·μ4+6ms

[0116] The probability density function of duration p j (d) The calculation formula is:

[0117]

[0118] Among them, j=1 represents that the heart sound is in the S1 state; j=2 represents that the heart sound is in the systolic state; j=3 represents that the heart sound is in the S2 state; j=4 represents that the heart sound is in the diastolic state.

[0119] The LR-HSMM model refers to a model in which an LR classifier is added to the HSMM model; the LR classifier refers to an LR classifier that has been trained; and the LR classifier training process refers to:

[0120] The four-class LR classifier is established using the Logistic function to fit the conditional probability P(q(m)=ε j |F(m)); The input of LR classifier is the observation vector of heart sound signal F(m)=[E(m),E(m-1),E(m-2),η(m),f(m)] T ;Set ω=[ω1,ω2,ω3,ω4,ω5] T and ω0 are the unknown parameters of the Logistic function;

[0121] The heart sound signals of 50 subjects were collected, and the heart sound signals of 4 cardiac cycles were manually cut out for each subject. One cardiac cycle represents one sample, and a training set of 200 samples was constructed; for each sample, according to Y2, the observation vector of the sample was constructed Where 1≤h≤200; at the same time, manually judge the state of the heart sound of the frame. If the heart sound is in the S1 state, then the state q h (m) = ε1; if the heart sound is in the systolic state, then state q h (m) = ε2; if the heart sound is in S2 state, then state q h (m) = ε3; if the heart sound is in the diastolic state, then state q h (m) = ε4.

[0122] Using all training samples {(F1(1),q1(1)),…,(F1(m),q1(m)),…,(F h (1),q h (1)),…,(F h (m),q h (m)),…} to train the undetermined parameters ω and ω0 of the LR classifier based on the gradient descent method to obtain the optimal weight parameter ω * and Statistics of F of 200 training samples h The mean vector μ of (m) F and the covariance matrix ∑ F .

[0123] Embodiment 2

[0124] This embodiment provides a storage medium, wherein the storage medium stores a computer program, and the computer program can be used by a processor to implement the heart sound segmentation method based on instantaneous damping and instantaneous frequency described in the first embodiment.

[0125] The above embodiments are merely examples of the spirit of the present invention, but the implementation methods of the present invention are not limited to the above embodiments. Any other modifications, substitutions, combinations and variations that do not deviate from the spirit and principle of the present invention are included in the scope of the claims of the present invention.

Claims

1. A heart sound segmentation method based on instantaneous frequency and instantaneous damping, characterized in that: The following steps are involved: Y1: For a single subject, collect the heart sound signal to be segmented; Y2: Preprocess the segmented heart sound signal to obtain the heart sound signal s(l); after the heart sound signal s(l) is framed and windowed, the power spectrum density envelope vector E(m), instantaneous damping η(m) and instantaneous frequency f(m) of the m-th frame of heart sound are calculated; construct the observation vector F(m) of the m-th frame of heart sound = [E(m), E(m-1), E(m-2), η(m), f(m)] T ; Y3: Arrange the observation vectors F(1),…,F(m) of the first m frames of heart sounds in time sequence to construct the observation sequence O m =[F(1),F(2),…,F(m)]; Y4: The observation sequence O m As the input of the LR classifier in the LR-HSMM model, the Viterbi algorithm is used to decode the LR-HSMM model to obtain the optimal segmentation sequence S of the heart sound signal. m ; Said Y2 comprises the following steps: Y2.1: Downsample the segmented heart sound signal to f s ; Y2.2: The downsampled heart sound signal is bandpass filtered to obtain the heart sound signal s(l); where l = 1, 2, ..., L s , L s is the number of sampling points of the heart sound signal after downsampling; Y2.3: Frame and window the heart sound signal s(l) to obtain the mth frame of windowed heart sound y w (m) and windowed heart sound z w (m): y w (m)=[s((m-1)Δ+1)·w y (1),…,s((m-1)Δ+L y )·w y (L y )]; z w (m)=[s((m-1)Δ+1)·w z (1),…,s((m-1)Δ+L z )·w z (L z )]; Among them, y w (m) is used to extract the power spectral density envelope, z w (m) is used to extract the instantaneous damping ratio and instantaneous natural oscillation frequency; L y is the frame length; L z is the frame length; w y (·) and w z (·) is the window function; Δ is the frame shift representing the time difference between the starting positions of the previous frame and the next frame; y w (m) and z w The frame shifts Δ of (m) are equal; Y2.4: Calculate the windowed heart sound y of the mth frame w The power spectral density envelope E of (m) E (m), and normalized to the power spectral density envelope vector E(m); Y2.5: Calculate the windowed heart sound z of the mth frame w (m) instantaneous damping ratio η η (m) and the instantaneous natural oscillation frequency f f (m), and normalized to the instantaneous damping η(m) and the instantaneous frequency f(m); Y2.6: Construct the observation vector of the mth frame heart sound F(m)=[E(m),E(m-1),E(m-2),η(m),f(m)] T ; The method for obtaining Y2.5, instantaneous damping η(m) and instantaneous frequency f(m) is: add the window heart sound z to the mth frame w (m) is subjected to a second-order linear prediction analysis to obtain a second-order full-pole model, and the two poles are z1(m) and z2(m); the instantaneous damping ratio η is calculated through the amplitude-angle relationship η (m) and the instantaneous natural oscillation frequency f f (m); Among them, |z1(m)| represents the modulus of the pole z1(m); ∠z1(m) represents the magnitude of the argument of the pole z1(m); f s is the sampling frequency of the heart sound signal s(l); Calculate the windowed heart sound z of the mth frame w The instantaneous damping η(m) and instantaneous frequency f(m) of (m): Among them, μ η and μ f are all instantaneous damping η of the first m frames respectively η (m) and instantaneous frequency f f The mean of (m); Σ η and Σ f are all instantaneous damping η of the first m frames respectively η (m) and instantaneous frequency f f The standard deviation of (m).

2. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 1, characterized in that: In Y2.4, the power spectrum density envelope vector E(m) is obtained by calculating the windowed heart sound y of the mth frame: w The power spectrum matrix Y(m,k) of (m), where k = 1, 2, ..., L y is the frequency number; calculate the power spectrum density envelope E E (m): in, f L and f H They are the set lower and upper frequency limits respectively; Indicates rounding to the nearest integer; Calculate the windowed heart sound y of the mth frame w The power spectral density envelope vector E(m) of (m): Among them, μ E is the power spectral density envelope E of all previous m frames E The mean of (m); E is the power spectral density envelope E of all previous m frames E The standard deviation of (m).

3. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 1, characterized in that: In Y4, the LR classifier refers to a trained LR classifier; The LR classifier training process refers to: The four-class LR classifier is established using the Logistic function to fit the conditional probability P(q(m)=ε j |F(m)); The input of LR classifier is the observation vector of heart sound signal F(m)=[E(m),E(m-1),E(m-2),η(m),f(m)] T ;Set ω=[ω1,ω2,ω3,ω4,ω5] T and ω0 are the unknown parameters of the Logistic function; Collect heart sound signals from P subjects, and manually cut out heart sound signals of N cardiac cycles for each subject. One cardiac cycle represents one sample, and construct a training set with P·N samples; for each sample, construct the observation vector of the sample according to Y2 Where 1≤h≤P·N; using the sample observation vector F h (m) Training is performed based on the gradient descent method to train the unknown parameters ω and ω0 of the LR classifier to obtain the optimal weight parameter ω * and 4. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 1, characterized in that: The parameters set and calculated for the LR-HSMM model include the state transition probability matrix A, the observation probability matrix B, the initial state probability vector π, and the probability density function p of the duration. j (d) Observation sequence O m The probability of generation P(O m ).

5. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 4, characterized in that: The LR-HSMM model, setting and calculating parameters include the following sub-steps: Z1: Calculate cardiac cycle T a The length of time T that the heart sound is in S1 state and systolic state b ; Z2: Set the initial state probability vector π, the state transition probability matrix A = {a ij }; Z3: Statistics of F of all training samples h The mean vector μ of (m) F and the covariance matrix ∑ F , calculate and generate the observation sequence O m The probability P(O m ) and the observation probability matrix B = {b j (O m )}; where b j (O m ) is the observation probability; Z4: Calculate the mean μ of the duration of the jth state of the heart sound j and standard deviation∑ j , calculate the probability density function p of the duration j (d); where j = 1, 2,…, 4.

6. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 5, characterized in that: In Z2, the initial state probability vector π is set to The state transition probability matrix A is set as: With q(m-1)=ε i As input, estimate the m-th frame heart sound state q(m) = ε j The conditional probability a ij : a ij =P(q(m)=ε j |q(m-1)=ε i ), 1≤i,j≤4, where q(m)∈{ε1,ε2,ε3,ε4} and q(m-1)∈{ε1,ε2,ε3,ε4} are the heart sound states at time m and time m-1 respectively; ε1,ε2,ε3,ε4 represent the four states of heart sound respectively; ε1 indicates that the heart sound is in the S1 state, ε2 indicates that the heart sound is in the systolic state, ε3 indicates that the heart sound is in the S2 state, and ε4 indicates that the heart sound is in the diastolic state.

7. The heart sound segmentation method based on instantaneous frequency and instantaneous damping according to claim 6, characterized in that: In Z3, the observation vectors F(1),…,F(m) of the first m frames of heart sounds are arranged in time sequence to construct the observation sequence O m =[F(1),F(2),…,F(m)]; Observation probability b j (O m )for: in, is the estimated heart sound state of the mth frame q(m) = ε j The conditional probability of generating the observation sequence O m The probability P(O m ) is estimated using a multivariate Gaussian model, and the corresponding mean vector and covariance matrix are μ F and∑ F ; The state is ε j The probability P(ε j ) is equal to the initial state probability of Z2 In Z4, the probability density function p of the duration j (d) is: Among them, j=1 represents that the heart sound is in the S1 state; j=2 represents that the heart sound is in the systolic state; j=3 represents that the heart sound is in the S2 state; j=4 represents that the heart sound is in the diastolic state.

8. A storage medium, characterized in that: The storage medium stores a computer program, and the computer program can be used by a processor to implement the heart sound segmentation method based on instantaneous frequency and instantaneous damping as described in any one of claims 1 to 7 above.

Images