Vital Sign Detection via Angular and Range Measurements With mmWave MIMO Radars: Algorithms and Trials

Vital signs such as respiration rate and pulse (heart rhythm) rates are very useful in determining health problems and managing well-being. Heart rate measurement techniques have evolved from counting pulses to contact-based approaches such as electrocardiography (ECG) and photoplethysmography (PPG). The current availability and proliferation of low-cost and easily integrable millimeter wave (mmWave) radar sensors make contact-free simultaneous detection of weak vital sign signals highly promising. However, radar echoes reflected from the mechanical behavior of the heart (i.e. heartbeat) are complex and mixed with other (usually stronger) mechanical signals (respiration, body movements, etc.). Such echoes, as captured by radar sensors, comprise a range of Doppler frequencies with varying densities. These frequencies can be directly related to the respiratory and heart rate patterns. Conventional spectral-based approaches adopted directly from radar signal processing do not always yield an accurate estimations of these frequencies. For example, the harmonics of the respiration signal, intermodulation, and cross-product terms complicate the robust detection of heart rate signals. In this study, we show that spatially scanning a human participant using digital beamforming methods provides a more reliable estimation technique for heart rate detection. Prior studies using radar sensors that can filter participants in the range use a range selection strategy based on either the maximum return signal or the highest variance. We developed a strategy to incorporate multiple range and angular bins to reduce the probability of erroneous heart rate estimation. The algorithm, tested on multiple participants, provides an accurate heart rate estimation compared to existing methods.


I. INTRODUCTION
Accurate and longitudinal measures of heart rate and breathing rate provide insights into the physiological and mental well-being of people. Vital sign analysis can detect the onset of physiological decline, hint at life-threatening events, The associate editor coordinating the review of this manuscript and approving it for publication was Larbi Boubchir . determine underlying causes, and allow physicians to make timely interventions to improve patients' health [1], [2], [3], [4], [5], [6].
In clinical settings, electrocardiography (ECG) and photoplethysmography (PPG) are typically used to extract heartbeat rates and information regarding cardiopulmonary signals [7], [8]. However, both methods are based on contactbased sensors. The inconvenience of contact-based methods VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ and the current application-driven need for contactless monitoring have opened avenues for research on noncontact measurements of heart and respiration rates. This inconvenience stems from discomfort and skin irritation caused to participants by the straps and electrodes usually used in contact-based systems, which is especially true in the case of neonates and burn victims in hospital environments [9], [10]. Application-driven requirements pertain to the demand for ubiquitous measurements of human physiological signals, and have surpassed the medical field. These signals, which were previously associated with hospitalized patients [11], are now used in a host of home, automotive, and security applications [12], [13], [14], [15], [16], [17], [18], [19]. The non-invasive nature of radar signals makes them suitable candidates to meet the requirements of short-range applications involving human subject detection, tracking, and the measurement of cardiopulmonary signals. Vital sign detection has been studied in the literature using different radar waveforms across various operating frequencies [11], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33]. The underlying theory governing most radar technologies in estimating vital signs is to measure the phase changes in the radar transmit signal due to chest and thoracic movements. These movements occur owing to a combination of mechanoacoustic signals, which are generated by natural physiological activities in the body. This leads to a complex vital signs signal that contains micro-vibrations of the chest wall surface further complicated by the slight motions of a subject [3], [34].
Given that the human body is a dispersed target, the movement of the heart muscle inside the chest causes the entire body to move slightly during each beat. The movement of the chest surface due to the heart is very small compared with that of the heart itself. However, because the attenuation of the electric field is high in body tissues, reflections from inside the body do not contribute significantly to signal modulations [35]. Hence, the measurement of heart rate using radars can be considered similar to the BallistoCardioGram (BCG) [36]. In essence, our objective is to measure vibrations in the area around the heart. These vibrations are about a fraction of a millimeter [37], [38], [39] and can be measured with a very sensitive radar directed towards the thorax. The sensitivity requirement is related to the radar carrier frequency. In [40], the authors compared the theoretical and empirical phase variations owing to chest displacement for different frequency signals. They noted that increased sensitivity at higher frequencies reduces the effect of noise and increases the accuracy of detecting the peaks of vital sign signals. In addition, systems at higher frequencies provide a better range resolution, have relatively smaller form factors, and are thus more suitable for ubiquitous applications.
One of the biggest challenges in estimating heart rates from chest vibrations is that the chest surface motions are dominated by movements due to breathing and are at least an order of magnitude greater than those due to the heartbeat. The frequencies corresponding to respiration range between 0.1-0.8 Hz and that due to heart rate between 0.8-4 Hz. In addition, the periodic influences of breathing and heart rate are not perfect sinusoids. As a result, spectral decomposition identifies significant variations around the integer multiples of the fundamental breathing frequencies [41]. These breathing harmonics are close to, and may interfere with, the detection of the heart rate signal. Accordingly, achieving reliable estimates of the heart rate becomes challenging.
The primary concern now is where the heart rate can be reliably monitored. Vital sign studies that involve resolving subjects in the range domain traditionally strive to detect an optimal range bin based on the maximum returned power or the variance of the signal [11], [29], [42], [43], [44], [45], [46], [47]. The idea of selecting the maximum power range bin is that, because a subject is still, most other reflections from the body can be removed using clutter removal in the algorithm, and the chest motion will be left as the dominant return. However, static clutter removal does not remove all other undesired returns because there are reflections from static body parts as well as body parts moving minutely, albeit intermittently. In addition, even if the highest power range bin corresponds to returns from a point around the chest, this may not be optimal for detecting the heart signal, which is often spectrally occluded by artifacts of respiration, intermodulation products, and noise. Variance-based approaches are also better at extracting respiration rates. The highest variance range bin is the one where the chest motion is high around the nominal range, but might not contain the best heart rate signal.
Physiologically, we are not attempting to detect the physical movement of the heart but the motion of the chest due to heart cycles [48]. Intuitively, given the relatively limited number of tissues between the right ventricle of the heart and sternum, the accuracy of the radar heart rate detection system can be improved by focusing the device on the sternum [49]. However, in real-life applications of vital signs, it is difficult to focus the signal at the sternum, given a number of factors. First, there are physical differences between the subjects, such as their height and size. Second, for applications similar to a patient lying on a hospital bed, the position of the subject changes over time, and it becomes challenging to track the optimal location. Hence, a position-independent vital signs detection algorithm is advantageous. For this, we leverage the fact that the heart rate can be detected at several points on the body, given that the signal is incident on the thorax.
In the present study, we do not attempt to determine a single optimal location. We first extend the spatial localization of subjects to include the angular domain by utilizing digital beamforming on an antenna array. Although beamforming techniques have traditionally been used to separate multiple subjects [24], [50], we propose a method that uses receive beamforming to virtually scan a single subject to make heart rate detection more robust. While multiple subjects can still be resolved in the range and angle domains, virtual scanning across multiple ranges and angles can be used to identify individual heart rates. Subsequently, we build an algorithm for heart rate estimation around this. Finally, we validate the existence of heart rate locations by presenting experimental results utilizing the algorithm for a number of participants and show that heart rate signals are indeed prevalent at multiple locations in the range-angle domain.
We are also interested in addressing several questions that impair heart rate detection. Because of the existence of harmonics and artifacts of breathing, typical filtering-based approaches to isolate heart rate signals may result in incorrect selection. Although methods such as adaptive noise cancellation [51], wavelet transform [52], and adaptive harmonic cancellation [53] have been proposed for isolating heart rate, they still suffer when the heart rate and respiration harmonics are adjacent. Other iterative approaches for harmonic correction, such as those utilized in [54], assume that the heart rate changes are slow. The phase changes owing to the heartbeat pulse are different in different parts of the body [55], and the respiration amplitude also varies at different points [56]. As a result, there is an inherent difference in the signal-tointerference-noise-ratio (SINR) at different locations. This implies that the impact of artifacts is not uniform across the spatial domain. Thus, the impact of artifacts and noise can be reduced by utilizing the spectral information from multiple locations. To ensure further reliability of the estimation, the algorithm includes a mechanism to avoid selecting artifacts of frequency peaks from the respiration region as heart rate candidates.
In vital signs detection, a less addressed problem is the impact of heart rate harmonics. Given a low reference heart rate, it is possible for an algorithm to select the second harmonic of the heart rate if its amplitude is larger than the fundamental frequency component. To correct this, we utilize the harmonic product spectrum (HPS), which is traditionally used in fundamental frequency measurement and pitch detection [57], [58], [59], to verify the fundamental frequency of the heart rate if the estimated frequency is high. To the best of our knowledge, HPS has not been used in vital signs detection.
In summary, our main contributions are • Introducing beamforming as a method to virtually scan a subject spatially in a single-subject vital-sign detection scenario.
• Presenting a heart rate estimation algorithm that utilizes information from multiple points on the body, using the said beamforming technique to make detection more robust.
• Providing an empirical study to show the existence of heartbeat signals across a human subject's torso.
• Validating the algorithm by conducting experiments on multiple subjects (nine in all) to assert that it works for subjects of a diverse demographic. The results show that the proposed algorithm significantly improves the accuracy of heart rate detection and tracking.
The remainder of this paper is organized as follows. Section II provides the data structure and details of the  algorithm used in this study. Section III details the setup and equipment utilized, as well as the radar chirp and frame parameters used in this study. Section IV presents the experimental results. Section V discusses further work and concludes the paper.

II. THEORY AND METHODS
In this section, we introduce the frequency modulated continuous wave (FMCW) transmission system used in the study, the consequent data structure relevant to the system employed, and the details of the proposed algorithm. It should be noted that even if an FMCW radar system is used in this study, the concepts discussed are not limited to FMCW radars, and are generalizable to all systems that are capable of differentiating objects in the spatial and Doppler domains.

A. FMCW RADAR AND DATA STRUCTURE
The FMCW radar enables us to localize subjects (or objects) in the range, angle, and Doppler domains. To meet the sensitivity requirements for heart rate estimation and deploy the spatial resolution of a subject in both the range and angle domains, we utilize a 77 GHz time division multiplexed (TDM) multiple-input-multiple-output (MIMO) FMCW radar sensor. The details of the radar used are discussed in Section III-B. Utilizing two transmit (Tx) and four receive (Rx) antenna elements, we obtain a virtual linear array of eight elements. The physical and virtual antenna structure of the sensor utilized in this study are shown in Fig. 1. The transmission scheme used in this study is shown in Fig. 2. In a single TDM MIMO frame or frame repetition interval (T FRI ), chirps are transmitted sequentially from each VOLUME 10, 2022 Tx element according to their natural spatial order. The value of T FRI , which relates to the slow time sampling rate, is chosen to satisfy the Nyquist criterion for vital sign frequencies.
When the delayed and attenuated versions of these transmitted chirps are received after reflections, they are mixed with their respective transmitted chirps to obtain an intermediate frequency (IF) signal. After analog-to-digital (ADC) sampling, we combine the discrete IF signal obtained from two chirps to generate 2D data, which contains fast time samples from one frame across the virtual channels. By stacking such 2D data collected over multiple frames (based on the observation time window), we represent the IF signal by a data cube of dimensions P × N × M , where P represents the number of virtual channels that observe the IF signal, N represents the number of fast time samples collected per chirp by the ADC, and M represents the number of slow time samples. The data cube is shown in Fig. 3. A detailed explanation of the FMCW signal and TDM MIMO transmission schemes can be found in [60], [61], [62], and [63].

B. VITAL SIGNS ESTIMATION ALGORITHMS VIA ANGULAR AND RANGE MEASUREMENTS
The proposed algorithm is divided into two parts. The first part localizes the subject in the range and angle domains. The second part utilizes the slow time signals along the localized range-angle bin pairs and processes them to obtain the vital signs. A flowchart of the first part of the algorithm is shown in Fig. 4a. For a continuous data stream, it is necessary to fix the observation time window to collect sufficient frames required for evaluation. In our study, we use a window of 25 s. A larger observation window corresponds to better Doppler resolution, as the corresponding signal contains more cycles of breathing and heart rate. Beginning with the IF signal, this step is referred to as the Time Gating block in the flowchart. Static clutter suppression is performed by the mean removal of the IF signal across both slow and fast time domains. This ensures that the reflections from the static clutter in the environment are suppressed, and the resulting IF signal contains returns from moving targets or human subjects.
Information regarding range and angle is then obtained by taking fast Fourier transforms (FFTs) in the range and array dimensions of the clutter removed data cube. These FFTs are referred to as the range FFT, and angle FFT, respectively. The size of the data cube along these two dimensions are reduced by the range and angle gating steps, depending on the size of the FFTs and the width of the range and angles that we wish to utilize. The number of slow-time samples remains the same. A holistic view of this idea is to locate a central region across the human radar cross-section and evaluate the signals across points around that center. This is shown in Fig. 4b and  4c. The radar signals first illuminate the entire environment, including the static clutter. After static clutter removal and range and angle gate processes, we restrict our evaluation to echoes emanating from the human subject.
The second part of the algorithm (after the subject localization) is summarized in Fig. 5. The phase information extraction block in the chain involves computing the phase for each slow time instant, unwrapping the phase to account for phase jumps, and computing the differential phase. Taking the differential enhances the heart rate signature (slightly higher frequency) relative to the respiration rate signal. The Doppler FFT is then conducted for each point. Starting from the 3D matrix A [u, v, m], which represents the data cube after FFTs in the range and angle dimensions, and also illustrated in Fig. 4a, the vital signs estimation part of the algorithm is as follows: 1) Select the range bin (v o ) corresponding to the highest power. If v denotes the indices of range bins, then Here, N a represents the size of the FFT along the array dimension. The index v o can be translated to physical distance d s (in meters) between the sensor and subject as where f f is the ADC sampling rate, c is the speed of light, K is the frequency deviation slope of the chirp signal, and N r represents the size of the FFT along the range dimension. This is the distance d s , as shown in Fig. 4c. 2) Select X 1 range bins centered around v o . The lowest index range bin is denoted by v 1 and the highest index range bin is denoted by v 2 . 3) Select the angle bin corresponding to the highest power.
If u o denotes the angle bin, then ) Select X 2 angle bins centered around u o . The lowest index angle bin is denoted by u 1 and the highest index angle bin is denoted by   Fig. 4a. The phase of the slow time signal contains information regarding the motion of the body due to breathing and heart rate. The heart rate region corresponds to selecting regions in the spectrum between 0.8 Hz and 2 Hz, where we expect the resting heart rate to typically be. Analyzing spectral peaks in different parts of the torso, a Pseudo Spectrum is formed and utilized to assign the heart rate and respiration rate.
restricted in range and array dimensions and is denoted byÂ[u, v, m]. 5) For each of the X 1 × X 2 range-angle bin pairs, extract the slow time signal. 6) For each slow time signal, extract the phase, and then unwrap it to incorporate phase jumps greater than ±π. 7) Compute the differential of the unwrapped phase signal.
8) Extract the Doppler spectrum by taking an FFT in the third dimension (slow time). 9) The heart rate estimate is computed for each of these bins by selecting the location of the maximum in the Doppler spectrum at the heart rate frequencies (48-120 bpm). Thus, we have a matrix H (of size X 1 × X 2 ), of heart rate estimates for each range-angle bin pair. Similarly, the respiration rate estimate is computed by considering the frequency region corresponding to the breathing frequencies (6-48 breaths per minute). Hence, we have another matrix R (of size X 1 × X 2 ) that contain the respiration rate estimates. 10) For each heart rate estimate in H, check if it is a multiple of the breathing rate (in R) and disregard those estimates. This corrects for the harmonics of respiration that may appear at the heart rate frequencies. 11) Count the number of occurrences of unique heart rate values in H. A visual representation with the estimated frequency on the horizontal axis and the number of occurrences on the vertical axis is referred to as the pseudo spectrum. 12) Similarly, count the number of occurrences of unique breathing rate values in R. 13) Assign the value that occurs most frequently in the respiration pseudo spectrum as the breathing rate. 14) Extract the value that occurs most frequently in the heart rate pseudo spectrum. If this estimate is less than 100 bpm, assign the value to the heart rate. Otherwise, generate a pseudo spectrum from the harmonic product spectrum (HPS) and assign the highest occurring estimate as the heart rate. The HPS ensures that the fundamental frequency of heart rate is selected. 15) Move the slow time window forward by the desired step duration and repeat steps 1-14. The number of range and angle bins selected, that is X 1 and X 2 , can be empirically set because, at greater distances, the number of angle bins that contain the subject is lower, and at lower distances, more angle bins are required to cover the entire subject. As a rule of thumb, the number of angle bins around the main peak in the angle FFT before the first zero is selected, and range bins within 15 cm around the main peak of the range FFT are selected. Thus, we can effectively scan the subject across the entire cross-section. In addition, further details of the pseudo spectrum and examples are illustrated in the results section IV-B. The details of the harmonic product spectrum (HPS) and its adaptation for verifying the fundamental frequency of the heart rate are provided in the Appendix.
The angular indices utilized here can be converted into the physical angle of arrival of the signals [64]. Assume that θ o , related to the angle index u o in the algorithm, represents the physical angle of the subject and corresponds to the maximum return signal. θ s , which is illustrated in Fig. 4c, represents the maximum deviation of the angle considered in the algorithm with respect to the angular location of the subject (θ o ), and is given by where θ 1 and θ 2 represent the physical angles corresponding to the lowest and highest angle indices (u 1 and u 2 ) considered in the algorithm.
In this approach, the participant is assumed to be within the same range bin during the slow time window. The step size of the new estimate can be set as desired. To account for the fact that the participants' range and angle locations may change over time, a recalculation and update of the nominal range and angle bin is performed at each sliding window. In our experiments, this update was performed every second.

A. PARTICIPANTS
The experiments were conducted in accordance with the approval and guidelines of the Institutional Review Board of the University of Texas at Dallas. Participants were a mixture of male and female university students. All participants completed a consent form that contained information regarding the experiments and provided consent to share the measured data anonymously.

B. DEVICES
The Texas Instruments (TI) AWR1843 BOOST [65] model was used for the vital signs experiments. Two of the three available transmit antennas at the same elevation level were chosen to obtain spatial scanning in the azimuth. All the four receiving channels were used. The radar data were transferred to a PC using a TI DCA1000EVM data capture board [66], and further processing was performed in MATLAB. Additional specifications and details are provided in [65]. Texas Instruments' AWR family of transceivers has been used for automotive startups [67] as well as vital signs detection research [43].
To validate the results obtained from the radar sensor, a commercially available chest strap, the Polar H10, was used as a reference sensor [68]. Polar H10 is considered a reliable sensor for measuring heart rate and is used as the gold standard for comparison with other wearable devices [69]. Both the AWR1843 sensor and wearable chest strap are shown in Fig. 6.

C. RADAR PARAMETERS FOR EXPERIMENTS
The transmitted chirp parameters are listed in Table 1. Based on the parameters listed, the maximum detectable range (R max ) is approximately 6 m (R max = c × f f 2K ). To extract motion or Doppler information, the environment was scanned continuously at a slow time sampling rate of 20 Hz.

D. PROCEDURE
The experiments were conducted in rooms containing regular clutter. The radar sensor was placed on a plane table, and a chair was placed approximately 1.3 m from the sensor. The measurement setup is illustrated in Fig. 7. The test participants were asked to wear the wearable chest strap to record the reference heart rate. The participants sat on a chair and remained still for 80 s per measurement duration. This duration was repeated between 4 and 10 times. The 80 s duration was chosen so that the participants could sit relatively still and did not move involuntarily due to tiring. The total duration of measurement for each participant is presented in Table 2. The sensor was placed at approximately the level of the chest.

A. EMPIRICAL SUBJECT LOCALIZATION AND VITAL SIGNS SIGNAL EXTRACTION
Subject localization in the range domain is shown in Fig. 8a. After static clutter removal, ideally, only echoes from nonstatic clutter, including the participant, are observed. The participant is indicated by an arrow in the range profile. This range information, however, may change with slow time owing to random motions of the subject, even in the case where one attempts to stay still. The range variation across slow time with minimal motion is shown in Fig. 8b. After the angle FFT, the participant can be localized in the range and angle domains, as shown in the heat map image in Fig. 8c.
We can observe that the human RCS covers multiple range and angle bins. After phase extraction and unwrapping, the signals can be filtered according to the frequencies generally associated with the respiration and heart rate. The filtered time-domain signals representing respiration and heart rate are shown in Fig. 9. They can be directly mapped to chest displacements; however, exact displacement measures are not extracted here, with frequency estimation being the goal.
After the differential phase of the slow time signal corresponding to a particular range-angle bin pair is converted to the frequency domain (without filtering), the Doppler spectrum ideally contains peaks for both respiration and heart rate. A typical Doppler spectrum with clear peaks at the respiration and heart rate frequencies is shown in Fig. 10a. The corresponding peaks (at the reference values) for breathing and heart rate are represented by a red diamond and a green triangular arrow, respectively.
However, in many scenarios, the heart rate is buried underneath other spurious peaks. For example, the spectrum in Fig. 10b shows respiration, its harmonics, and other spurs that overwhelm the heart rate peak. The diamond marker again represents the respiration peak. The dotted line in the middle indicates the reference heart rate location. In this case, it is difficult to distinguish the heart rate peak from the other peaks. In such a scenario, even if visual inspection and knowledge of the reference heart rate might allow us to deduce the existence of the heart rate peak, deriving it through a peak selection algorithm becomes difficult.

B. MULTIPLE POINTS PSEUDO SPECTRUM GENERATION AND TRACKING OF HEART RATE
In the algorithm, Doppler information dispersed across multiple range-angle bin pairs was extracted. From the Doppler FFT derived for each range-angle bin pair, we selected the maximum peak for the heart rate signal that lies within 48-120 bpm. Simultaneously, the maximum peak in the respiration region was selected.
Then, a histogram was formed to count the number of peaks represented in those range-angle bin combinations, and a pseudo spectrum was formed. This pseudo spectrum represents the count of the maximum-amplitude frequency component. An example pseudo spectrum for respiration and heart rate is shown in Fig. 11. The peak in the heart rate region corresponds to a reference heart rate of 80 bpm. Similarly, the pseudo spectrum for the respiration region shows a peak at 14 breaths per minute. The respiration peak is also enhanced, with tiny peaks appearing at the harmonic frequencies of the breathing rate. To reduce the impact of respiration harmonics, all range-angle bin combinations in which the peak in the heart rate region is a multiple of the peak in the VOLUME 10, 2022  respiration region, are not included in the generation of the pseudo spectrum. The effect of this harmonic correction step is illustrated in Fig. 12. The incidence of the heart rate remained the same. However, the respiration harmonic was substantially reduced.
In some instances, it so happens that the highest spectral component in the heart rate region is a harmonic of the heart rate itself. Hence, if the estimated heart rate is greater than 100 bpm (greater than twice the lowest heart rate frequency), we use the HPS to validate whether the peak is indeed the fundamental or its second harmonic. An example is shown in Fig. 13. The initial Doppler spectrum, shown in Fig. 13a, has a larger peak at the second harmonic of the heart rate, as opposed to the fundamental frequency. The HPS generated by taking the product of the decimated spectra is shown in Fig. 13b, where the highest peak is at the fundamental frequency (denoted by f h ).
These pseudo spectra were generated for each sliding window, and the peak was tracked to obtain heart rate tracking over time. The heart rate tracking of one measurement dataset for nine different participants sitting in front of a radar, obtained using a measurement window of 25 s and shifts of 1 s, is shown in Fig. 14. This is the heart rate tracking for a period of 80 s and is shown in conjunction with the reference heart rate values in beats per minute obtained from the wearable chest strap. The use of the pseudo spectrum, generated by combining spectral results across different spatial locations, greatly reduces the number of spurs and enhances the strongest spectral components. Hence, the proposed algorithm correctly tracked the reference heart rate.

C. LOCALIZING HEART SIGNAL
To empirically validate the existence of heart rate spectral peaks across multiple range and angle locations, the information generated from the algorithm (Section II-B) was used in conjunction with the reference values obtained from the wearable sensor.
The estimation algorithm uses spectral information from multiple range and angle bins, with each bin proposing a heart rate candidate based on the Doppler spectrum at that point. This method can be used to form a 2-D matrix with heart rate candidates. H is the matrix of the heart rate estimates for each range-angle bin pair. Let H g be a matrix of reference heart rate values and be of the same size as H. Each cell in H g contains the same reference value. We used these two matrices to assign a correct detection, as opposed to an incorrect detection, as follows: 1) Subtract the two matrices -M true = H g − H 2) For each element in M true , where the absolute difference between the estimated and ground value is greater than a certain threshold, set that value to a 0. This represents locations of incorrect estimations. In our study, the threshold is set to 3 bpm. 3) For all range-angle indices, where the difference is lower than the threshold, set the corresponding value to a 1. 4) Store the binary matrix formed for that particular measurement block.
We obtained the binary matrices for multiple measurement blocks over time. If we sum these matrices, we obtain a tentative estimate of the spectral location of the correct heart rate values. To validate this, we used data from two participants. For the first participant, approximately 500 estimations were summed, whereas for the second participant, approximately 800 were summed. The resulting heat map of the spatial location of correct detection is shown in Fig. 15. The heat map intensity represents the number of instances in which a  particular range-angle bin pair provided a correct estimate of the heart rate.
For both participants, there is a dispersed region in the range and angle domains, where we can expect to find the heart rate signal. For any single measurement instance, this location will change. This is because the precise return signal depends on various factors including the orientation, distance, and size of the participant. Hence, it is not ideal to attempt to find the perfect location to estimate the heart rate of a participant. However, across multiple measurements, there is a region represented by the bright spots in the figures, where we can expect good estimates to be located. The range locations were approximately 30 cm to 40 cm across. The angular locations were much wider and restricted by the angular resolution of the system as well as the distance of the sensor and subject. We can expect this to become narrower when the subject is farther away and in the case of a larger MIMO antenna structure, rather than the eight-channel linear array used in this study.

D. PERFORMANCE ANALYSIS
The mean absolute errors and correct detection percentages for each participant are listed in Table 3. The accuracy of the proposed algorithm was high across all participants, ranging from approximately 84% to 97%. Overall, 4320 heart rate estimates were obtained. Among all estimates and across all participants, the mean absolute error was 1.31 bpm and the detection accuracy was 93.54%. VOLUME 10, 2022 FIGURE 12. Effect of the harmonic correction process on the resulting pseudo spectrum. (a) shows the respiration peak. In (b), the third harmonic of the respiration rate is initially seen as the highest peak in the heart region pseudo spectrum (diamond marker on the dotted line). After harmonic correction, the harmonic peak is removed, and the heart rate is correctly selected (triangle marker on the solid line). Several metrics have been employed in the literature to characterize the performance of vital signs estimation algorithms. The values shown in Table 3 and the corresponding accuracy numbers relate to accuracy being defined as the proportion of instances where the radar estimate of the heart rate is within 3 bpm of the reference sensor. The authors in [70] achieved an overall accuracy of 55% when they allowed for the estimated values to lie within 5% of the reference device and 60% for deviations of up to 10%. The proposed algorithm achieves accuracies of 93.7% and 97.3% when deviations of 5% and 10% are allowed, respectively. In either case, the achieved accuracy is significantly higher.
The relationship between the heart rate estimates from the radar sensor and the reference values is shown in the scatter plot in Fig. 16a. A reference correlation line that represents 100% correlation was drawn through the plot, and all measurement instances of all participants were included.
A high Pearson correlation coefficient (R = 94.4%) was observed. In terms of correlation, the authors in [43] report an 80% correlation for heart rate. This was for a single measurement scenario with a single participant in a supine position. The proposed algorithm performed better across multiple participants and measurements in a seated position.
However, a correlation plot may not be the best method to justify the agreement between different measurements [71]. When comparing two measurement methods, the Bland-Altman plot may be better suited to identify an interval of agreement within which a proposed method, compared to the reference method, falls. The Bland-Altman plot corresponding to the measurements is shown in Fig. 16b. The horizontal axis represents the mean of the heart rate estimates as measured by the radar and the reference sensor. The vertical axis represents the difference between the two measurements. The mean of the differences is denoted by the solid line. This mean difference of 0.46 beats per minute shows a very small bias. The dotted lines on either side of the mean represent the 95% confidence interval and hence the limits of agreement (LoA) at −6.3 bpm and 7.2 bpm. Fig. 16c shows the Bland-Altman plot with the differences (vertical axis) represented as percentages of the mean of the reference and radar estimates. Finally, the histogram in Fig. 16d shows the proportion of deviation of the radar estimates from reference measurements. The bin width is 1 bpm and approximately 96% of the estimates are within 5 bpm or less of the reference measurements.
The authors in [42] devise a range bin selection method by utilizing the coherency between the magnitude and phase in each range bin. They used data from one channel of an FMCW radar to estimate heart rates and reported 70.7% accuracy by allowing a deviation of 5 bpm from the reference values. Utilizing multiple channels in our proposed algorithm, as well as information from multiple range-angle bins, we obtained 95.9% accuracy by allowing the same deviation from the reference values. This represents a significant improvement compared to the latter study. A comparison of the accuracy measures from the literature and the proposed algorithm is summarized in Table 4. Overall, the accuracy of the proposed algorithm is significantly higher than that of conventional methods.
Worth bearing in mind is that the proposed algorithm can estimate vital signs based on information from a single range-angle bin pair. This will reduce the computational complexity of the algorithm but will come at a cost of reduced accuracy. Firstly, we define complexity in terms of only the number of complex multiplications and arctangent operations involved, given their dominance over other operations. As defined in Section II-B, X 1 and X 2 denote the number of range and angle bins to be utilized respectively. The number VOLUME 10, 2022 FIGURE 15. Heart rate correct detection heat map across the range and angle measurements for two participants. of complex multiplication operations (C CM ) involved is then C CM = PMN r log 2 (N r ) + X 1 MN a log 2 (N a ) + X 2 X 1 N d log 2 (N d ), (1) where P is the number of virtual channels, and M is the number of slow time samples. N r , N a and N d are the sizes of range, angle, and Doppler FFTs respectively. The number of arctangent operations (C AT ) is When only a single range-angle bin pair is utilized, then X 1 = X 2 = 1. Using multiple range-angle bins, the number of arctangent operations increases by a factor of X 1 X 2 . However, the number of complex multiplications is controlled by the range gating and angle gating steps. Based on the parameters utilized in this study, there are approximately 2.5 times more complex multiplications while using the proposed algorithm (X 1 X 2 range-angle bin pairs) as opposed to using a single range-angle bin pair. However, the improvement in accuracy is significant. Fig. 17 shows the Bland-Altman plot for the same dataset, generated by estimating the heart rate using a single range-angle bin pair. Unlike the proposed method, the limits of agreement (LoA) were larger. Only approximately 69% of the radar estimates were within 5 bpm of the reference sensor. This illustrates a significant improvement in accuracy when employing the proposed method.

V. CONCLUSION
In this study, we utilized an 8-channel MIMO FMCW radar sensor to extract the heart rates of the human participants. Our novel approach differs from the notion of treating the human body as a small reflector in space, and hence, attempting to find an optimal location for the heart rate signal. An algorithm for dividing the space into range-angle bins and utilizing multiple numbers of them to extract the heart rate was proposed. A harmonic correction method to account for respiration harmonics and utilizing the harmonic product spectrum (HPS) to correct the spectral influence of heart rate harmonics was included. The fact that heart rate signals are found at multiple locations in the body was empirically validated by showing heat maps of the correct heart rate values across the range and angle domains. The proposed algorithm advances vital signs detection and tracking, demonstrating superior performance compared with previous research. Testing on multiple participants yielded results with a low mean absolute error of 1.31 bpm and a high accuracy of 93.5% with respect to the heart rate. In addition, the proposed method can be used in conjunction with other signal processing techniques to enhance existing methods.
It should be noted that while the antennas utilized in this study were aligned in the azimuth, extending this work to include elevation bins to create a 3-D search algorithm could yield interesting results. The proposed algorithm can be extended to include the elevation domain as well. Such an approach would require a 2-D antenna array with a sufficient number of channels along the elevation domain, thereby improving the angular resolution in that domain. Other methods to improve angular resolution are to utilize either a synthetic aperture radar approach or a cascaded radar board. This represents a potential avenue for future research. With a lower angular resolution, improvements in the estimation results may not be as significant.
A further challenge in utilizing radar sensors as the tool of choice for vital signs measurement applications is the problem of motion. Robust detection of vital signs in the presence of motion will be pursued in the future. While designing applications such as sleep monitoring or monitoring   the vital signs of patients in hospital environments, the subjects will not be in constant motion. The occurrence of motion, other than that of the torso, is intermittent. As such, cases in which the vital signs signal is obscured by motion can be addressed by either removing the duration of motion corruption or by attempting to detect vital signs signal from a motion-corrupted signal. The ideas presented herein can be used in conjunction with any motion correction algorithm and will enhance the detection and tracking of heart rates, regardless of the presence or absence of movement.

APPENDIX HARMONIC PRODUCT SPECTRUM (HPS)
The harmonic product spectrum (HPS) is used to detect the fundamental frequency of a periodic signal [57], [58], [59]. The HPS is generated by multiplying the decimations of the original magnitude spectrum and the peak in the HPS is determined as the fundamental frequency. Provided that harmonics exist at integer multiples of the fundamental frequency, HPS enhances the fundamental frequency component. Given a discrete Fourier transform (DFT) represented by Y (f ), the HPS is obtained as where Q is the number of harmonics to be considered and determines the number of decimations of the original VOLUME 10, 2022 FIGURE 18. The harmonic product spectrum (HPS) generation method. The amplitude spectrum is decimated up to a certain rate and then multiplied. The decimation is denoted by /q above, with q representing the decimation rate (1 to Q).

FIGURE 19.
From the top -The decimated amplitude spectra for decimation rates of 1, 2, 3 and the harmonic product spectrum (bottom). The fundamental frequency component is enhanced relative to harmonic frequencies.
amplitude spectrum. The decimated amplitude spectra for a decimation rate q is denoted by |Y (qf )|. The fundamental frequency (f 1 ) is then determined by In the application of HPS to our vital signs spectra, because the respiration amplitudes are much larger than the heart rate amplitudes, the peak search in the HPS is performed at the heart rate frequencies. That is, if the fundamental frequency of the heart rate is denoted by f h , then: The HPS generation method is illustrated in Fig. 18. Spectral representations of a simulated sinusoid (with a fundamental frequency of 1 Hz) and it's harmonics are shown in Fig. 19. It also shows the decimated spectra and the corresponding HPS obtained from the product. With each decimation, the harmonic component aligns with the fundamental frequency, and their product enhances the amplitude of the fundamental frequency component. The choice of the number of harmonics to be considered is dependent on the system sampling rate and the incidence of harmonics. In our evaluations, we considered a maximum decimation rate Q of 3.