The relationship between the frequency of the Chirp signal emitted by the FMCW millimeter-wave radar and time is expressed by Equation 1: f t = f c + B T c t (1) Where f c is the frequency at the initial moment when the chirp signal is transmitted, T c is the duration of the chirp signal frequency sweep, and B is the frequency modulation bandwidth. ∅ 1 is initial phase. Therefore, the FMCW millimeter wave radar transmission signal can be expressed as Equation 2: x T t = A T ⁡ cos 2 π f c t + π B T c t 2 + ∅ 1 (2)

Assuming the distance between the radar and the target person is R , the radar emits a chirp signal that reaches the chest of the person and then reflects back. The round trip time is expressed by Equation 3: t d = 2 R c (3)

Where c is the speed of light, then the expression of the reflected echo signal is expressed by Equation 4: x R t = x T t − t d = A R ⁡ cos 2 π f c t − t d + π B T c t − t d 2 + ∅ 2 (4)

The internal circuit of the radar mixes the transmitted signal and the reflected echo signal to obtain an IF signal in expression (Equation 5): x I F t = x T t · x R t (5)

From the product and difference formula of trigonometric functions, we can get: x I F t = A T · A R 2 cos 4 π B R c T c t + 4 π f c R c + 4 π B R 2 c 2 T c t 2 + ∅ , ∅ = ∅ 1 − ∅ 2 (6)

Since 4 π B R 2 is much smaller than c 2 T c , and the chirp signal period is very short, t 2 is also close to zero. Therefore, the quadratic term 4 π B R 2 c 2 T c t 2 in expression (Equation 6) is almost equal to zero and can be ignored. At the same time, let A I F = A T · A R / 2 , then the simplified IF signal expression is Equations 7, 8: x I F t = A I F ⁡ cos 4 π B R c T c t + 4 π f c R c + ∅ (7) x I F t = A I F ⁡ cos f I F t + φ I F + ∅ (8) f I F = 2 B c · T c · R (9) φ I F = 4 π f c R c = 4 π λ · R (10) Where f I F is the frequency of the IF signal in Equation 9. φ I F is the phase of the IF signal. λ is the wavelength of the electromagnetic wave emitted by the millimeter wave radar. After the phase expression (Equation 9) is differentiated, the phase change of the IF signal with respect (Equation 10) to the distance is expressed as: Δ φ I F = 4 π λ · Δ R (11) Where Δ φ I F is the phase change of the IF signal, and Δ R is the change in chest vibration displacement. From Equation 11, it is evident that when the initial frequency of the millimeter-wave radar is 77 GHz, the wavelength λ ≈ 4 mm, resulting in a phase shift of 0.314 radians for a chest displacement of 0.1 mm induced by the heartbeat. Thus, we can capture the minute chest vibrations caused by heartbeat and respiration through significant phase variations. This principle forms the foundation for using FMCW millimeter-wave radar to measure human vital sign signals.

When the radar transmits chirp signals toward the human body from a certain angle, the reflected waves from the thoracic cavity may scatter in multiple directions, and only a portion of these echoes are received by the radar antennas. Beamforming can effectively enhance the echo signals from the desired direction while suppressing clutter and interference from other directions. Beamforming relies on a multiple-input multiple-output (MIMO) antenna array. The AWR1642 mm-wave radar used in this study can be configured as a 2-transmit, 4-receive (2T4R) antenna array, which is virtually extended to a 1-transmit, 8-receive (1T8R) linear array using TDM. To enhance the echo signal from the thoracic cavity and reduce the output noise and interference power, we adopt the Capon beamforming method (Capon, 1969). Capon beamforming, also known as minimum variance distortionless response (MVDR) beamforming, can effectively suppress clutter signals from non-target directions while enhancing signals from the desired direction. It provides superior angular resolution performance compared to conventional beamforming methods (Marino and Chau, 2005), and although its resolution is lower than that of the multiple signal classification (MUSIC) algorithm (Schmidt, 1986), Capon beamforming offers significantly lower computational complexity. In the 1T8R uniformly spaced linear antenna array, the IF signal received at each antenna element experiences a time delay relative to the adjacent previous antenna. Let the time delay be denoted as τ , d be the inter-element spacing, and θ be the angle of arrival (AoA) of the incoming signal, then τ is expressed by Equation 12: τ = d s i n θ c (12)

The slow-time complex sequence corresponding to the range bin of the human target can be reasonably modeled as a narrowband signal. Due to the different propagation paths between the target and individual receiving antennas, the signal phase exhibits variations across antenna elements within the inter-antenna propagation delay range. Let the IF signal at the target range bin have a frequency of f h , an amplitude of A h , and a phase of φ h . At a given time instant t , assuming the presence of zero-mean Gaussian noise n m t , and denoting the signal amplitude received by the m-th antenna as A h m , the signal received by the m-th antenna can be expressed as Equation 13: x m t = A h m · ⁡ cos 2 π f h t − m τ + φ h + n m t = A h m · ⁡ cos 2 π f h t + φ h · ⁡ cos − 2 π f h · m τ + n m t   m = 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 . (13)

Let s t = A h m · ⁡ cos 2 π f h t + φ h , a m θ = cos − 2 π f h · m d s i n θ c a θ = a 0 θ , a 1 θ , a 2 θ , … … … . . , a 7 θ T , n t = n 0 t , n 1 t , n 2 t , … … … . . , n 7 t T

The target range gate signals (as snapshot signals) from the eight virtual antennas are arranged into a vector form as Equation 14: X i = s i a θ + n i , i = 0 , 1 , 2 , 3 , … . . N (14)

Where N is the number of snapshots. The calculation formula of Capon beamforming weight vector is as Equation 15: w θ = R X X − 1 θ a θ a H θ R X X − 1 a θ (15)

The spatial spectrum function is expressed by Equation 16: P θ = 1 a H θ R X X − 1 a θ (16) w θ is the optimal weight vector, R X X is the snapshot signals covariance matrix, and R X X − 1 represents the inverse matrix of the snapshot signal covariance matrix. There is expressed by Equation 17: R X X = 1 N ∑ i = 1 N   X i X H i (17)

The weighted sum of the target range gate slow-time complex signals of the eight virtual antennas is performed, and the expression of the enhanced slow-time complex signal is obtained as Equation 18: Y = w H X (18)

As shown in Formula 18, the slow-time complex sequence signals of the target range gates from all eight antenna elements can be weighted summed to enhance the signal from the direction of the human body while suppressing clutter signals from other directions. However, the direction corresponding to the peak of the spatial spectrum function in Equation 16 does not necessarily coincide with the direction of the highest SNR of the heartbeat signal. In some cases, the spectral peak may occur in the direction of clutter or other objects. Furthermore, even if the spectral peak appears in the direction of the human body, it does not guarantee that this is the direction with the highest SNR for the heartbeat signal. For most individuals, whose body height range is about 1.55–1.9 m, only the region around the heart yields heartbeat signals with the highest SNR. Additionally, in some individuals who exhibit abdominal breathing, significant motion occurs in the abdominal region. Conventional beamforming methods combined with respiration detection (Wang et al., 2021b) may misidentify the abdominal region as the target direction, which still does not correspond to the direction of the strongest heartbeat signal. To address these limitations, we propose a new azimuth search method to replace the traditional spatial spectrum function approach, aiming to identify the direction with the highest SNR of the SCG signal. At the beginning of the measurement, we collect the first 20 s of data and process them using the procedure outlined in Figure 1 to extract the SCG signal. To assess the quality of SCG signals, we propose a method based on dynamic time warping (DTW) (Sakoe and Chiba, 1978). Specifically, a 20 s segment with high-quality SCG signal, selected from all radar-measured SCG recordings, is designated as the template signal. The DTW distance is then calculated between the SCG signal under evaluation and the template signal. The beamforming angle is selected by minimizing the DTW distance between the extracted SCG signal and the reference template. A smaller DTW distance indicates a higher similarity between the two signals, thereby reflecting better signal quality. Conversely, a larger DTW distance suggests poorer signal quality. Let the selected template signal be denoted as S C G t e m p l a t e , and the SCG signal to be evaluated as S C G . Prior to DTW computation, both S C G and S C G t e m p l a t e are standardized. Here, σ represents calculating the variance of the SCG sequence as shown in (Equations 19, 20): S C G = S C G − S C G ¯ σ S C G (19) S C G t e m p l a t e = S C G t e m p l a t e − S C G t e m p l a t e ¯ σ S C G t e m p l a t e (20) S Q I = 1 D T W S C G , S C G t e m p l a t e (21)

To identify the direction with the highest SNR of the SCG signal, we perform a directional search over the angular range of [-60°,60°] with a step size of 1°. For each angle, a corresponding beamforming weight vector w is computed, and based on this, a SQI value is calculated, as shown in Equation 21. This process yields a total of 121 angle-SQI pairs. Among all scanned directions, the angle θ o p t corresponding to the maximum SQI value is selected, which we consider to be the direction with the highest SCG signal SNR. Upon determining the cardiac orientation, the SCG template is no longer utilized, and the search for the optimal orientation via improved Capon beamforming is terminated. The resulting beamforming azimuth is subsequently fixed for all further processing. After determining θ o p t , the optimal beamforming weight vector w o p t is calculated according to Formula 15, and the slow-time complex signal in the direction of the heart is obtained using Formula 18. This slow time series is then demodulated using a method called MDACM (Xu et al., 2021) to produce a phase series containing respiratory and cardiac motion components, i.e., the original vital sign signal.

The original vital sign signal, primarily contains the heartbeat signal, respiratory signal, and other noise. Due to the greater displacement of the respiratory signal compared to that of the heartbeat signal, the heartbeat signal is often overshadowed by the respiratory signal. Additionally, the third harmonic frequency band of the respiratory signal overlaps with the frequency band of the heartbeat signal, complicating the separation of these two signals. Designing an algorithm to effectively separate the heartbeat signal from the respiratory signal remains a challenge. Commonly used methods, such as the BPF, typically set parameters of 0.1 Hz–0.5 Hz for the respiratory signal and 0.7 Hz–3 Hz for the heartbeat signal. However, these methods yield only moderate results and fail to address the issue of frequency band overlap caused by the respiratory signal’s third harmonic. Standard binary wavelet analysis decomposes signals layer by layer along the low-frequency direction, which limits its ability to finely segment the frequency bands, particularly in the high-frequency region. The heartbeat signal contains both high-frequency and low-frequency components, and binary wavelet analysis is insufficient for capturing the entire frequency spectrum of the heartbeat signal. Wavelet packet analysis offers a more detailed and nuanced approach to signal analysis. Unlike standard binary wavelet analysis, which analyzes low and high-frequency components separately, wavelet packet decomposition simultaneously considers both, providing a more accurate local analysis. This method allows for finer division of the time-frequency plane and offers superior resolution for the high-frequency components of the signal. Moreover, wavelet packet analysis introduces the concept of optimal basis selection, allowing for adaptive selection of the most suitable basis functions that align with the characteristics of the signal, thereby significantly enhancing the effectiveness of the signal analysis. In this study, we employ WPT using the “db45” wavelet to break down the original vital sign signal into a series of finely segmented wavelet packet coefficients within specific frequency ranges. The “db45″wavelet has a 45th-order vanishing moment. A higher vanishing moment means a higher concentration in the frequency domain, which can better capture the high-frequency details of the signal by selecting the wavelet packet coefficients corresponding to the SCG signal frequency range, we reconstruct the SCG signal. The original vital sign signal, sampled at a rate of f s = 200 Hz over a duration of 50 s, undergoes six levels of decomposition using the WPT. The first 16 wavelet coefficients from the sixth level are illustrated in Figure 2. According to the WPT theory, the frequency range of the i-th wavelet coefficient C i of the sixth layer is expressed by Equation 22 C i ∼ i − 1 · f s / 2 2 6 , i · f s / 2 2 6 (22)

Where f s = 200 Hz is the sampling rate of the original vital signs signal. Further analysis of the wavelet packet coefficient energy percentage was conducted, with results illustrated in Figure 3. The energy distribution shows that the first five wavelet packet coefficients account for a significantly higher energy percentage compared to the remaining coefficients. According to Equation 22, the frequency range for the first wavelet packet coefficient spans from 0 Hz to 1.5625 Hz, primarily containing respiratory components and motion artifacts. The second through fourth wavelet packet coefficients mainly consist of harmonic components of respiratory and motion artifacts. From the sixth wavelet packet coefficient to approximately the 12th, a peak in energy appears, indicating that these six coefficients primarily capture the heartbeat signal components. Although the first five wavelet packet coefficients may contain some heartbeat signal components, the heartbeat signal content in these coefficients is significantly lower than that of respiratory components, motion artifacts, and their harmonics. Literature (Taebi et al., 2019) reports that the primary frequency range of the SCG signal is below 25 Hz. Additionally, literature Equation 22, the frequency range for these coefficients is 7.8125 Hz–18.750 Hz, which closely aligns with the SCG frequency range described in (Zhang et al., 2023). Therefore, the selected wavelet packet coefficients effectively reconstruct the SCG signal. The time-domain waveform and the synchronized ECG are shown in Figure 4. Figure 5 shows the SCG signal extracted from the radar signal and the SCG signal measured from the synchronously acquired cardiac acceleration signal. As shown in Figure 5, the SCG signal measured by radar exhibits a high degree of similarity to the SCG signal measured by an accelerometer attached near the skin overlying the heart.

The SCG is generated by micro-vibrations of the thoracic cavity caused by cardiac contractions. It consists of several peaks and troughs that encompass physiological information related to the heart, including the following features: Atrial Systole (AS), Mitral Closure (MC), Isovolumetric Contraction (IM), Aortic Opening (AO), Isovolumetric Relaxation (IC), Rapid Ejection (RE), Aortic Closure (AC), Mitral Opening (MO), and Rapid Filling (RF), as illustrated in Figure 6. Among these, the AO point is characterized by its distinct features, making it easier to detect. We have designed a method for accurately locating the AO point, which is both simple and efficient. The detailed steps of the algorithm are outlined in Algorithm 1. The algorithm was developed in the MATLAB 2021b environment. Initially, the Hilbert Transform (HT) is employed to calculate the upper and lower envelope signals of the SCG. The MATLAB function findpeaks is then utilized to detect the peaks of the upper envelope and the troughs of lower envelope, with a minimum spacing parameter set to delta = 0.6 * f s . To address potential noise factors that may lead to false cardiac cycle segments in the SCG signal, we apply a criterion based on the absolute difference between the x-coordinates of the peaks of the upper envelope and the troughs of the lower envelope. If this absolute difference exceeds delta 1 = 0.15 * f s , we discard the detection of the AO point within that cardiac cycle segment. We observed that the IM and IC points are distinctly characterized as the lowest and second-lowest points within each valid cardiac cycle segment. For each qualified trough of the lower envelope, we first identify the lowest and second-lowest points, which correspond to the IM and IC points, respectively. The AO point is always located within the interval between the IM and IC points. Therefore, we continue to search for the maximum value within the x-coordinate range defined by the IM and IC points, which represents the AO point. This method effectively prevents the erroneous detection of the MC and RE points as AO points, as well as mitigating the risk of falsely identifying peaks from noisy cardiac cycle segments as AO points.

Algorithm 1

Input : f s , SCG k , k = 1 , 2 , 3 , . . L , L = length SCG

HRV analysis metrics can be derived from the IBI sequence. In this study, we utilize the three most commonly used time-domain metrics for HRV analysis: SDNN, RMSSD, and pNN50. SDNN represents the standard deviation of all normal sinus beats’ IBI. A higher SDNN indicates greater HRV, reflecting overall autonomic nervous system activity. RMSSD is the square root of the mean squared differences of successive IBI and is primarily used to reflect vagal (parasympathetic) activity; higher RMSSD values indicate increased vagal activity. pNN50 denotes the percentage of successive IBI differences greater than 50 ms, which, similar to RMSSD, reflects vagal activity. A higher pNN50 percentage suggests a stronger regulatory influence from the vagal nerve. The time-domain statistical expressions for SDNN, RMSSD, and pNN50 in HRV analysis are as Equations 23–25: S D N N = 1 N ∑ i = 1 N   I B I i − I B I ¯ 2 (23) R M S S D = 1 N − 1 ∑ i = 2 N   I B I i − I B I i − 1 2 (24) pNN 50 = ∑ i = 2 N   IBI i − IBI i − 1 > 50 ms N (25) Where IBI ¯ represents the mean of the IBI sequence, and N represents the length of the IBI sequence.

We utilized the Texas Instruments AWR1642Boost commercial millimeter-wave radar development kit and the DCA1000EVM data acquisition board to collect reflected echo data from the human thorax. The raw IF data from the radar was transmitted to a laptop via Ethernet. Additionally, we employed the Shimmer ECG sensor to collect ECG signals synchronously as a gold standard reference. A synchronization data acquisition program for the AWR1642 mm-wave radar sensor and Shimmer ECG sensor was developed using MATLAB R2021b. AWR1642 millimeter wave radar antenna uses TDM method to set up two transmit four receive mode, with the Chirp parameters specified in Table 1. A total of 13 participants were involved in the experiment, each person collected data for 10 min. The ages of the participants ranged from 22 to 34 years, and they were dressed in regular t-shirts while lying on a bed. The radar sensor was positioned 0.6 m above the participants, with ECG electrodes attached. Upon starting the MATLAB program and setting the acquisition time, data could be collected for the predetermined duration. The configurations of the radar sensor, Shimmer ECG sensor, and the experimental setup are illustrated in Figure 7. All data processing and algorithm development were conducted within the MATLAB R2021b environment. The primary computer hardware included an Intel i5-10400 CPU, 16 GB of RAM, a 500 GB solid-state drive, and a 1 TB mechanical hard drive.