A Respiratory Rate Detection Method Based on Time-Frequency Analysis of Hilbert-Huang Transform

Abstract

The present invention proposes a breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform. First, a WiFi link with strong environmental perception sensitivity is selected to construct a channel state information (CSI) ratio model. Secondly, the filtered CSI ratio time series is projected, and candidate sets of different breathing mode signals are generated by combining amplitude and phase information. Thirdly, the selected breathing modes are subjected to signal variational mode decomposition (VMD) and time-frequency analysis of Hilbert-Huang transform, so as to remove non-human breathing frequency components. On this basis, reconstruction is carried out, and the reconstructed signals are fused using principal component analysis (PCA). Finally, the breathing frequency is calculated through a false peak detection algorithm. The breathing frequency detection algorithm designed by the present invention not only expands the sensing range and eliminates "blind spots" while removing time-varying phase offsets and high-frequency noise and selecting the optimal subcarriers, but also improves the breathing frequency detection accuracy, providing a more accurate detection method for non-contact breathing frequency detection.

Description

Technical Field

[0001] The present invention belongs to the technical field of realizing perception by using commercial WiFi devices, and specifically relates to a breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform, which is used to realize non-contact breathing frequency detection by processing Channel State Information (CSI) signals. Background Art

[0002] With the improvement of living standards, the monitoring of daily breathing frequency has attracted people's attention. In the medical field, breathing frequency is an important physiological index to measure a person's health status. Breathing-related diseases include chronic obstructive pulmonary disease, sleep apnea syndrome, nocturnal hypoventilation syndrome, etc. If the abnormal changes in the breathing frequency of the human body can be detected earlier, early medical advice can be provided for people, thereby improving the health level. The existing research methods are divided into two categories: wearable detection systems and non-contact detection systems. However, commonly used wearable devices, such as polysomnographs, breathing belts, etc., are inconvenient to carry and not suitable for daily use. Therefore, the non-contact-based solution is a feasible option. Non-contact solutions include video-based and wireless technology-based. Video-based solutions are greatly limited in the usage environment due to potential privacy issues and susceptibility to light. While wireless technology can well protect privacy and is insensitive to light. Therefore, wireless technology-based breathing frequency detection has been favored by researchers.

[0003] In the field of wireless technology, radar-based solutions typically use Doppler radar to measure periodic movements generated by breathing, heartbeat, etc. However, the drawback is that dedicated hardware needs to be deployed. Radio Frequency Identification (RFID)-based solutions usually detect breathing by utilizing the radio wave information reflected by tags. However, the tag information of RFID is easily illegally read and even tampered with without permission. Since commercial WiFi devices are inexpensive and currently very popular, detecting the breathing rate based on WiFi signals has become a research hotspot. Moreover, WiFi-based solutions do not require the deployment of additional hardware devices as they are built on existing indoor WiFi devices. Initially, researchers used the Received Signal Strength Indicator (RSSI) information of WiFi to detect the breathing rate. Limited by the coarse-grained nature of the received signal strength index itself and the uncertainty of signal fluctuations, the sensing accuracy and application scope of early WiFi RSSI-based methods were greatly restricted. In addition, there is another available information on WiFi devices, which is CSI. Different from RSSI, CSI can provide more fine-grained wireless channel information at the physical layer and is thus considered an alternative solution for precise sensing. CSI contains channel amplitude and phase information on different subcarriers and can distinguish multipath characteristics. This information reveals the signal scattering, reflection, and power attenuation phenomena of the carrier as the transmission distance changes. By analyzing and studying the changes in CSI, the changes in the physical environment that cause the channel state changes are inferred, thereby realizing non-contact breathing perception.

[0004] However, existing methods require close-range detection. The main reason is that CSI sensing relies on weak reflected signals, and the subtle signal changes caused by breathing are easily masked by noise. In addition, there is no strict clock synchronization for WiFi transceivers, which can lead to random phase offsets. CSI data is affected by the transmission rate, transmission power, and fluctuations in the internal CSI reference level of the Wi-Fi network card, introducing outliers and high-frequency noise. Moreover, existing work also has the problem of blind spots, that is, when a person approaches the sensing device, breathing cannot be effectively detected at certain positions. The short sensing distance and the "blind spot" constraint greatly limit the practical application of existing methods.

[0005] To solve the above problems, the present invention provides a breathing frequency detection method based on Hilbert-Huang transform time-frequency analysis. Aiming at the problems of low accuracy, time-varying phase shift, and blind spots, a WiFi link with strong environmental perception sensitivity is selected. While eliminating the time-varying phase shift using the CSI ratio, the CSI ratio is projected, making full use of amplitude and phase information, and selecting "good" subcarriers. For the problem of high-frequency noise introduced by devices and network cards, special signal decomposition and Hilbert-Huang transform time-frequency analysis are used to remove non-human breathing frequency components and then reconstruct, and principal component analysis is used for fusion, so as to achieve the purpose of expanding the perception range and improving the accuracy of breathing frequency detection, and having strong robustness. Summary of the Invention

[0006] The object of the present invention is to provide a breathing frequency detection method based on Hilbert-Huang transform time-frequency analysis to improve the accuracy of breathing frequency detection in the case of not wearing any devices or cameras and the presence of relative interference of human targets.

[0007] A breathing frequency detection method based on Hilbert-Huang transform time-frequency analysis according to the present invention specifically includes the following steps:

[0008] Step 1: The transmitting device sends a radio frequency signal to the receiving device, and the data is transmitted through orthogonal frequency-division multiplexing (OFDM) so that the signal can be modulated to achieve parallel transmission of multiple subcarriers. The CSI information of the WiFi signal is obtained using an Intel 5300 network card. If the number of WiFi device antennas is I, the number of subcarriers is K, and the number of data packets is N. Then the CSI data matrix can be written as:

[0009]

[0010] Among them, H i represents the CSI data of the i-th (1≤i≤I) WiFi link, represents the CSI data of the k-th (1≤k≤K) subcarrier of the i-th (1≤i≤I) WiFi link, represents the CSI data of the k-th (1≤k≤K) subcarrier of the i-th (1≤i≤I) antenna in the n-th (1≤n≤N) data packet.

[0011]

[0012] Among them, and are respectively the real part and the imaginary part of the i-th antenna and the k-th subcarrier in the n-th data packet, and represents the amplitude and phase of the i-th antenna and the k-th subcarrier in the n-th data packet.

[0013] Step 2: Use the two-way variance to set a threshold to select "sensitive" WiFi links. The specific algorithm process is as follows:

[0014] First, calculate the CSI amplitude variance of each subcarrier on different links. The amplitude variance of I antennas and K subcarriers can be expressed as:

[0015]

[0016] where V i is the vector of the amplitude variances of the K subcarriers of the i-th (1 ≤ i ≤ I) antenna, is the amplitude variance of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as:

[0017]

[0018] where μ i,k is the mean of the amplitudes of the N data packets of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as:

[0019]

[0020] Secondly, set the threshold to ε i = 0.7·max{V i}(1 ≤ i ≤ I), where max{·} represents taking the maximum value of {·}, select the subcarriers on each WiFi link whose amplitude variance is greater than the threshold, and the number of subcarriers is K i , K i = card{k|{V i}> ε i}, where card{·} represents calculating the number of elements of {·}. And take the average value of the amplitude variances of the subcarriers in the set. Then the average amplitude variance of the subcarriers selected by the I antennas is denoted as: v = [v1,…,v i ,…,v I . Among them, v i represents the average amplitude variance of the subcarriers selected by the i-th (1 ≤ i ≤ I) WiFi link, which is a numerical value. Expressed as:

[0021]

[0022] Finally, select the WiFi links with the largest and the second largest average amplitude variances, and rename them as H′1 and H′2.

[0023] Step 2: For the two "sensitive" WiFi links obtained in Step 1, first, the CSI ratio method is used to eliminate the time-varying phase offset. Divide the CSI readings of the two receiving antenna subcarriers. The CSI ratio calculation formula is as follows:

[0024]

[0025] where x k represents the CSI ratio data of the k-th subcarrier, and are the CSI data of the k-th subcarrier of the two selected WiFi links H′1 and H′2 respectively.

[0026] Secondly, the CSI ratio is still a complex number. The CSI data collected through the Monitor mode may have packet loss, and the packet loss rate is about 0.1%-0.05%, which will affect the accuracy of the original data. Further, in order to make the collected data more accurate. First, perform linear interpolation on the original data. Secondly, each subcarrier is processed by Hampel to remove outliers, and the data after removing outliers is smoothed again by the Savitzky-Golay filter. The present invention reconstructs the processed CSI ratio data as:

[0027] X′ = [x′1,…,x′ k ,…,x′ K

[0028] where x k ′ represents the CSI ratio data of the k-th subcarrier after filtering x k .

[0029]

[0030] where, and respectively represent the real part and the imaginary part of the k-th subcarrier x′ k of the CSI ratio data, and respectively represent the real part and the imaginary part of the n-th data packet of the k-th subcarrier of the CSI ratio data.

[0031] Step 3: Project the CSI ratio data, extract the short-term breathing feature Breathing-to-Noise Ratio (BNR), and select the subcarriers that meet the BNR requirements. The specific algorithm process is as follows:

[0032] Project the filtered CSI ratio signal obtained in Step 2, and combine the amplitude and phase information to generate a candidate set representing different breathing mode signals.

[0033] ​First, project the time series x′ of the filtered CSI ratio k (1 ≤ k ≤ K) onto the rotation axis [cosθ sinθ], gradually increase the parameter θ from 0, and generate different candidate sequences:

[0034]

[0035] where y k represents the candidate sequence generated when the rotation angle of the projection coordinate axis of the subcarrier x′ k is θ. Set the range of the projection angle θ to [0, π], and the step size to π / P. A total of P candidate time series are generated for each subcarrier. The different candidate sequences generated after the projection of all subcarriers are expressed as follows:

[0036]

[0037] where represents the p-th (0 ≤ p ≤ P - 1) candidate sequence of the k-th subcarrier, that is: the candidate time series obtained when the projection angle of the k-th subcarrier is θ = (p - 1)·π / P.

[0038] Secondly, calculate the short-term respiratory noise ratio of the signal to extract respiratory features to measure the performance of the combined candidates. The periodicity of the signal pattern within a period of time represents its ability to sense respiration. The specific calculation steps are as follows:

[0039] First, if the sampling rate f s = 100 Hz, use a window length of 12, corresponding to 1200 samples, pad with 6992 zero-value samples to get 8192 samples, and perform Fourier transform on all samples. Secondly, screen out the maximum energy within the range of human respiration (10 - 37 bpm), which is the energy of the respiratory signal. Take the ratio of this energy to the total energy in the frequency domain to obtain the short-term respiratory noise ratio.

[0040] Find the respiratory feature BNR of all candidate time series of K subcarriers, which is expressed as follows:

[0041]

[0042] where represents the BNR value of the candidate time series obtained when the projection angle of the k-th subcarrier is θ = (p - 1)·π / P, which is a numerical value.

[0043] Then, find the maximum BNR value of the k-th subcarrier under different projection axes, that is to form a set {b k}. Take the maximum value in the set {b k} as b, and set the threshold to δ = 0.7b, and select the set {bk b greater than the threshold value in {...} k The corresponding sub - carrier candidate sequence. The number of selected sub - carriers is: S = card{k|{b k}>δ}, and find the candidate time sequence of the k - th sub - carrier in the set {k|{b k}>δ}. The selection result is as follows:

[0044] Y′ = [y′1,…,y′ s ,…,y′ S

[0045] where y′ s represents the candidate time sequence corresponding to the s - th (1≤s≤S) selected sub - carrier.

[0046] Step 4: Signal Variational Mode Decomposition (VMD) and Hilbert - Huang transform time - frequency analysis are performed to remove non - human breathing frequency components for reconstruction. The specific algorithm flow is as follows:

[0047] First, perform variational mode decomposition on the S sub - carriers selected in Step 3 respectively:

[0048]

[0049] where y′ s (t) is the s - th sub - carrier in Y′, M is the number of modal components, and u m (t) is the m - th (1≤m≤M) modal component signal.

[0050] Secondly, construct the analytic signal for each modal component. The analytic signal is based on the combination of the original signal and its Hilbert transform, and the instantaneous frequency and instantaneous energy are obtained according to the analytic signal:

[0051]

[0052] where z m represents the analytic signal, represents the Hilbert transform, a m and θ m respectively represent the instantaneous amplitude and instantaneous phase of the m - th modal component, |a m (t)| 2 and w m respectively represent the instantaneous energy and instantaneous frequency of the m - th modal component.

[0053] ​Then, perform Hilbert-Huang transform time-frequency analysis on each modal component using the instantaneous frequency and instantaneous energy obtained from the above formula. Remove the first two high-frequency components that are not related to the breathing frequency, and reconstruct the modal components within the remaining breathing range. The reconstructed time series of the sth (1 ≤ s ≤ S) subcarrier is expressed as:

[0054]

[0055] All the reconstructed subcarrier combinations are denoted as Y″:

[0056] Y″ = [y″1, …, y″ s , …, y″ S

[0057] Step Five: Signal fusion. The specific algorithm is as follows:

[0058] Fuse the signals reconstructed in Step Four using Principal Component Analysis (PCA). PCA transforms the original feature vectors into a set of linearly independent principal components. This algorithm first centrally calculates the covariance matrix of all samples, and then solves for the eigenvalues and eigenvectors through singular value decomposition. Finally, the main breathing signal is obtained using the largest eigenvalue and eigenvector as the principal component analysis signal:

[0059]

[0060] where, a n represents the nth (1 ≤ n ≤ N) data packet of the fused principal component analysis signal .

[0061] Step Six: Perform peak detection on the principal component breathing signal obtained in Step Five and extract the breathing rate. The specific algorithm process is as follows:

[0062] First, obtain the local peak set for the principal component breathing signal : Maxset = {τ j , 1 ≤ j ≤ J}. The verification window is set to: M;

[0063] Secondly, for each peak point τ j find the location and peak locs := location(τ j ); amp := amplitude(τ j );

[0064] Thirdly, for each point in this range, first determine whether n is within the range of 1 < n < N. If it is, by comparing ... ​For each point in this range, determine whether there is a point larger than the current peak point. If: amp < x(n), delete τ from MarSet k ; and update MarSet;

[0065] Finally, from the updated Maxset = {τ l , 1 ≤ l ≤ L}, obtain the total time interval between adjacent true peaks of true peaks, denoted as sum; from the formula: obtain the breathing rate, f s is the sampling rate.

[0066] Advantageous Effects

[0067] This method first selects two WiFi links with larger fluctuation degrees from the obtained three WiFi links through double variance, and eliminates the random phase offset of commercial WiFi over time by taking their ratio. Secondly, project the CSI ratio, extract the breathing feature BNR, and select subcarriers that meet the threshold for VMD signal decomposition and Hilbert-Huang transform for time-frequency analysis, thereby removing non-human breathing frequency components. On this basis, reconstruction is performed, and the reconstructed signals are fused using principal component analysis. Finally, peak detection is used to extract the breathing rate.

[0068] The present invention provides a breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform. While removing time-varying phase offset and a large amount of high-frequency noise and selecting the optimal subcarriers, it not only expands the sensing range, but also eliminates "blind spots" and improves the detection accuracy of breathing frequency, providing a more accurate detection method for non-contact breathing frequency detection. Description of the Drawings

[0069] Figure 1 is a block diagram of the breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform of the present invention Specific Embodiment

[0070] The object of the present invention is to provide a breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform without the need to wear any devices or cameras and in the presence of relative interference of the human target, so as to improve the detection accuracy of breathing frequency.

[0071] A breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform according to the present invention specifically includes the following steps:

[0072] Step 1: The transmitting device sends a radio frequency signal to the receiving device, and the data is transmitted through Orthogonal Frequency-Division Multiplexing (OFDM) so that the signal can be modulated to achieve parallel transmission of multiple subcarriers. The CSI information of the WiFi signal is obtained using an Intel 5300 network card. If the number of antennas of the WiFi device is I, the number of subcarriers is K, and the number of data packets is N. Then the CSI data matrix can be written as:

[0073]

[0074] Among them, H i represents the CSI data of the i-th (1 ≤ i ≤ I) WiFi link, represents the CSI data of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) WiFi link, represents the CSI data of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna in the n-th (1 ≤ n ≤ N) data packet.

[0075]

[0076] Among them, and are respectively the real part and the imaginary part of the i-th antenna and the k-th subcarrier in the n-th data packet, and represent the amplitude and phase of the i-th antenna and the k-th subcarrier in the n-th data packet.

[0077] Step 2: Use the two-way difference to set a threshold to select "sensitive" WiFi links. The specific algorithm process is as follows:

[0078] First, calculate the variance of the CSI amplitude of each subcarrier on different links. The variance of the amplitude of I antennas and K subcarriers can be expressed as:

[0079]

[0080] Among them, V i is the vector of the variance of the amplitude of K subcarriers of the i-th (1 ≤ i ≤ I) antenna, is the variance of the amplitude of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as:

[0081]

[0082] Among them, μ i,k is the mean value of the amplitudes of N data packets of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as:

[0083]

[0084] Secondly, set the threshold to ε i = 0.7·max{V i (1 ≤ i ≤ I), where max{·} represents taking the maximum value of {·}. Select the subcarriers whose amplitude variances of each WiFi link are greater than the threshold, and the number of subcarriers is K i , K i = card{k|{V i} > ε i}, where card{·} represents calculating the number of elements of {·}. And take the average value of the amplitude variances of the subcarriers in the set. Then, the average amplitude variances of the subcarriers selected by I antennas are denoted as: v = [v1,…,v i ,…,v I . Among them, v i represents the average amplitude variance of the subcarriers selected by the i-th (1 ≤ i ≤ I) WiFi link, which is a numerical value. It is expressed as:

[0085]

[0086] Finally, select the WiFi links with the largest and the second largest average amplitude variances, and rename them as H′1 and H′2.

[0087] Step 2: For the two "sensitive" WiFi links obtained in Step 1, first, use the CSI ratio method to eliminate the time-varying phase offset. Divide the CSI readings of the subcarriers of the two receiving antennas. The CSI ratio calculation formula is as follows:

[0088]

[0089] where x k represents the CSI ratio data of the k-th subcarrier, and are the CSI data of the k-th subcarrier of the two selected WiFi links H′1 and H′2 respectively.

[0090] Secondly, the CSI ratio is still a complex number. The CSI data collected through the Monitor mode may have packet loss, and the packet loss rate is about 0.1% - 0.05%, which will affect the accuracy of the original data. Further, in order to make the collected data more accurate. First, perform linear interpolation on the original data. Secondly, pass each subcarrier through the Hampel to remove outliers, and then pass the data after removing outliers through the Savitzky-Golay filter for smoothing. The present invention reconstructs the processed CSI ratio data as:

[0091] X′ = [x′1,…,x′k , …, x′ K

[0092] Among them, x k ′ represents the CSI ratio data of the k-th subcarrier after filtering for x k .

[0093]

[0094] Among them, and respectively represent the real part and the imaginary part of the k-th subcarrier x′ of the CSI ratio data k , and respectively represent the real part and the imaginary part of the n-th data packet of the k-th subcarrier of the CSI ratio data.

[0095] Step 3: Project the CSI ratio data, extract the short-term breathing feature breathing-to-noise ratio (BNR), and select the subcarriers that meet the BNR requirements. The specific algorithm process is as follows:

[0096] Project the filtered CSI ratio signal obtained in Step 2, and generate a candidate set representing different breathing mode signals by combining the amplitude and phase information.

[0097] First, project the time series x′ of the filtered CSI ratio k (1 ≤ k ≤ K) onto the rotation axis [cosθ sinθ], gradually increase the parameter θ from 0, and generate different candidate sequences:

[0098]

[0099] Among them, y k represents the candidate sequence generated when the rotation angle of the subcarrier x′ k on the projection coordinate axis is θ. Set the projection angle θ range to [0, π], and the step size is set to π / P. A total of P candidate time series are generated for each subcarrier. The different candidate sequences generated after projecting all subcarriers are represented as follows:

[0100]

[0101] Among them, represents the p-th (0 ≤ p ≤ P - 1) candidate sequence of the k-th subcarrier, that is: the candidate time series obtained when the projection angle of the k-th subcarrier is θ = (p - 1)·π / P.

[0102] ​Secondly, the short-term breathing noise ratio of the signal is calculated to extract the breathing feature to measure the performance of the combination candidate. The periodicity of the signal pattern over a period of time represents its ability to detect breathing. The specific calculation steps are as follows:

[0103] First, if the sampling rate f s =100Hz, using a window length of 12, corresponding to 1200 samples. 6992 zero-valued samples were padded with zeros, resulting in 8192 samples. All samples were Fourier transformed. Next, the maximum energy within the human breathing range (10-37bpm) was filtered out, representing the respiratory signal energy. This energy was then compared to the total energy in the frequency domain to obtain the short-term respiratory noise ratio.

[0104] Find the breathing feature BNR of all candidate time series of K subcarriers, which is expressed as follows:

[0105]

[0106] in, It represents the BNR value of the candidate time series obtained when the projection angle of the k-th subcarrier is θ=(p-1)·π / P, which is a numerical value.

[0107] Then, find the maximum BNR value of the kth subcarrier under different projection axes, that is, Composition set {b k}. Take the set {b k} is b, and the threshold is set to δ=0.7b, select the set {b k} is greater than the threshold b k The corresponding subcarrier candidate sequence. The number of selected subcarriers is: S = card {k | {b k}>δ}, and find the corresponding set {k|{b k The candidate time sequence of the kth subcarrier in}>δ}. The selection results are as follows:

[0108] Y′=[y′1,…,y′ s ,…,y′ S ]

[0109] Among them, y′ s Indicates the candidate time sequence corresponding to the selected s-th (1≤s≤S) subcarrier.

[0110] Step 4: Perform signal variational mode decomposition (VMD) and Hilbert-Huang transform time-frequency analysis to remove non-human respiratory frequency components for reconstruction. The specific algorithm flow is as follows:

[0111] First, perform variational mode decomposition on the S subcarriers selected in Step 3 respectively:

[0112]

[0113] Among them, y′ s (t) is the sth subcarrier in Y′, M is the number of modal components, and u m (t) is the mth (1 ≤ m ≤ M) modal component signal.

[0114] Secondly, construct the analytic signal for each modal component. The analytic signal is based on the combination of the original signal and its Hilbert transform, and the instantaneous frequency and instantaneous energy are obtained according to the analytic signal:

[0115]

[0116] Among them, z m represents the analytic signal, represents the Hilbert transform, a m and θ m respectively represent the instantaneous amplitude and instantaneous phase of the mth modal component, |a m (t)| 2 and w m respectively represent the instantaneous energy and instantaneous frequency of the mth modal component.

[0117] Then, perform Hilbert-Huang transform time-frequency analysis on each modal component using the instantaneous frequency and instantaneous energy obtained above, remove the first two high-frequency components unrelated to the breathing frequency, and reconstruct the modal components within the remaining breathing range. The reconstructed time series of the sth (1 ≤ s ≤ S) subcarrier is expressed as:

[0118]

[0119] The combination of all the reconstructed subcarriers is denoted as Y″:

[0120] Y″ = [y″1,…,y″ s ,…,y″ S

[0121] Step 5: Signal fusion. The specific algorithm is as follows:

[0122] Fuse the signals after reconstruction in Step 4 using principal component analysis (PCA). PCA transforms the original feature vectors into a set of linearly independent principal components. This algorithm first centrally calculates the covariance matrix of all samples, and then solves the eigenvalues and eigenvectors through singular value decomposition. Finally, the main breathing signal is obtained using the largest eigenvalue and eigenvector as the principal component analysis signal:​

[0123]

[0124] where a n represents the nth (1 ≤ n ≤ N) data packet of the fused principal component analysis signal of

[0125] Step Six: Perform peak detection on the principal component respiration signal obtained in Step Five to extract the respiration rate. The specific algorithm flow is as follows:

[0126] First, for the principal component respiration signal obtain the local peak set: Maxset = {τ j , 1 ≤ j ≤ J}. The verification window is set to: M;

[0127] Second, for each peak point τ j find the position and peak locs: = location(τ j ); amp: = amplitude(τ j );

[0128] Third, for each point in this range first determine whether n is within the range of 1 < n < N. If so, by comparing each point in this range to determine whether there is a point larger than the current peak point. If: amp < x(n), delete τ from MarSet k ; and update MarSet;

[0129] Finally, from the updated Maxset = {τ l , 1 ≤ l ≤ L}, obtain the total time interval between adjacent true peaks of the true peaks, denoted as sum; from the formula: obtain the respiration rate, f s is the sampling rate.

Claims

1. A breathing frequency detection method based on time-frequency analysis of Hilbert-Huang transform, comprising the following steps: Step 1: The transmitting device sends a radio frequency signal to the receiving device, and the data is transmitted through orthogonal frequency-division multiplexing (OFDM) so that the signal can be modulated to achieve parallel transmission of multiple subcarriers; the WiFi signal CSI information is obtained by using an Intel 5300 network card. If the number of antennas of the WiFi device is I, the number of subcarriers is K, and the number of data packets is N; then the CSI data matrix is written as: Among them, H i represents the CSI data of the i-th (1 ≤ i ≤ I) WiFi link, represents the CSI data of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) WiFi link, represents the CSI data of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna in the n-th (1 ≤ n ≤ N) data packet; Among them, and are the real part and the imaginary part of the \(i\)-th antenna, \(k\)-th subcarrier in the \(n\)-th data packet respectively, and represent the amplitude and phase of the \(i\)-th antenna, \(k\)-th subcarrier in the \(n\)-th data packet; Step 2: Use the two-sided variance to set a threshold to select "sensitive" WiFi links; the specific algorithm flow is as follows: First, calculate the CSI amplitude variance of each subcarrier on different links; the amplitude variances of I antennas and K subcarriers are expressed as: Among them, V i is the vector of the amplitude variances of K subcarriers of the i-th (1 ≤ i ≤ I) antenna, is the amplitude variance of the k-th (1 ≤ k ≤ K) subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as: where, μ i,k is the mean of the amplitudes of N data packets of the k-th subcarrier of the i-th (1 ≤ i ≤ I) antenna, expressed as: Secondly, set the threshold to ε i = 0.7·max{V i}(1 ≤ i ≤ I), where max{·} represents taking the maximum value of {·}, select the subcarriers whose amplitude variance of each WiFi link is greater than the threshold, and the number of subcarriers is K i , K i = card{k|{V i}} > ε i}}, where card{·} represents calculating the number of elements of {·}; and take the average value of the amplitude variances of the subcarriers in the set, then the average amplitude variance of the subcarriers selected by I antennas is denoted as: v = [v1, …, v i , …, v I ; where v i represents the average amplitude variance of the subcarriers selected by the i-th (1 ≤ i ≤ I) WiFi link, which is a numerical value; it is expressed as: Finally, select the WiFi links with the largest and the second largest average amplitude variances, and rename them as H1′ and H2′; Step 3: For the two "sensitive" WiFi links obtained in Step 2; first, use the CSI ratio method to eliminate the time-varying phase offset, divide the CSI readings of the subcarriers of two receiving antennas, and the CSI ratio calculation formula is as follows: where x k represents the CSI ratio data of the k-th subcarrier, and are the CSI data of the k-th subcarrier of the two selected WiFi links H1′ and H2′ respectively; Secondly, the CSI ratio is still a complex number. In order to make the collected data more accurate; first, perform linear interpolation on the original data; secondly, remove outliers of each subcarrier through Hampel, and the data after removing outliers is smoothed again through a Savitzky-Golay filter; the processed CSI ratio data is reconstructed as: X′ = [x′1, …, x′ k , …, x′ K ​ where x k ′ is represented as x k the CSI ratio data of the k-th subcarrier after filtering; Among them, and respectively represent the real part and the imaginary part of the k-th subcarrier x′ k of the CSI ratio data, and respectively represent the real part and the imaginary part of the n-th data packet of the k-th subcarrier of the CSI ratio data; Step 4: Project the CSI ratio data, extract the breathing feature short-term breathing-to-noise ratio (BNR), and select subcarriers that meet the BNR requirements; the specific algorithm flow is as follows: Project the filtered CSI ratio signal obtained in Step 3, and combine the amplitude and phase information to generate a candidate set representing different breathing mode signals; First, project the time series \(x'\) of the filtered CSI ratio k (1 ≤ k ≤ K) onto the rotation axis \([\cos\theta\ \sin\theta]\), and gradually increase the parameter \(\theta\) from 0 to generate different candidate sequences: where y k represents the candidate sequence generated when the rotation angle of the projection coordinate axis is θ for subcarrier x′ k ; setting the range of the projection angle θ to be [0, π], with a step size of π / P, a total of P candidate time sequences are generated for each subcarrier; the different candidate sequences generated after projection of all subcarriers are represented as follows: Among them, represents the p-th (0 ≤ p ≤ P - 1) candidate sequence of the k-th subcarrier, that is, the candidate time sequence obtained when the projection angle of the k-th subcarrier is θ = (p - 1)·π / P; Secondly, calculate the signal short-term breathing-to-noise ratio to extract the breathing feature to measure the performance of the combined candidates; the periodicity of the signal mode within a period of time represents its ability to sense breathing, and the specific calculation steps are as follows: First, if the sampling rate f s = 100 Hz, a window length of 12 is used, corresponding to 1200 samples. Zero-padding is performed with 6992 zero-valued samples to obtain 8192 samples, and Fourier transform is performed on all samples. Secondly, the maximum energy is selected within the human breathing frequency range of 10 - 37 bpm, which is the energy of the breathing signal. The ratio of this energy to the total energy in the frequency domain is calculated to obtain the short-term breathing noise ratio; Obtain the breathing feature BNR of all candidate time series of K subcarriers, which is expressed as: Among them, represents the BNR value of the candidate time series obtained when the k-th subcarrier is at the projection angle θ = (p - 1)·π / P, which is a numerical value; Then, find the maximum BNR value of the k-th subcarrier under different projection axes, that is to form a set {b k}; take the maximum value in the set {b k} as b, and set the threshold as δ = 0.7b. Select the subcarrier candidate sequences corresponding to the b k greater than this threshold in the set {b k}; the number of selected subcarriers is: S = card{k|{b k} > δ}, and find the candidate time sequence of the k-th subcarrier corresponding to the set {k|{b k} > δ}; the selection result is as follows: Y′ = [y1′, …, y s ′, …, y′ S ​ where y s ′ represents the candidate time series corresponding to the selected sth (1 ≤ s ≤ S) subcarrier; Step 5: Signal variational mode decomposition (VMD) and Hilbert-Huang transform time-frequency analysis are performed to remove non-human breathing frequency components for reconstruction; the specific algorithm flow is as follows: First, perform variational mode decomposition on the S subcarriers selected in Step 4 respectively: where y s ′(t) is the s-th subcarrier in Y′, M is the number of modal components, and u m (t) is the m-th (1 ≤ m ≤ M) modal component signal; Secondly, construct the analytic signal for each modal component. The analytic signal is based on the combination of the original signal and its Hilbert transform, and the instantaneous frequency and instantaneous energy are obtained according to the analytic signal: where z m represents the analytic signal, represents the Hilbert transform, a m and θ m respectively represent the instantaneous amplitude and instantaneous phase of the m-th mode component, |a m (t)| 2 and w m respectively represent the instantaneous energy and instantaneous frequency of the m-th mode component; Then, the instantaneous frequency and instantaneous energy obtained by the above formula are used to perform Hilbert-Huang transform time-frequency analysis on each modal component. The first two high-frequency components unrelated to the breathing frequency are removed, and the modal components within the remaining breathing range are reconstructed. The reconstructed time series of the sth (1 ≤ s ≤ S) subcarrier is expressed as: All the reconstructed subcarrier combinations are denoted as Y″: Y″ = [y″1, …, y″ s , …, y″ S ​ Step Six: Signal fusion; the specific algorithm is as follows: Use principal component analysis (PCA) to fuse the signals after reconstruction in Step Five. PCA transforms the original feature vectors into a set of linearly independent principal components. This algorithm first centrally calculates the covariance matrix of all samples, then solves the eigenvalues and eigenvectors through singular value decomposition. Finally, the main breathing signal is obtained using the largest eigenvalue and eigenvector as the principal component analysis signal: Among them, a n represents the nth (1 ≤ n ≤ N) data packet of the fused principal component analysis signal ; Step Seven: Perform peak detection on the principal component respiratory signal obtained in Step Six to extract the respiratory rate. The specific algorithm flow is as follows: First, for the principal component respiration signal Obtain the local peak set: Maxset = {τ j , 1 ≤ j ≤ J}, where J is the number of peak points; the verification window is set to: M; Secondly, for each peak point τ j find the location and peak locs := location(τ j ); amp := amplitude(τ j ); Again, for each point in this range, first determine whether n is within the range of 1 < n < N; if so, by comparing the values of each point in this range to determine whether there is a point larger than the current peak point. If: Delete τ from Maxset j ; and update Maxset; Finally, from the updated Maxset = {τ l , 1 ≤ l ≤ L}, the total time interval between adjacent true peaks of true peaks is obtained, denoted as sum; where L is the number of updated peak points; according to the formula: the respiratory rate f is obtained, s where f is the sampling rate.

Images