Quantitative evaluation method for influence of rhythm audio pulse wave synchronization on cardiovascular system

Abstract

The invention discloses a quantitative evaluation method for the influence of rhythm audio pulse wave synchronization on heart and blood vessels. The method comprises the following steps: processing a single-tone audio of a tone library to obtain a rhythm audio; pulse waves before and after audio listening of the listening intervention group, before and after audio playing of the playing intervention group and pulse waves of the quiet sitting control group are collected respectively; then the pulse waves and the rhythm audio are preprocessed; extracting features related to cardiovascular health in the pulse waves; calculating synchronization features between the pulse waves collected each time and the rhythm audio used correspondingly; and finally, constructing a regression formula taking the synchronization characteristic as an independent variable and cardiovascular health as a dependent variable by utilizing a machine learning method so as to quantitatively analyze the influence of the synchronization degree on the cardiovascular system. According to the influence of the synchronization degree of the rhythm audio and the pulse wave on the cardiovascular system, personalized rhythm audio intervention method recommendation is provided for the user.

Description

Technical Field

[0001] The present invention relates to the field of signal processing, and in particular to a quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization. Background Art

[0002] In recent years, with the deepening integration of music and biosignal research, a growing number of studies have begun to explore the relationship between music and human physiological signals, particularly in the field of cardiovascular health management. Research on the synchronization effect between rhythmic audio and pulse wave signals aims to explore the synchronization relationship between rhythmic audio and human physiological signals. Pulse wave signals contain many important indicators reflecting the human cardiovascular state, and their changes can intuitively reflect an individual's physiological condition. By analyzing the synchronization relationship between rhythmic audio and pulse wave signals, synchronization features can be obtained and correlated with cardiovascular health characteristics such as pulse rate variability. This process requires not only signal processing techniques to extract pulse wave features, but also machine learning algorithms to model the synchronization effect and quantify the impact of audio on cardiovascular health.

[0003] Analyzing the synchronization between rhythmic audio and pulse wave signals can provide new approaches for personalized health interventions. However, existing music intervention methods, which mostly rely on traditional music therapy, are unable to accurately recommend interventions based on an individual's physiological state in real time. Therefore, achieving more precise and personalized health interventions has become an urgent problem. Summary of the Invention

[0004] The purpose of the present invention is to solve the above-mentioned defects in the prior art and to propose a quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization.

[0005] The purpose of the present invention can be achieved by taking the following technical solutions:

[0006] A method for quantitatively evaluating the cardiovascular effects of rhythmic audio pulse wave synchronization includes the following steps:

[0007] S1. Perform endpoint detection on the monophonic audio of Chinese drums and then generate multiple sets of rhythmic audio based on music theory.

[0008] S2. Divide the subjects into three groups: a listening intervention group, a playing intervention group, and a sitting meditation control group, and then conduct an experiment; wherein the subjects in the listening intervention group are required to listen to the rhythm audio, and the pulse wave signals of the subjects are collected before and after listening; the subjects in the playing intervention group are required to play the rhythm audio, and the pulse wave signals of the subjects are collected before and after playing; and the pulse wave signals of the subjects in the sitting meditation control group are collected while sitting meditation;

[0009] S3, preprocessing each collected pulse wave signal;

[0010] S4, preprocessing the rhythm audio used in the experiment;

[0011] S5. extracting pulse rate variability features from the pre-processed pulse wave signal;

[0012] S6. Calculating synchronization characteristics between the pre-processed pulse wave signal and the pre-processed rhythm audio;

[0013] S7. Using the synchronization feature as the independent variable and the pulse rate variability feature as the dependent variable, a machine learning method is used to obtain a quantitative evaluation formula for cardiovascular health based on the degree of synchronization between the rhythm audio and the pulse wave signal, and a cardiovascular health prediction value is calculated based on the quantitative evaluation formula.

[0014] Furthermore, step S2 includes the following steps:

[0015] S201. The subjects were divided into three groups: a listening intervention group, a performance intervention group, and a meditation control group;

[0016] S202, using a dual-wavelength photoelectric capacitance sensor module to collect pulse wave signals, including:

[0017] The participants in the listening intervention group were required to collect pulse wave signals for 2 minutes, and then listen to different rhythmic audio clips and collect pulse wave signals again for 2 minutes.

[0018] The performance intervention group also needed to collect pulse wave signals for 2 minutes, and then play different rhythm audio clips and collect pulse wave signals again for 2 minutes;

[0019] The subjects in the meditation control group only needed to sit quietly without any music intervention, and their pulse wave signals were directly collected for 2 minutes.

[0020] Furthermore, the process of pre-processing the pulse wave signal in step S3 is as follows:

[0021] S301, performing VMD processing on the pulse wave signal;

[0022] S302, performing FFT processing on the pulse wave signal after VMD processing to obtain a spectrum diagram of each IMF component;

[0023] S303, performing wavelet threshold denoising on the high-frequency IMF component, and then performing signal reconstruction;

[0024] S304, after signal reconstruction, use FIR filter to further remove residual noise;

[0025] S305 , performing baseline removal processing on the pulse wave signal processed by the FIR filter.

[0026] Furthermore, the process of pre-processing the rhythm audio used in the experiment in step S4 is as follows:

[0027] S401, resampling the rhythm audio used in the experiment in step S2;

[0028] S402, performing bandpass filtering on the resampled rhythm audio;

[0029] S403 , performing DC offset removal and normalization processing on the rhythm audio after bandpass filtering.

[0030] Furthermore, step S5 includes the following steps:

[0031] S501, performing peak detection on the pre-processed pulse wave signal using a method based on global peak detection and local peak detection;

[0032] S502: Calculate multiple pulse rate variability features related to cardiovascular health in the pulse wave signal based on the detected peak value.

[0033] Furthermore, step S6 includes the following steps:

[0034] S601, calculating the Pearson correlation coefficient between the pulse wave signal and the rhythm audio, which is used to measure the strength and direction of the linear relationship between the pulse wave signal and the rhythm audio;

[0035] S602, calculating the mutual information between the pulse wave signal and the rhythm audio, which is used to measure the information sharing amount between the pulse wave signal and the rhythm audio;

[0036] S603, calculating the cross-spectral density between the pulse wave signal and the rhythm audio, which is used to measure the similarity between the pulse wave signal and the rhythm audio in the frequency domain;

[0037] S604: Calculate the phase lock value between the pulse wave signal and the rhythm audio signal to measure the phase synchronization between the two.

[0038] Furthermore, step S7 includes the following steps:

[0039] S701. Using a multi-factor repeated measures ANOVA method to examine the significance of changes in each pulse rate variability feature under different experimental conditions; wherein the different experimental conditions refer to: different groups or different rhythm audio used in experiments of the same group;

[0040] S702, using the synchronization feature as an independent variable and the pulse rate variability feature as a dependent variable, and using the ridge regression model to train a regression model for each pulse rate variability feature;

[0041] S703. Using different weighting methods, weight the regression models obtained from each pulse rate variability feature to obtain a quantitative evaluation formula with the cardiovascular health prediction value as the dependent variable and the synchronization feature as the independent variable, which is called the multi-feature synchronization analysis formula;

[0042] S704: Calculate the subject's cardiovascular health prediction value based on the multi-feature synchronization analysis formula.

[0043] The present invention has the following advantages and effects compared to the prior art:

[0044] Compared to traditional music intervention methods, this invention combines the synchronization of rhythmic audio and pulse wave signals to provide personalized rhythmic music intervention methods for subjects through real-time data acquisition and analysis, further improving the targetedness and effectiveness of intervention. By combining modern signal processing technology with machine learning algorithms, this invention not only enables quantitative assessment of cardiovascular health but also provides customized audio intervention plans for individuals, thereby promoting personalized health management.

[0045] This paper innovatively designs a pulse wave acquisition system based on a cloud-edge-device architecture. This system uses hardware devices to collect pulse wave signals in real time, then transmits the data to the cloud. Finally, the cloud data is downloaded on a PC for subsequent processing. Unlike traditional pulse wave acquisition methods, this cloud-edge-device architecture significantly reduces latency, alleviates cloud computing pressure, and improves the system's real-time performance and stability.

[0046] 2) This study, for the first time, combines the degree of synchronization between rhythmic audio and pulse wave signals with the PRV characteristics of the pulse wave signal to further explore the regulatory effects of rhythmic audio intervention on the cardiovascular system. This multi-dimensional, multi-feature analysis method provides a new perspective for understanding the effects of rhythmic music on the human body and provides a theoretical basis for the subsequent development of personalized audio intervention programs.

[0047] 3) The present invention proposes a multi-feature synchronization analysis method, which jointly analyzes the changing trends of synchronization features and PRV features, and weights the synchronization features through machine learning and significance testing methods, so as to quantify the effects of different rhythmic audio and different intervention methods on cardiovascular health. Finally, based on this, personalized rhythmic audio intervention methods are recommended to the subjects. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is a flow chart of a method for quantitatively evaluating the cardiovascular effects of rhythmic audio pulse wave synchronization in an embodiment of the present invention;

[0049] Figure 2 This is the architecture diagram of the cloud-edge-based pulse wave collection system;

[0050] Figure 3 Schematic diagram of endpoint detection based on short-time energy zero ratio;

[0051] Figure 4 This is the experimental flow chart;

[0052] Figure 5 This is the decomposition result graph of the VMD algorithm;

[0053] Figure 6 is the spectrum diagram of each IMF;

[0054] Figure 7 is a graph of the audio signal during preprocessing;

[0055] Figure 8 A visualization diagram of the synchronization features;

[0056] Figure 9 The F-value histogram of variance analysis of each PRV feature. DETAILED DESCRIPTION

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0058] A quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization, the detection method comprising the following steps:

[0059] S1. After performing endpoint detection on single audio of three Chinese drum tube sounds with different intensities, multiple sets of rhythm audio are produced based on music theory.

[0060] S2. The experiment was conducted after the subjects were divided into three groups: a listening intervention group, a playing intervention group, and a meditation control group. The listening intervention group required the subjects to listen to rhythmic audio, and their pulse waves were collected before and after the experiment; the playing intervention group required the subjects to play rhythmic audio, and their pulse waves were collected before and after the experiment; the meditation control group only required the subjects' pulse waves to be collected while they were meditating.

[0061] S3. Preprocess the pulse wave signal collected in each experiment;

[0062] S4. Preprocess the rhythm audio used in the experiment;

[0063] S5. performing feature extraction on the preprocessed pulse wave signal, wherein the extracted feature is pulse rate variability (PRV) feature related to cardiovascular health;

[0064] S6. Calculating a synchronization feature between the pre-processed pulse wave signal and the rhythm audio (pre-processed) used in the experiment for obtaining the pulse wave signal;

[0065] S7. Using synchronization features as independent variables and PRV features as dependent variables, and utilizing machine learning methods, we obtain a quantitative evaluation formula for cardiovascular health based on the degree of synchronization between rhythmic audio and pulse waves, thereby providing personalized recommendations for rhythmic audio intervention methods for experimenters.

[0066] The processing method in step S1 is as follows:

[0067] S101. Perform noise suppression on the three single-tone audios required for the experiment. The three single-tone audios are the Chinese drum drum sound - strong dynamic audio, drum drum sound - medium dynamic audio, and drum drum sound - weak dynamic audio. The core of noise suppression is to set a noise threshold value and suppress the part of the signal that is less than the threshold value (usually set to zero). This is a simple and effective noise reduction method. The following is the specific process and implementation:

[0068]

[0069] Among them, x is the sample value of the signal, and NoiseLevel is the set noise threshold value;

[0070] S102, using the short-time energy zero ratio method to perform endpoint detection on the single-tone audio. The short-time energy zero ratio is the ratio of the short-time energy to the short-time zero-crossing rate, and the short-time energy zero ratio EZR i The calculation formula is as follows:

[0071]

[0072] Where b is a constant, LE i is the short-time energy, and the calculation formula is:

[0073]

[0074] ZCR n is the short-time average zero-crossing rate, and the calculation formula is:

[0075]

[0076] S103. Apply weighted sliding average smoothing to the short-term energy and zero-crossing rate. The purpose is to reduce fluctuations and improve the accuracy of signal trend analysis. The calculation formula is as follows:

[0077]

[0078] S104: Use baseline noise to dynamically calculate energy and short-time energy zero ratio threshold to achieve endpoint detection of audio. This is used to distinguish valid signals from background noise. The calculation formula for baseline noise is:

[0079]

[0080] The calculation formula of dynamic energy threshold is as follows:

[0081] Energy thr =BaselineNoise+k·max(LE i )

[0082] Where k is the scaling factor, usually set to k = 0.2. max(LE i ) is the maximum value of all short-time energies in the signal.

[0083] The calculation formula for the dynamic EZR threshold is as follows:

[0084] EZR thr =c·max(EZR)

[0085] Where c is the scaling factor, usually set to c = 0.1. max(EZR) is the maximum value of all EZRs in the signal. By determining the short-time energy and the short-time energy zero ratio, the starting point and end point of the signal are found. The calculation formula is:

[0086] StartFrame=min{t|LE i [t]>Energy thr ∧EZR[t]>EZR thr}

[0087] EndFrame=max{t|LE i [t]>Energy thr ∧EZR[t]>EZR thr}

[0088] Among them, StartFrame and EndFrame are the frame indexes of the starting point and end point.

[0089] StartTime=StartFrame·FrameStride

[0090] EndTime=EndFrame·FrameStride

[0091] Among them, FrameStride is the frame shift, that is, the time interval between each frame, StartTime and EndTime are the starting point and end point of the signal after endpoint detection;

[0092] S105. According to music theory, the monophonic audios after endpoint detection are combined to create 2 / 4, 3 / 4, and 4 / 4 beat rhythm audios with beats per minute (BPM) of 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, and 150, respectively. 2 / 4 beat means that a quarter note is used as a beat, each measure has 2 beats, and each measure is composed of an audio sequence of "strong dynamics-weak dynamics"; 3 / 4 beat means that a quarter note is used as a beat, each measure has 3 beats, and each measure is composed of an audio sequence of "strong dynamics-weak dynamics-weak dynamics"; 4 / 4 beat means that a quarter note is used as a beat, each measure has 4 beats, and each measure is composed of an audio sequence of "strong dynamics-weak dynamics-less strong dynamics"

[0093] Furthermore, the experiment process after dividing the subjects into three groups in step S2: a listening intervention group, a performance intervention group, and a meditation control group is as follows:

[0094] S201. Divide the subjects into three groups: a listening intervention group, a performance intervention group, and a meditation control group. Subjects in the listening intervention group were required to first collect pulse wave signals for 2 minutes, then listen to different rhythmic audio clips, and then collect pulse wave signals for another 2 minutes; subjects in the performance intervention group were also required to first collect pulse wave signals for 2 minutes, then play different rhythmic audio clips, and then collect pulse wave signals for another 2 minutes; subjects in the meditation group only needed to sit quietly, without music intervention, and their pulse wave signals were directly collected for 2 minutes;

[0095] S202: Use a dual-wavelength photoelectric capacitance sensor module to collect pulse wave signals.

[0096] Furthermore, the process of pre-processing the pulse wave signal in step S3 is as follows:

[0097] S301. Perform variational mode decomposition (VMD) on the pulse wave signal. VMD decomposes the complex signal into several intrinsic mode functions (IMFs) and a residual term. The specific steps are as follows:

[0098] Extract all local maxima and local minima of the original pulse wave signal x(t), connect all maxima and minima respectively using cubic spline interpolation to generate the upper envelope e max (t) and the lower envelope e min (t):

[0099] Calculate the average value m1(t) of the upper and lower envelopes:

[0100]

[0101] Extract IMF candidate classification from the original pulse wave signal x(t)

[0102]

[0103] generally It is not a stationary signal and does not meet the two conditions defined by IMF. Repeat the above process. Assume that after k times (k is generally less than 10) If the definition of IMF is satisfied, the first-order IMF component IMF1(t) of the original signal x(t) is:

[0104]

[0105] Subtract c1(t) from the original signal x(t) to get a new signal r1(t) without high-frequency components:

[0106] r1(t)=x(t)-c1(t)

[0107] Repeat the process of obtaining c1(t) for r1(t) to obtain the second IMF component c2(t), and repeat this process until the nth order IMF component c n (t) or the remainder r n (t) is less than the preset value; or when the residual component r n When (t) is a monotonic function or a constant, the VMD decomposition process stops. Finally, the original signal can be expressed as follows after VMD decomposition:

[0108]

[0109] Among them, the IMF i (t) is the i-th order IMF component;

[0110] S302, perform FFT processing on the pulse wave signal after VMD processing to obtain the spectrum of each IMF component:

[0111] The peak frequency f of the spectrum of each IMF component p Make a judgment, if f p >f threshold , then the IMF component is divided into high-frequency IMF components and enters the wavelet threshold denoising process; if f p ≤f threshold , then the IMF component is retained without wavelet threshold denoising, and is subsequently directly reconstructed with the high-frequency IMF component after wavelet denoising to form a new pulse wave signal;

[0112] S303: Perform wavelet threshold denoising on the high-frequency IMF components, and then reconstruct the signal. Wavelet threshold denoising is a signal processing technique based on wavelet transform. By decomposing the signal into sub-bands of different frequencies, it is possible to effectively distinguish between signal and noise. First, perform wavelet decomposition to decompose the high-frequency IMF components into approximation coefficients cA and detail coefficients cD:

[0113]

[0114] This is followed by thresholding. This method uses a hard thresholding method for denoising because the pulse wave signal has a low amplitude, and hard thresholding can better preserve important details such as heart rate fluctuations. A global threshold is used to balance denoising effectiveness with signal fidelity. The global threshold α is calculated by the median σ of the highest-frequency component (i.e., cD1) in the detail coefficients, based on the assumption that the high-frequency component is primarily dominated by noise:

[0115] σ=median(|cD1|)

[0116]

[0117] Where N is the signal length. The calculated α can remove small coefficients caused by high-frequency noise in the detail coefficients while retaining important signal components. The hard threshold formula is:

[0118]

[0119] For each layer detail coefficient cD i Apply hard thresholding layer by layer (starting from i=1). Set the coefficients whose absolute value is less than the threshold α to 0, and keep the other coefficients.

[0120] After processing the detail coefficients, the high-frequency IMF signal is reconstructed. The processed coefficients are combined into a denoised high-frequency IMF signal through inverse wavelet transform. At the same time, in order to ensure that the signal length is consistent with the original signal, the result needs to be cropped. Each denoised high-frequency IMF signal is expressed as:

[0121] IMF i (t)=Inverse(cA n ,{cD1,cD2,...,cD n})

[0122] After VMD and wavelet denoising, the low-frequency and denoised high-frequency IMF components are reconstructed again to obtain the denoised signal. The reconstruction formula is as follows:

[0123]

[0124] IMF i(t) is each IMF component after denoising;

[0125] S304: After signal reconstruction, use an FIR filter to further remove residual noise. The FIR filter used in this method is designed as a window function method, which is to obtain a practical filter by designing the impulse response of an ideal filter and weighting it with a window function.

[0126] First, for a low-pass filter, the ideal frequency response H d (f) is:

[0127]

[0128] where f c is the cutoff frequency of the filter. The impulse response h of the ideal low-pass filter is d [n] is the ideal frequency response H d The inverse Fourier transform of (f) is:

[0129]

[0130] This is a signal with sinc function characteristics, which has an infinite time response. In order to change the ideal impulse response into a finite length, this method uses the window function ω[n] for windowing processing.

[0131] h[n]=h d [n]·ω[n]

[0132] The impulse response h[n] after weighting by the window function is the coefficient of the FIR filter. FIR [n] is reconstructed from the filter coefficients h[n] and the signal x 重构 (t) is obtained by convolution:

[0133]

[0134] S305: Perform baseline removal on the pulse wave signal after FIR filtering. The purpose of baseline removal is to remove baseline drift from the signal, preserve true heartbeat fluctuations and other physiological changes, and improve the accuracy of signal analysis. This baseline removal process consists of three steps: reference point detection, baseline fitting, and baseline removal.

[0135] First, benchmark detection is divided into peak detection and starting point detection. Peak point detection uses a sliding window to detect the peak of the signal. Specifically, for each point in the signal, the values ​​in the left and right windows are checked. If the value of the point is greater than the values ​​of the left and right points, it is considered a peak.

[0136] Starting point detection uses the results of peak detection to further determine the reference point (i.e. the starting point of the heartbeat fluctuation) by calculating the difference of the signal. First, the signal difference is processed to detect the direction of signal change and find the starting point of the upward trend.

[0137] Then comes the cubic spline interpolation, which is the core step of baseline removal. The baseline of the signal is fitted by cubic spline interpolation and then removed from the original signal. Cubic spline interpolation is a method commonly used for smoothing signals. It generates a smooth curve to approximate the trend of the signal through continuous second-order derivatives and higher-order fitting accuracy. For a set of discrete data points (x1, y1), (x2, y2), ..., (x n ,y n The goal of cubic spline interpolation is to find a function S(x) that passes through the interpolation point at every data point and remains smooth between those points. Cubic spline interpolation ensures that: the interpolated curve passes through every data point; and the derivative of the curve is continuous. Using cubic spline interpolation, we can obtain the fitted baseline f(t).

[0138] Finally, the signal x after FIR filtering is FIR Subtract the fitted baseline f(t) from [n] to obtain the signal x after removing the baseline drift 去基线 (t). At the same time, in order to maintain the amplitude of the signal, the signal after removing the baseline is added to the mean of the original signal.

[0139]

[0140] Furthermore, the process of pre-processing the rhythm audio used in the experiment in step S4 is as follows:

[0141] S401. Resample the rhythm audio used in the experiment in S2. This method uses linear interpolation for resampling. When calculating each new sampling point, the new sample value is estimated from the original sample through a simple linear relationship. The specific interpolation steps are as follows: Given two data points (x1, y1) and (x2, y2), then for any given x value (between x1 and x2), the result y of linear interpolation is:

[0142]

[0143] Where x1≤x≤x2, and x is the new sampling point.

[0144] And for the new sampling point time position t of the resampled signal new It is calculated according to the following formula:

[0145]

[0146] Among them, f s is the target sampling frequency, f original is the sampling frequency of the original signal, and t is the time position corresponding to a sample of the original signal;

[0147] S402: Bandpass filter the resampled rhythm audio to remove irrelevant frequency components to prevent them from interfering with subsequent analysis and calculations.

[0148] S403: Perform DC offset removal and normalization on the bandpass filtered audio. First, removing the DC component facilitates subsequent signal analysis by preventing it from affecting signal periodicity. Second, normalization ensures that the amplitudes of different signals are within the same range, facilitating comparison and analysis.

[0149] Furthermore, the process of extracting pulse rate variability (PRV) features from the pre-processed pulse wave signal in step S5 is as follows:

[0150] S501. Perform peak detection using a method based on global peak detection and local peak detection. Since the peak point detection method in S305 may miss detection and make false detections at the front end of the signal, this method proposes a peak detection method based on global peak detection and local peak detection. First, calculate the mean and standard deviation of the signal and set a global threshold:

[0151] global_threshold=μ(x)+k·σ(x)

[0152] Where μ(x) is the mean of the signal, σ(x) is the standard deviation of the signal, and k is the amplification factor. This threshold is used to detect the global peak.

[0153] Calculate the local threshold using the mean and standard deviation of the local signal within the first 2 seconds of the signal:

[0154] local_threshold=μ(x local )+k·σ(x local )

[0155] where μ(x local ) and σ(x local ) are the mean and standard deviation of the local signal, respectively, to detect local peaks. By merging local and global peaks, we obtain a final peak set and remove duplicate peak indices. This method effectively identifies peaks in the signal and avoids false detections due to noise.

[0156] The starting point is then determined by detecting the minimum value between each pair of adjacent peaks. The starting point is typically located at the lowest point within each cycle, representing the beginning of the signal waveform. For each pair of adjacent peaks, the signal segment between them is extracted, and the minimum point in that segment is found as the starting point.

[0157] S502: Calculate multiple PRV features related to cardiovascular health in the pulse wave signal. The calculation process is based on the RR interval, which is the time difference between adjacent peaks. The RR interval is calculated by the difference between the indices of two adjacent peak points.

[0158] The first thing to calculate is PNN50, which refers to the proportion of adjacent RR intervals with a difference of more than 50ms. The calculation formula is as follows:

[0159]

[0160] The RangeNN is then calculated, which represents the difference between the maximum and minimum values ​​of the RR interval. The calculation formula is as follows:

[0161] RangeNN = max(RR) - min(RR)

[0162] Then we calculate SDNN, which represents the standard deviation of the RR interval:

[0163] SDNN=std(RR)

[0164] Then the frequency domain characteristics of the pulse wave signal are calculated. The present invention uses the Lomb-Scargle (LS) method for power spectrum analysis. The LS method is suitable for short-time signal analysis. By calculating the power of the signal at different frequencies, it reveals the characteristics of the signal in the frequency domain. The Lomb-Scargle method can effectively extract the periodic components in the signal. By analyzing the spectrum of the signal, the Lomb-Scargle method can reveal the power distribution of low and high frequencies in the signal, and then analyze the relative intensity of the LF (low frequency) and HF (high frequency) regions. Lomb-Scargle power spectrum analysis can quantify the main frequency components in the signal by analyzing the power spectrum, revealing the energy distribution of the signal in a specific frequency band. The calculation formula is as follows:

[0165]

[0166] where t i is the sampling time of the signal, y i is the amplitude of the signal, and f is the frequency.

[0167] By calculating the power of the signal in the low frequency and high frequency range, the LF, HF and The frequency range of low-frequency power is 0.04-0.15Hz, and the frequency range of high-frequency power is 0.15-0.4Hz. The calculation formula is as follows:

[0168]

[0169]

[0170] Then we can find the ratio of the two and get ratio.

[0171] Finally, the nonlinear characteristics of the pulse wave signal are calculated. The present invention selects the SD1 parameter and the sample entropy (SampEn).

[0172] First, let's introduce the Poincaré plot. The Poincaré plot is a nonlinear method for analyzing time series data, particularly widely used in the study of heart rate variability (HRV) and pulse rate variability (PRV). It graphically displays the changing pattern of heartbeat intervals by comparing the value of each RR interval (interbeat interval) with the value of the previous RR interval. The plot uses RR intervals (the time interval between consecutive heartbeats) as data points, typically pairing each RR interval with its previous RR interval. The position of each point is determined by two parameters: the horizontal axis represents the current RR interval, and the vertical axis represents the previous RR interval. By plotting all data points on the graph, a scatter plot is formed. This graph illustrates the changing patterns of heartbeat intervals, helping us understand the regularity and complexity of cardiac rhythms. In the Poincaré plot, the SD1 parameter is the standard deviation of the points in the plot, indicating the short-term variability along the unit line. It reflects short-term fluctuations in heartbeat intervals. The calculation formula is as follows:

[0173]

[0174] Where ΔRR is the difference sequence between adjacent RR intervals, and Var(ΔRR) is the variance of the sequence.

[0175] For a time series of length N {x1,x2,...,x n}, the embedding dimension is m, and the tolerance is r. Set r to 0.2 times the standard deviation of the signal. First, extract all continuous subsequences of length m from the time series to form a set of template vectors. For each pair of template vectors, calculate the similarity between them. This method uses the maximum difference (Chebyshev distance) metric. Calculate the number of template vector pairs that satisfy a similarity less than r, denoted as B(m,r). Then, construct a template vector of length m+1, calculate the similarity, and count the number of template vector pairs that satisfy a similarity less than r, denoted as A(m+1,r). The calculation formula of SampEn is as follows:

[0176]

[0177] Furthermore, the process of calculating the synchronization feature between the pre-processed pulse wave signal and the rhythm audio (pre-processed) used in the experiment to obtain the pulse wave signal in step S6 is as follows:

[0178] S601. Calculate the Pearson correlation coefficient (PCCs). This measures the strength and direction of the linear relationship between the pulse wave signal and the rhythm audio signal. The value of the correlation coefficient ranges from -1 to 1, where:

[0179] 1 indicates a perfect positive correlation (i.e., when one variable increases, the other variable also increases).

[0180] -1 indicates a perfect negative correlation (i.e., when one variable increases, the other decreases completely).

[0181] 0 means there is no linear relationship between the two variables.

[0182] The calculation formula of Pearson correlation coefficient is as follows:

[0183]

[0184] Among them, r xy is the Pearson correlation coefficient, x i and y i Represents the i-th data point of two signals x and y respectively, and Represents the sample means of the two signals x and y respectively, and n is the total number of data points;

[0185] S602. Calculate mutual information (MI). Mutual information measures the amount of information shared between the pulse wave signal and the rhythm audio signal and is applicable to relationships between nonlinear and non-Gaussian signals. In this method, mutual information is used to evaluate the dependency between two signals, and is calculated as follows:

[0186] I(X,Y)=H(X)+H(Y)-H(X,Y)

[0187] Where H(X) and H(Y) are the entropies of signals X and Y, and H(X,Y) is their joint entropy. In this method, the signals are first converted to discrete values, the appropriate bin boundaries are calculated, the signals are converted to discrete integer representations, and then the mutual information between them is calculated based on the discretized signals to measure the dependence of the two signals.

[0188] S603. Calculate the cross spectral density (CSD). The cross spectral density measures the similarity between the pulse wave signal and the rhythm audio signal in the frequency domain. It reveals the synergistic relationship between the two signals at different frequencies, especially whether they have synchronous changes or frequency dependence at specific frequencies. CSD is a complex value, and its modulus (amplitude) represents the strength of the correlation between the signals. Frequencies with large amplitudes indicate that there is a stronger relationship between the signals in this frequency band. First, perform Fourier transform on the signals x(t) and y(t) to convert them from the time domain to the frequency domain. Its calculation formula is:

[0189]

[0190] Next, calculate the complex conjugate Y of Y(f) * (f). The definition of complex conjugate is:

[0191]

[0192] In the frequency domain, complex conjugation reverses the sign of the imaginary part of a complex number.

[0193] Then, calculate the expected value E[X(f)Y * (f)], which is the mutual covariance of signals X(f) and Y(f) in the frequency domain. For a signal of finite duration, the expected value can be approximated by estimating the mean. Typically, the expected value is calculated from data in multiple windows.

[0194] Finally, by calculating S for each frequency component XY (f), we can get the cross spectral density of the whole signal:

[0195]

[0196] Among them, S XY (f) is the cross spectral density, N is the number of sliding windows or segments, X n (f) and Y n (f) is the frequency domain representation of the nth segment signal, Y n * (f) is the complex conjugate. This method calculates the cross-spectral density value of the frequency value corresponding to the BPM of the rhythmic audio;

[0197] S604. Calculate the Phase Lock Value (PLV). This method measures the phase synchronization between the pulse wave signal and the rhythm audio signal by calculating the instantaneous phase and the average of the phase difference between them. In this method, the instantaneous phase of each signal is extracted using a Hilbert transform. The phase difference between the two signals is then calculated, and the phase consistency of the two signals is measured using the PLV formula.

[0198] Specifically, for two signals x1(t) and x2(t), their analytical signals are first obtained by Hilbert transform;

[0199] Then, the instantaneous phases θ1(t) and θ2(t) are extracted:

[0200] θ1(t)=arg(z1(t)), θ2(t)=arg(z2(t))

[0201] Finally, calculate the phase lock value PLV:

[0202]

[0203] Four synchronization features are obtained: Pearson correlation coefficient, mutual information, phase locking value and cross spectral density.

[0204] Furthermore, in step S7, using the synchronization feature as the independent variable and the PRV feature as the dependent variable, a machine learning method is used to obtain a quantitative evaluation formula for the degree of synchronization between rhythm audio and pulse wave for cardiovascular health, thereby providing the subject with personalized rhythm audio intervention method recommendations as follows:

[0205] S701. Use a multi-factor repeated measures ANOVA method to examine the significance of changes in each PRV feature under different experimental conditions. The different experimental conditions refer to different groups or different rhythm audio used in experiments within the same group. Screen out PRV features that change significantly under different experimental conditions, proving that the current experimental conditions have a greater impact on these PRV features related to cardiovascular health. PRV features that do not change significantly indicate that the current experimental conditions have little impact on them and can be screened out.

[0206] Before starting the analysis, it is necessary to perform normality tests and homogeneity tests on the changes in the PRV eigenvalues ​​under each experimental condition to ensure that the data meet the conditions for using repeated measures ANOVA. The two assumptions made in multi-factor repeated measures ANOVA are:

[0207] Null hypothesis (H0): There is no significant difference in the mean of a certain PRV feature under different experimental conditions.

[0208] Alternative hypothesis (H1): There are significant differences in the mean of a certain PRV feature under different experimental conditions.

[0209] Hypothesis testing is performed by calculating the F value.

[0210] Consult the F distribution table. If the calculated F value is greater than the critical F value, the null hypothesis can be rejected, and the change in PRV characteristics under the experimental conditions is significant. Otherwise, the null hypothesis is accepted, and the change in PRV characteristics under the experimental conditions is not significant. For PRV characteristics that do not change significantly, it proves that the current experimental conditions have little effect on them, so they are screened out.

[0211] S702 , using the synchronization feature as an independent variable and the PRV feature as a dependent variable, and using the ridge regression model to train a regression model for each PRV feature.

[0212] First, the PRV features need to be aligned for directional consistency. Literature indicates that LF and LF / HF within the PRV feature are inversely proportional to cardiovascular health; that is, lower LF and LF / HF values ​​indicate better cardiovascular health. Since the other PRV features studied in this method are all directly proportional to cardiovascular health, the resulting LF and LF / HF values ​​need to be inverted before further processing. This adjusts the PRV features for consistent direction.

[0213] Then, in order to eliminate individual differences, before building the model, the four synchronization eigenvalues ​​calculated by the experimenter after the experiment were subtracted from the four synchronization eigenvalues ​​calculated in the sitting state, and the various PRV eigenvalues ​​obtained after the experiment were subtracted from the various PRV eigenvalues ​​in the sitting state to complete the data de-individualization processing.

[0214] The dedifferentiated values ​​are then normalized using Min-Max normalization, converting the four synchronized eigenvalues ​​and each PRV eigenvalue to the range [0, 1]. Using normalized values ​​for data processing not only reflects the degree of change but also eliminates the effects of different dimensionality of the data.

[0215] This method uses the ridge regression algorithm to build a model. The normalized synchronized feature value is used as the independent variable, and the normalized PRV feature value is used as the dependent variable. A regression model is trained for each PRV feature. Each PRV feature corresponds to a regression formula. In the formula, a weight is assigned to each synchronized feature, indicating the contribution of the synchronized feature to the change of the PRV feature. Hyperparameter tuning and cross-validation are performed through grid search to evaluate the performance of the model and ensure the stability and generalization ability of the model. The multi-feature synchronized regression formula corresponding to each PRV feature follows the following form:

[0216]

[0217] in, is the dependent variable, representing the normalized predicted value of a PRV feature. x is a column vector of independent variables, representing the normalized values ​​of the synchronized features. w is a column vector of regression coefficients, representing the contribution of each synchronized feature to the change in the PRV feature. b is the bias constant, representing the bias of the model.

[0218] Ridge Regression is a regression method with L2 regularization, which is mainly used to deal with multicollinearity problems and prevent overfitting. When using ridge regression, we not only want to minimize the residual sum of squares, but also want to introduce a regularization term to penalize the complexity of the model (that is, the size of the regression coefficient), so that the regression coefficient is smoother and avoids overfitting. Ridge regression obtains the optimal regression coefficient by minimizing the objective function, and its objective function is:

[0219] J(w)=‖y-Xw|| 2 +λ||w|| 2

[0220] Among them, the dependent variable y is an n×1 column vector, n is the number of subjects, representing the normalized value (true value) of the PRV feature of all subjects. The independent variable X is an n×p matrix, p is the number of synchronized features, each column is the normalized value of a different synchronized feature, and each row is the normalized value of the synchronized feature of each subject under different experimental conditions. w is the p×1 regression coefficient column vector to be learned, which represents the weight of each synchronized feature and the contribution of each synchronized feature to the amplitude of the change in the PRV feature. λ is the regularization parameter, which controls the complexity of the model and avoids overfitting. The computational goal of ridge regression is to find the regression coefficient w so that the objective function J(w) is minimized:

[0221]

[0222] Among them, ||y-Xw|| and ||w|| refer to the norms of the two vectors. In the objective function, it represents the L2 norm, also called the Euclidean norm.

[0223] The regression coefficient w can be calculated using the following formula:

[0224] w=(X T X+λI) -1 X T y

[0225] Among them, X T X is the product of the transpose of the independent variable matrix and itself, which represents the relationship between the synchronized features. λI is the regularization term, where I is the identity matrix and λ controls the regularization strength. T y is the inner product of the transpose of the independent variable matrix and the dependent variable (PRV feature change).

[0226] To ensure the stability and generalization ability of each model, this method uses GridSearch for hyperparameter tuning to find the most appropriate regularization parameter λ. At the same time, k-fold cross-validation is used to evaluate the performance of each model to ensure that the model is not overfitting and has good generalization ability.

[0227] S703, using equal weight, standardized effect size Cohen's d value, standardized model explanatory power R 2 , standardized mean square error MSE and comprehensive consideration of Cohen's d, R 2 The regression formula obtained by each PRV feature is weighted using five different weighting methods, including PRV and MSE, to obtain a regression formula with cardiovascular health as the dependent variable; the regression formula with four synchronized features as independent variables is called the multi-feature synchronized analysis formula.

[0228] Since no relevant research has objectively quantified the correlation between PRV characteristics and cardiovascular health, this method attempts to weight each PRV characteristic regression formula using the five different weighting methods mentioned above, thereby deriving a multi-characteristic synchronization analysis formula. The prediction error of the integrated model under the five different weighting methods is calculated as a performance indicator to measure the performance of the five weighting methods.

[0229] First, equal weight is to give the same weight to each PRV feature regression formula, the weight is q represents the number of PRV features with significant changes after screening.

[0230] We then attempted to use the standardized effect size, Cohen's d, to weight the regression formula. Calculating Cohen's d quantitatively measures the magnitude of change in PRV characteristics before and after the experiment. Cohen's d is a measure of the effect size between two sets of data and is well-suited for comparing the magnitude of change in PRV characteristics before and after an experiment. By calculating the Cohen's d value for each PRV characteristic before and after each experiment, we can determine the magnitude of change in a particular PRV characteristic under a given experimental condition.

[0231] The obtained Cohen's d is normalized. This method performs the following normalization: the sum of the Cohen's d of the same PRV feature calculated under all experimental conditions is divided by the sum of the Cohen's d of all PRV features calculated under all experimental conditions. This normalization method can simultaneously consider the changes in the same PRV feature under different experimental conditions and the overall changes in all PRV features under all conditions. It can comprehensively consider the effect size of each experimental condition and PRV feature, thereby achieving a more comprehensive normalization. The calculation formula is as follows:

[0232]

[0233] Among them, β i and d i ' represents the Cohen's d value after the normalization of the i-th PRV feature, that is, the weight of each regression formula, d i,j represents the Cohen's d value calculated for a PRV feature under certain experimental conditions, m represents the total number of experimental conditions, and q represents the number of PRV features. Therefore, the corresponding synchronization feature weight calculation formula after the regression formula is weighted is:

[0234]

[0235] Similarly, we try to use the standardized model explanatory power R 2 Perform weighting of the regression formula. Each PRV feature regression model obtained above corresponds to a model explanatory power R 2 , for R 2 Perform the following standardization processing, and the calculation formula is as follows:

[0236]

[0237] Among them, β i and R i 2 ' represents the model explanatory power of the standardized PRV feature, that is, the weight of each regression formula, and q represents the number of PRV features. The weight of the kth synchronized feature corresponding to the final regression formula is:

[0238]

[0239] Similarly, we try to use the standardized mean square error (MSE) to weight the regression formula. Each PRV feature regression model obtained above corresponds to a mean square error (MSE), which is consistent with Cohen's d and R 2The difference is that the smaller the mean square error (MSE), the smaller the prediction error, and the greater the weight should be given. Therefore, when performing normalization, this method first takes the inverse of the MSE and then performs the following normalization. The calculation formula is as follows:

[0240]

[0241] Among them, β i and (1 / MSE i )' represents the mean square error of the normalized PRV feature i, i.e., the weight of each regression formula, and q represents the number of PRV features. The k-th synchronized feature weight corresponding to the final regression formula is:

[0242]

[0243] Finally, we tried to use the comprehensive standardized Cohen's d value, standardized R 2 The regression formula is weighted according to the standardized Cohen's d value and the standardized R 2 and the standardized 1 / MSE, Cohen's d value, R 2 and 1 / MSE to assign weight proportions, and assign weight proportions a to Cohen's d value, R 2 Assign a weight ratio of b, and 1 / MSE a weight ratio of c, satisfying a+b+c=1. Then the weight corresponding to the i-th PRV regression formula is:

[0244] β i =a·d′ i +b·R i 2 ′+c·(1 / MSE i )′

[0245] The kth synchronized feature weight corresponding to the final regression formula is:

[0246]

[0247] The synchronized feature regression coefficient column vector w'=[w1',w2',···,w p '] T . p is the number of synchronized features.

[0248] The final multi-feature synchronized regression formula based on cardiovascular health is as follows:

[0249]

[0250] Among them, the dependent variable The cardiovascular health prediction value obtained by weighting the normalized PRV feature. The independent variable x is the normalized value of the synchronization feature. b' is the bias constant, b i is the bias constant in each regression formula.

[0251] The predicted value of the regression formula corresponding to the i-th PRV feature defined by this method is So the dependent variable It can also be written as:

[0252]

[0253] The normalized true value of the i-th PRV feature is y i , so the true value Y of cardiovascular health is:

[0254]

[0255] By calculating the prediction error between the predicted value and the true value under 5 different weighting methods, the optimal weighting method is selected;

[0256] S704: Based on the obtained multi-feature synchronization analysis formula, four synchronization features are calculated between the subject's pulse wave collected after each experiment and the corresponding rhythmic audio. These features are then substituted into the multi-feature synchronization analysis formula to obtain a cardiovascular health prediction value for that experimental condition. The experimental condition that yields the maximum cardiovascular health prediction value is recommended to the subject, indicating that this experimental condition is likely to have the most positive impact on the subject's cardiovascular health, thereby achieving personalized recommendations for rhythmic audio intervention.

[0257] Example 1

[0258] Figure 1 Flowchart of the quantitative evaluation method of the cardiovascular effect of rhythmic audio pulse wave synchronization in an embodiment of the present invention. Figure 1 The heart rate variability feature in is the pulse rate variability feature.

[0259] Figure 2 This is the architecture diagram of the pulse wave acquisition system based on cloud-edge hardware devices;

[0260] S1. Obtain the three Chinese drum tube sound monophonic audios of different intensities required for the experiment from the sound library, perform noise suppression and endpoint detection on them, and synthesize them into the rhythm audio required for the experiment based on music theory knowledge.

[0261] Noise suppression requires setting the noise threshold NoiseLevel. This method sets NoiseLevel to 2% of the maximum signal amplitude.

[0262] Figure 3Schematic diagram of endpoint detection based on short-time energy ratio, where the high threshold is selected as 0.5EZR m , the lower threshold is selected as 0.3EZR m , EZR m It is the average value after median filtering of the short-time energy zero ratio;

[0263] S2, Figure 4 The experimental flow chart for obtaining pulse wave signals of three groups of subjects: the listening intervention group, the playing intervention group, and the meditation control group;

[0264] S3. Preprocess the pulse wave signal obtained in step S2:

[0265] S301, such as Figure 5 The VMD (Variational Mode Decomposition) decomposition diagram shown in the figure is used to process the obtained pulse wave signal. This method decomposes the pulse wave signal into 10 IMF components.

[0266] S302, such as Figure 6 The spectrum diagram of each IMF component is shown, and the spectrum peak frequency f of each IMF component is p Make a judgment, if f p >f threshold , then the IMF component is divided into high-frequency IMF components and enters the wavelet threshold denoising process; if f p ≤f threshold , then the IMF component is retained and wavelet threshold denoising is not required. According to the frequency characteristics of the pulse wave signal, the effective signal is mainly concentrated in 0-3Hz, so this method will f threshold Set to 3Hz.

[0267] S303, for f p The IMF components exceeding 3 Hz are defined as high-frequency components and subjected to wavelet threshold denoising.

[0268] The high-frequency IMF component is decomposed into the approximate coefficient cA and the detail coefficient cD. This method uses the db8 wavelet basis and decomposes the signal into 5 layers (i.e., cA5 and cD1~cD5). Among them, cA5 represents the highest layer approximate coefficient, which represents the low-frequency component of the signal; cD1~cD5 represent the detail coefficients of each layer, which represent the high-frequency detail component of the signal:

[0269]

[0270] Then threshold processing is performed. This method uses the hard threshold method for denoising. The global threshold α is calculated by the median σ of the highest frequency part (i.e. cD1) in the detail coefficient, which is based on the assumption that the high frequency part is mainly dominated by noise:

[0271] σ=median(|cD1|)

[0272]

[0273] Where N is the signal length. The calculated α can remove small coefficients caused by high-frequency noise in the detail coefficients while retaining important signal components. The hard threshold formula is:

[0274]

[0275] For each layer detail coefficient cD i Apply hard thresholding layer by layer (starting from i=1). Set the coefficients whose absolute value is less than the threshold α to 0, and keep the other coefficients.

[0276] After processing the detail coefficients, the high-frequency IMF signal is reconstructed. The processed coefficients are combined into a denoised high-frequency IMF signal through inverse wavelet transform. At the same time, in order to ensure that the signal length is consistent with the original signal, the result needs to be cropped. Each denoised high-frequency IMF signal is expressed as:

[0277] IMF i (t)=Inverse(cA5,{cD1,cD2,...,cD5})

[0278] After VMD and wavelet denoising, the low-frequency and denoised high-frequency IMF components are reconstructed again to obtain the denoised signal. The reconstruction formula is as follows:

[0279]

[0280] IMF i (t) is each IMF component after denoising.

[0281] S304. The reconstructed signal is further subjected to an FIR filter to remove residual noise and smooth the signal. The FIR filter used in this method is designed as a window function method, and the window function used is a Hamming window. The Hamming window can significantly reduce the width of the main lobe and has low side lobes, which can effectively reduce the interference of the filter on the signal. It is suitable for low-pass filter design, especially in situations where stability and moderate computational complexity are required. It is very suitable for processing pulse wave signals and is defined as:

[0282]

[0283] Where M is the order of the filter, which determines the length of the window.

[0284] The hardware sampling frequency of the pulse wave signal is 100Hz, so the sampling frequency of the FIR filter is set to 100Hz; the cutoff frequency is based on the frequency range of the pulse wave signal, and 3Hz is selected as a reasonable cutoff frequency, which can effectively remove noise above 3Hz while retaining the useful components of the signal; the filter order is selected as 64, which can achieve better filtering effects and avoid excessive computational complexity; the boundary rejection parameter is selected as 100, which means removing the first 100 data points and the last 100 data points generated during the filter calculation. Usually, this part contains a large boundary effect. Removing these data points can ensure that the filtered signal is more stable and accurate. The impulse response h[n] after weighting by the window function is the coefficient of the FIR filter. Signal x FIR [n] is reconstructed from the filter coefficients h[n] and the signal x 重构 (t) is obtained by convolution:

[0285]

[0286] S305: Perform baseline removal processing on the denoised pulse wave signal. The baseline removal processing is divided into the following three steps: reference point detection, baseline fitting, and baseline removal operation.

[0287] Peak point detection uses a sliding window to detect the peak of the signal. The width of the sliding window, 2W (i.e., W samples on each side), is set to 0.3s. This window size ensures that most signal peaks are detected without over-smoothing the signal and causing loss of details.

[0288] Starting point detection uses the peak detection results to further determine the reference point by calculating the signal's difference. First, the signal's difference is processed to detect the direction of change and identify the starting point of the upward trend. Similarly, a sliding window size of 0.3s is selected to smooth the detection.

[0289] Cubic spline interpolation is the core step of baseline removal. It fits the baseline of the signal through cubic spline interpolation and removes the baseline from the original signal. The interpolation form used in this method is:

[0290] S(x)=a+bx+cx 2 +dx 3

[0291] The parameters a, b, c, and d are determined by given data points and boundary conditions. The fitted baseline f(t) can be obtained by cubic spline interpolation.

[0292] Finally, the signal x after FIR filtering is FIR Subtract the fitted baseline f(t) from [n] to obtain the signal x after removing the baseline drift 去基线(t). At the same time, in order to maintain the amplitude of the signal, the signal after removing the baseline is added to the mean of the original signal.

[0293]

[0294] S4. Pre-process the rhythm audio used in S2 by first resampling it. Since the degree of synchronization between the rhythm audio and the pulse wave signal needs to be calculated later, the rhythm audio is resampled to the sampling frequency of the pulse wave signal, 100 Hz.

[0295] The resampled rhythm audio is then bandpass filtered to retain a frequency range close to the preprocessed pulse wave signal. This method uses a Butterworth filter with a passband range of 0.1Hz-5Hz and a 4th order.

[0296] The bandpass filtered rhythmic audio is subjected to DC offset removal and normalization to prevent the DC component from affecting the signal's periodicity analysis. Normalization ensures that the amplitudes of different signals are within the same range, facilitating comparison and analysis. Figure 7 The pulse waves obtained from each step of the preprocessing process are shown.

[0297] S5. Extract pulse rate variability (PRV) features from the pre-processed pulse wave signal. First, use the method of S305 to detect the starting point and peak point.

[0298] Then, according to the calculation formula of each PRV feature, features related to cardiovascular health such as PNN50, RMSSD, SDNN, LF, HF, LF / HF, SD1, and SampleEn were extracted respectively.

[0299] S6. Calculate four synchronization features that measure the degree of synchronization between the pulse wave signal and the rhythmic audio signal using the formulas provided in the manual: Pearson correlation coefficients (PCCs), mutual information (MI), cross spectral density (CSD), and phase lock value (PLV). This method uses a sliding window approach to calculate the synchronization features for each pulse wave signal window and each rhythmic audio signal window. To obtain comprehensive analysis results, the sliding window size is set to 0.1s, with a step size of 0.02s. Figure 8 Four synchronization features are shown between a 4 / 4 beat audio with a BPM of 100 and a pulse wave signal measured after a subject played the audio.

[0300] S7. Using the synchronization feature as the independent variable and the PRV feature as the dependent variable, a ridge regression algorithm was used to derive a quantitative evaluation formula for the degree of synchronization between rhythmic audio and pulse wave and cardiovascular health. This formula was then used to provide the experimenter with personalized recommendations for rhythmic audio intervention methods. The following process was used:

[0301] S701. Use multi-factor repeated measures ANOVA to test whether each PRV feature changes significantly under different experimental conditions. The different experimental conditions refer to different groups or different rhythm audio used in experiments of the same group. Screen out PRV features that change significantly under different experimental conditions, such as Figure 9 The F value histogram of each PRV feature is shown.

[0302] S702: A ridge regression model was used to model each PRV feature, yielding a regression formula with the PRV feature as the dependent variable and the synchronization feature as the independent variable. The effect size of the PRV feature and indicators such as the explanatory power and mean square error of the regression model were calculated, and the indicators were standardized. The results are shown in Table 1.

[0303] Table 1 Parameters of each PRV feature regression model and corresponding normalized indicators

[0304]

[0305] S703. Since there is no relevant research to objectively quantify the correlation between PRV characteristics and cardiovascular health, this method uses equal weight, standardized effect size Cohen's d value, standardized model explanatory power R 2 , the inverse of the standardized mean square error 1 / MSE and the comprehensive consideration of Cohen's d, R 2 The regression formula obtained by each PRV feature is weighted by five different weighting methods, namely, 1 / MSE, to obtain a regression formula with cardiovascular health as the dependent variable and the four synchronized features as independent variables, which is called the multi-feature synchronized analysis formula.

[0306] The performance of the five weighting methods was measured by calculating the prediction error of the integrated model under the five different weighting methods. This method used cross-validation to calculate the mean square error (CV-MSE) between the predicted value and the actual value, and used this as the criterion to select the optimal weighting method. The calculation results are shown in Table 2.

[0307] Table 2 Model parameters and CV-MSE values ​​corresponding to different weighting methods

[0308]

[0309] For the fifth weighting method - comprehensive consideration of Cohen's d, R 2 and 1 / MSE, this method traverses Cohen's d, R 2All possible weight combinations of the predicted value and the true value were calculated by grid search and cross validation, and the comprehensive standardized Cohen's d, R 2 The regression formula is weighted with 1 / MSE (a=0.1, b=0.6, c=0.3).

[0310] The final regression formula with cardiovascular health as the dependent variable and synchronization characteristics as the independent variable is as follows:

[0311] Y=0.1566·PCCs+0.3924·MI+0.7063·PLV+0.5747·CSD+0.1079

[0312] PCCs, MI, PLV, and CSD are the four synchronized eigenvalues ​​calculated for each subject after the experiment, after individual differences were removed and normalized using Max-Min normalization. A larger Y indicates better cardiovascular health, so the rhythmic audio and intervention method corresponding to the maximum Y were recommended to the subjects as personalized music intervention methods.

[0313] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be considered as equivalent replacement methods and are included in the scope of protection of the present invention.

Claims

1. A quantitative evaluation method for the cardiovascular effects of rhythmic audio pulse wave synchronization, characterized in that: The following steps are involved: S1. Perform endpoint detection on the monophonic audio of Chinese drums and then generate multiple sets of rhythmic audio based on music theory. S2. Divide the subjects into three groups: a listening intervention group, a playing intervention group, and a sitting meditation control group, and then conduct an experiment; wherein the subjects in the listening intervention group are required to listen to the rhythm audio, and the pulse wave signals of the subjects are collected before and after listening; the subjects in the playing intervention group are required to play the rhythm audio, and the pulse wave signals of the subjects are collected before and after playing; and the pulse wave signals of the subjects in the sitting meditation control group are collected while sitting meditation; S3, preprocessing each collected pulse wave signal; S4, preprocessing the rhythm audio used in the experiment; S5. extracting pulse rate variability features from the pre-processed pulse wave signal; S6. Calculating synchronization characteristics between the pre-processed pulse wave signal and the pre-processed rhythm audio; S7. Using the synchronization feature as the independent variable and the pulse rate variability feature as the dependent variable, a machine learning method is used to obtain a quantitative evaluation formula for cardiovascular health based on the degree of synchronization between the rhythm audio and the pulse wave signal, and a cardiovascular health prediction value is calculated based on the quantitative evaluation formula.

2. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The step S2 comprises the following steps: S201. The subjects were divided into three groups: a listening intervention group, a performance intervention group, and a meditation control group; S202, using a dual-wavelength photoelectric capacitance sensor module to collect pulse wave signals, including: The participants in the listening intervention group were required to collect pulse wave signals for 2 minutes, and then listen to different rhythmic audio clips and collect pulse wave signals again for 2 minutes. The performance intervention group also needed to collect pulse wave signals for 2 minutes, and then play different rhythm audio clips and collect pulse wave signals again for 2 minutes; The subjects in the meditation control group only needed to sit quietly without any music intervention, and their pulse wave signals were directly collected for 2 minutes.

3. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The process of pre-processing the pulse wave signal in step S3 is as follows: S301, performing VMD processing on the pulse wave signal; S302, performing FFT processing on the pulse wave signal after VMD processing to obtain a spectrum diagram of each IMF component; S303, performing wavelet threshold denoising on the high-frequency IMF component, and then performing signal reconstruction; S304, after signal reconstruction, use FIR filter to further remove residual noise; S305 , performing baseline removal processing on the pulse wave signal processed by the FIR filter.

4. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The process of pre-processing the rhythm audio used in the experiment in step S4 is as follows: S401, resampling the rhythm audio used in the experiment in step S2; S402, performing bandpass filtering on the resampled rhythm audio; S403 , performing DC offset removal and normalization processing on the rhythm audio after bandpass filtering.

5. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The step S5 comprises the following steps: S501, performing peak detection on the pre-processed pulse wave signal using a method based on global peak detection and local peak detection; S502: Calculate multiple pulse rate variability features related to cardiovascular health in the pulse wave signal based on the detected peak value.

6. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The step S6 comprises the following steps: S601, calculating the Pearson correlation coefficient between the pulse wave signal and the rhythm audio, which is used to measure the strength and direction of the linear relationship between the pulse wave signal and the rhythm audio; S602, calculating the mutual information between the pulse wave signal and the rhythm audio, which is used to measure the information sharing amount between the pulse wave signal and the rhythm audio; S603, calculating the cross-spectral density between the pulse wave signal and the rhythm audio, which is used to measure the similarity between the pulse wave signal and the rhythm audio in the frequency domain; S604: Calculate the phase lock value between the pulse wave signal and the rhythm audio signal to measure the phase synchronization between the two.

7. The quantitative evaluation method for the cardiovascular effect of rhythmic audio pulse wave synchronization according to claim 1, characterized in that: The step S7 comprises the following steps: S701. Using a multi-factor repeated measures ANOVA method to examine the significance of changes in each pulse rate variability feature under different experimental conditions; wherein the different experimental conditions refer to: different groups or different rhythm audio used in experiments of the same group; S702, using the synchronization feature as an independent variable and the pulse rate variability feature as a dependent variable, and using the ridge regression model to train a regression model for each pulse rate variability feature; S703. Using different weighting methods, weight the regression models obtained from each pulse rate variability feature to obtain a quantitative evaluation formula with the cardiovascular health prediction value as the dependent variable and the synchronization feature as the independent variable, which is called the multi-feature synchronization analysis formula; S704: Calculate the subject's cardiovascular health prediction value based on the multi-feature synchronization analysis formula.

Images