Signal Separation Using a Mathematical Model of Physiological Signals for the Measurement of Heart Pulse Wave Propagation With Array Radar
aa r X i v : . [ ee ss . SP ] S e p Signal Separation Using a Mathematical Model ofPhysiological Signals for the Measurement of HeartPulse Wave Propagation With Array Radar
Takuya Sakamoto
Abstract —The arterial pulse wave, which propagatesalong the artery, is an important indicator of variouscardiovascular diseases. By measuring the displacementat multiple parts of the human body, pulse wave ve-locity can be estimated from the pulse transit time.This paper proposes a technique for signal separationusing an antenna array, so that pulse wave propagationcan be measured in a non-contact manner. The bodydisplacements due to the pulse wave at different bodyparts are highly correlated, and cannot be accuratelyseparated using techniques that assume independent oruncorrelated signals. The proposed method formulatesthe signal separation as an optimization problem, basedon a mathematical model of the arterial pulse wave. Theobjective function in the optimization comprises fourterms that are derived based on a small-displacementapproximation, unimodal impulse response approxima-tion, and a causality condition. The optimization pro-cess was implemented using a genetic algorithm. Theeffectiveness of the proposed method is demonstratedthrough numerical simulations and experiments.
Index Terms —Array radar, mathematical model,pulse wave velocity, signal separation
I. Introduction A CCORDING to reports from the National Center forHealth Statistics of the Centers for Disease Controland Prevention in the United States, the number of adultsdiagnosed with heart disease was 30.3 million in 2018[1]. That accounts for 12.1% of all adults, resulting inmore than 600 000 deaths/year, making heart disease theleading cause of death in the United States. Constantmonitoring of health status is important for the preventionand treatment of heart disease. Some signs of heart diseasecan be detected from pulse wave propagation along thearteries caused by the pulsation of the heart. In particular,pulse wave velocity (PWV), which is the velocity of thepulse wave, is related to blood pressure and vascularstiffness, and is important for early detection of signs ofheart disease such as hypertension, arteriosclerosis, andmyocardial infarction. Conventionally, a cuff-type contactsensor has been used to measure pulse waves. For example,
This study was supported in part by JSPS KAKENHI 19H02155,JST PRESTO JPMJPR1873, and JST COI JPMJCE1307.T. Sakamoto is with the Department of Electrical Engineering,Graduate School of Engineering, Kyoto University, Kyoto, Kyoto615-8510, Japan. [email protected]. Sakamoto is also with Japan Science and Technology Agency,PRESTO, Kawaguchi, Saitama 332-0012, Japan.
PWV can be calculated by wearing multiple cuffs on thefour limbs, measuring the pulse transit time (PTT) be-tween the upper arm and ankle, and dividing the distancebetween these parts by the PTT. However, in contact-typemeasurements, it is necessary to attach multiple sensorsto multiple body parts at the same time. That makes themeasurement inconvenient and time-consuming, resultingin discomfort and a sense of restraint, which makes thiscontact measurement unsuitable for long-term, continuousmonitoring.Instead of such contacting sensors, if non-contactingradar-based measurement is introduced, PWV can becontinuously measured for a long time. We review existingradar-based measurement of pulse waves. Buxi et al. [2]–[4] attached a 1-GHz continuous wave (CW) radar antennato the neck/chest and also attached an impedance cardio-graphy sensor on the waist/shoulders to measure the PTTfrom the time difference between the signals. Ebrahim etal. [5], [6] attached a 900-MHz CW radar antenna to thechest, and also attached a photoelectric plethysmography(PPG) sensor to the ear; they measured the PTT from therise time difference between the radar and PPG signals.Kuwahara et al. [7] measured PTT and PWV from thepeak time difference between the 2.4-GHz CW radar anda piezoelectric pulse wave sensor attached to the finger.In these studies, radar systems and contact sensors wereused together, to measure pulse wave propagation.Next, we discuss the measurement of pulse waves usingonly radar, without other sensors. Lauteslager et al. [8]used ultra-wideband (UWB) radar with a center frequencyof 3.8 GHz and a bandwidth of 1.0 GHz; they placed asingle antenna on six body parts sequentially, assumingstationarity of pulse wave propagation during the measure-ment. Tao et al. [9] used two wristwatch-type impulse radioUWB radars with a bandwidth of 4 GHz. They attachedthem to the arm and leg and calculated the PTT andPWV from the delay time of the radar signals. Tang etal. [10] used a wrist-worn device with a CW radar and aself-injection-locked radar. The user is supposed to keepthe device in front of their chest to measure the pulsewave at the chest and wrist simultaneously. In that study,the displacement of the chest was measured using non-contacting radar, but the displacement of the wrist wasmeasured using contact-type radar. None of the abovethree studies reported non-contact measurements of pulsewaves because they attached radar antennas to the body surface.Lu et al. [11] measured the chest and calf simultane-ously using two C-band radars and calculated the PTT.Vasireddy et al. [12] installed two 24-GHz radar systems15-cm away from the chest and legs, and calculated thePTT. Although they installed radar antennas withouttouching the body, realizing a non-contact measurement,the installation of antennas near multiple body partsmeans that the positions of the antennas must be adjusteddepending on the position and posture of the person.By contrast, Michler et al. [13] used phased-array radarat a distance of 10 cm from the abdomen, to performsimultaneous measurements at two abdominal points 7.3cm apart, using a single radar system. Because the phased-array radar transmits signals in a specific direction, setin advance, the position and posture of the person beingevaluated must be restricted, which is also unsuitable forlong-term measurements.In this research, non-contacting measurement of pulsewave propagation is performed using a single radar systemwith array antennas, without imposing restrictions on theposition and orientation of the person. Non-contact andunrestrained pulse wave measurement is realized by usingarray signal processing, to form appropriate beam patternsaccording to the position and posture of the person beingtested. The main problem is that the body displacementsdue to the pulse wave in multiple parts of the bodyare highly correlated, although there is a time differencecorresponding to the PTT. For this reason, techniquesfor non-correlated signals are not suitable. In the presentresearch, we propose a method to separate the radarechoes from multiple body parts using a mathematicalmodel of physiological signals based on prior knowledge. Inparticular, the skin displacement signals due to the pulsewave are modeled to formulate the measurement of PTT asan optimization problem. The proposed method is appliedto simulated and experimental data to demonstrate itsperformance.In this paper, we use the following notations: Lowercase a for scalars, uppercase A for matrices, and boldfacelowercase a for vectors. The complex conjugate is denotedby a superscript asterisk a ∗ . ℜ and ℑ denote the realand imaginary parts of a complex number, respectively. ∠ a denotes the phase of a complex value a . C m denotesan m -dimensional complex coordinate space. The matrixtranspose is denoted by a superscript T, as in A T ; thematrix conjugate transpose is denoted by a superscriptH, as in A H . The inverse of matrix A is denoted as A − .Symbol ◦ denotes the element-wise multiplication opera-tor. diag {· · · } denotes a diagonal matrix whose entries areall zero except for its main diagonal, where the diagonalentries starting in the upper left corner are given in theargument.The definitions of ‘uncorrelated’ and ‘independent’ areas follows. Let f ( t ) and g ( t ) be stationary ergodic signals.If f ( t ) and g ( t ) are uncorrelated, R ∞−∞ f ∗ ( t ) g ( t )d t = 0holds. Otherwise, the signals are correlated. If f ( t ) and g ( t ) are independent, p f ( f ) p g ( g ) = p fg ( f, g ) holds, where p f ( f ) and p g ( g ) are the probability distribution functionsof f and g , and p fg ( f, g ) is the joint probability distri-bution function of f and g . Otherwise, the signals aredependent. II. Measurement of Pulse Wave Propagation
The typical PWV in healthy young people is 3 to 5m/s, and if it exceeds 10 m/s, treatment is required.In a simplified model, PWV is represented as v PWV = p βP/ ρ , where β is a stiffness parameter, P is arterialpressure, and ρ is the blood density. Because this formulaincludes β and P , which are related to arteriosclerosisand hypertension, PWV is an important index that re-flects various signs of cardiovascular disease. Widely usedmeasures include carotid-femoral PWV, which measuresthe carotid and femoral arteries, and brachial-ankle PWV(baPWV), which measures the brachial and ankle arteries.Because the baPWV is measured by wearing cuffs on thearms and legs, restraining the subject during measurementis unavoidable. To perform non-contact measurement us-ing a radar system instead of the conventional contacttype measurement, the displacement of the skin surfacecaused by pulse wave propagation is measured. If theskin displacement at two body parts on the pulse wavepropagation path can be measured at the same time, thePTT can be measured from the time difference betweenthe pulse wave signals. The PWV is obtained simply bydividing the distance between the two body parts by thePTT. As described in the previous section, in the past,the antenna position had to be adjusted according to theposition and orientation of the person being tested, sothe method was non-contacting but not unconstrained.In this paper, a non-contacting and non-constraint pulsewave propagation measurement is realized, by combiningan array radar system and array signal processing. III. Pulse Wave Measurement Using Radar
A. Signal Model in Radar Measurement of the Pulse Wave
A radar system with an antenna array is used to mea-sure physiological signals. In particular, we use a multiple-input multiple-output (MIMO) array comprising M and M elements for transmitting and receiving, respectively.Following [14]–[16], we assume that the distance fromthe array to the targets is much greater than the arrayaperture; array elements are closely spaced so that thedirection of arrival of each echo can be approximated asthe same for all elements. In addition, we also assumethat mutual coupling between array elements is negligible.Under these conditions, the system is modeled using avirtual array of M = M M elements. In particular, weassume that the virtual array is a uniform linear arraywith a spacing of λ/
2, where λ is the wavelength.Let us assume that the number of body parts (hereaftercalled "targets") contributing to the reflection is N , andthat N ≤ M is satisfied. The displacement of the j -thtarget is d j ( t ) as a function of time t . Strictly speaking, d j ( t ) is not the actual displacement, but the line-of-sight component of the displacement. The skin displacementvector is denoted by d ( t ) = [ d ( t ) , d ( t ) , · · · , d N ] T . Weassume that the transmitted signal is approximated as anarrow-band signal and that the time delay correspond-ing to the propagation path length can be expressed asa phase shift. The echoes are phase-modulated by thedisplacement as s j ( t ) ∝ e j2 kd j ( t ) , where k = 2 π/λ isthe wave number. The echo vector is denoted by s ( t ) =[ s ( t ) , s ( t ) , · · · , s N ( t )] T . Let the propagation channel be-tween the i -th element and the j -th target be a i,j , andthe propagation channel matrix be an M × N matrix A = ( a i,j ). We assume a stationary propagation path; A does not change during the measurement. The signal x i ( t )is received at the i -th element, which forms a signal vector x ( t ) = [ x ( t ) , x ( t ) , · · · , x M ( t )] T . Then x ( t ) is expressedas x ( t ) = A s ( t ) + n ( t ) , (1)where n ( t ) is an additive noise component. The propa-gation channel matrix A is also called a ’mixing matrix.’For example, when assuming free-space propagation andtargets located in the far-field, propagation channel a i,j can be approximated as a i,j ≃ ξ ( θ j )e j2 kl j e j ku i sin θ j /l j ,where ξ ( θ j ) is a complex coefficient that depends on thedirection of arrival θ j , and l j is the distance between the j -th target and the center of the array. In addition, u i isthe x coordinate value of the i -th element, where the arraybaseline is on the x -axis.The purpose of this study is to find an ’unmixing matrix’ W , with which we can obtain an estimate of the echo ˆ s ( t )as ˆ s ( t ) = W x ( t ) (2)and an estimate of the displacement ˆ d ( t ) asˆ d ( t ) = 12 k ∠ ˆ s ( t ) . (3)For example, in a noiseless case with A being nonsin-gular and square ( M = N ), we can set W = A − sothat the signal s ( t ) = A − x ( t ) can be completely restored,including the order and amplitude of the signal. It is notnecessary to obtain s ( t ) itself; ambiguity of permutationand constant multiplication is allowed. In other words, if α i e j kd i ( t ) is obtained with an arbitrary complex constant α i , important parameters such as PTT and PWV areestimated. For simplicity, the norm of each echo s j ( t ) isassumed to be 1, i.e., | s j ( t ) | = 1. In addition, displacement d i ( t ) is defined so that ∠ s j ( t ) = kd i ( t ), which meansthat the phase is zero when there is no displacement. Ifthese conditions are not satisfied, the complex coefficientis incorporated into the mixing matrix. For example, if theecho vector s ′ is not normalized, the mixing matrix A ′ canbe replaced by A as x ( t ) = A ′ s ′ ( t ) + n ( t )= A ′ diag {| s ′ ( t ) | , · · · , | s ′ N ( t ) |} s ( t ) + n ( t )= A s ( t ) + n ( t ) , (4) where A = A ′ diag {| s ′ ( t ) | , · · · , | s ′ N ( t ) |} and s j ( t ) = s ′ j ( t ) / | s ′ j ( t ) | . Then, the echo can be expressed using onlyits phase as s j ( t ) = e j2 kd j ( t ) .Note that the received signal includes not only echoesfrom the target person, but also clutter from surroundingobjects (e.g., the floor and bed). Because the objects otherthan the target person are assumed to be stationary,the Doppler shift of the clutter is zero. Exploiting thischaracteristic, the zero-Doppler (DC) components of thereceived signal are removed before processing. B. Signal Separation Using Adaptive Beamforming
The minimum variance distortionless response (MVDR)[17]–[20], which is also known as the directionally con-strained maximization of power method [21], is an adap-tive beamforming technique that minimizes the power ofthe weighted sum of the signals y ( t ) = w H x ( t ) whilemaintaining a constant response in the desired direction,e.g., w H a = 1 and a is a mode vector for the de-sired direction. The MVDR method can suppress inter-fering signals if the source signals are uncorrelated, i.e.,E[ s i ( t ) s j ( t )] = 0 for i = j . If the estimate of a mode vectorˆ a j = [ˆ a ,j , ˆ a ,j , · · · , ˆ a M,j ] T for the j -th target is given, thesignal from the j -th target is estimated as ˆ s j ( t ) = w H j x ( t )using a MVDR weight w j = (cid:16) ˆ a H j R − ˆ a j (cid:17) − R − ˆ a j , (5)where R = E[ xx H ] is a correlation matrix. The unmixingmatrix W is obtained as W = [ w , w , · · · , w N ] H .As stated above, the MVDR method is simple to im-plement and is effective for suppressing interference if thesource signals are uncorrelated. Pulse waves measured atmultiple body parts are, however, highly correlated, whichdoes not satisfy the conditions for the MVDR method. Forcorrelated signals, the spatial-averaging technique [22] isoften used with the MVDR method, to un-correlate thesignals by reducing the size of the array. However, usingthis technique, the effective array size is diminished anddegrees of freedom lowered. In addition, the techniquecannot be applied to arbitrary arrays. For these reasons,we do not discuss the spatial-averaging technique in thefollowing sections. C. Signal Separation Using Independent Component Anal-ysis
Independent component analysis (ICA) is a technique ofblind signal separation, which decomposes a multivariatesignal into multiple non-Gaussian components (except forone component) that are statistically independent of eachother. In particular, JADE (Joint Approximation Diago-nalization of Eigenmatrices) [23] is a type of ICA; JADEuses the fourth-moment signals to decompose signals intoindependent components. Let us assume that the inputsignals have been whitened in preprocessing; their mean iszero (E { x i } = 0 (1 ≤ i ≤ M )) and they are uncorrelatedE { x i x j } = 0 (1 ≤ i = j ≤ M ). We also assume that the probability distribution function of each component s i (1 ≤ i ≤ N ) is symmetric; its odd moments are zero.Because the components are statistically independent, thefourth cross-cumulantcum( s i , s j , s k , s l ) = E { s i s j s k s l } − E { s i s j }{ s k s l }− E { s i s k }{ s j s l } − E { s i s l }{ s j s k } (6)is obtained ascum( s i , s j , s k , s l ) = (cid:26) c i ( i = j = k = l )0 (otherwise) (7)Using a matrix M = ( m ij ), the cross-cumulant isreduced to a matrix F i,j ( M ) = M X k,l =1 m k,l cum( x i , x j , x k , x l ) . (8)If the unmixing matrix is unitary W = ( u , · · · , u N ), weobtain F ( M ) = M X k,l =1 ( c i u T i M u i ) u i u T i (9)because cum( s i , s j , s k , s l ) = 0 unless the indices are allequal. By defining a diagonal matrix D as D ( M ) =diag( c i u T1 M u , · · · , c M u T M M u M , · · · ), we can express F ( M ) as F ( M ) = W D ( M ) W T . From this equation, it isfound that the diagonalization of F ( M ) gives the unmix-ing matrix W . By obtaining an eigenmatrix expansion ofthe tensor cum( x i , x j , x k , x l ), we obtain eigenmatrices thatare used to obtain multiple reduction matrices. We thenfind a unitary matrix W that simultaneously diagonalizesthe reduction matrices generated from the eigenmatricesusing the Jacobi method. IV. Signal Separation Using a MathematicalModel of Physiological Signals
In this section, we propose a mathematical model ofa body displacement due to a pulse wave. Using themathematical model, we formulate the signal separationprocess as an optimization problem. The objective func-tion is proposed to estimate an unmixing matrix W =[ w w · · · w N ] T . For simplicity, each weight vector is as-sumed to be normalized as | w i | = 1. A. Small-Displacement Approximation
The displacement caused by pulse waves is, typically, atmost 100 µ m, although its value depends on individual dif-ferences and body parts. We assume that the displacementis sufficiently small, compared with the wavelength λ = 3 . kd j ( t ) = 4 πd j ( t ) /λ ≪ π . Therefore, an echo can beapproximated as s j ( t ) = e j2 kd j ( t ) ≃ kd j ( t ) , (10)where a displacement can be estimated without the ∠ operation. In addition, Eq. (10) indicates that the I-Q plotof the echo on the complex plane can be approximated by a line segment, which can be expressed using the covariancematrix of the real and imaginary parts of the signal.For the i -th echo, the covariance matrix R i of thereal and imaginary parts of the complex-valued signal iswritten as R i = Z ∞−∞ (cid:20) ℜ [ˆ s i ( t )] ℜ [ˆ s i ( t )] ℑ [ˆ s i ( t )] ℑ [ˆ s i ( t )] ℜ [ˆ s i ( t )] ℑ [ˆ s i ( t )] (cid:21) d t, (11)where ˆ s i ( t ) = w T i x ( t ) and we should note that R i dependson w i . Let κ ( R i ) = | λ max /λ min | be the condition numberof the matrix R i , where λ max and λ min are the largestand smallest eigenvalues of R i , respectively. To obtain echoestimates that can be approximated by Eq. (10), we find w i ∈ C m that maximizes κ ( R i ) for i . Because R i is a 2 × κ ( R i ) = ( λ ( i ) + λ ( i )) / ( λ ( i ) − λ ( i )) , (12)where λ ( i ) = Z ∞−∞ (cid:0) ℜ [ˆ s i ] + ℑ [ˆ s i ] (cid:1) d t (13)and λ ( i ) = (cid:12)(cid:12)(cid:12)(cid:12)Z ℜ [ˆ s i ] −ℑ [ˆ s i ] d t (cid:12)(cid:12)(cid:12)(cid:12) +4 (cid:12)(cid:12)(cid:12)(cid:12)Z ℜ [ˆ s i ] ℑ [ˆ s i ]d t (cid:12)(cid:12)(cid:12)(cid:12) , (14)where the integration interval [ −∞ , ∞ ] is omitted. Weshould note that λ ( i ) and λ ( i ) are both positive and λ ( i ) > λ ( i ). Therefore, when the echo power λ isconstant, w i that maximizes κ ( R i ) also maximizes λ .The same discussion is valid for all i ( i = 1 , , · · · , n );the optimum unmixing matrix W should maximize anobjective function F ( W ) defined as F ( W ) = min ≤ i ≤ n λ ( i ) . (15)By increasing F ( W ), the estimated echo ˆ s i ( t ) on the com-plex plane becomes flat and elongated, which is consistentwith the small-displacement approximation in Eq. (10).However, this model alone cannot estimate multiple echoesfrom different parts of a human body. In the next section,we extend the objective function so that the problem canbe avoided. B. Simplified Pulse Wave Propagation Model
In this section, we consider the relationship between thedisplacements at multiple body parts. We introduce a sim-plified model of pulse wave propagation; the displacementscaused by the pulse wave are assumed to have similarwaveforms with a time delay corresponding to the pulsetransit time.First, we assume that each displacement d i ( t ) ( i =1 , , · · · , N ) is approximated by a constant multiple of atime-shifted template waveform d ( t ), i.e., d i ( t ) = α i d ( t − τ i ). We also assume that delay τ i > i , and the delays are sorted in an ascending order as τ < τ < · · · < τ N without losing the generality of themodel. With this assumption, body parts 1 , , · · · , N are arranged in order from the heart side to the terminal sideif the PWV is constant. The proposed model is formulatedas d ( t ) ≃ α d ( t − τ ) α d ( t − τ )... α N d ( t − τ N ) = d ( t ) ∗ α δ ( t − τ ) α δ ( t − τ )... α N δ ( t − τ N ) , (16)where δ () is the Dirac delta function and ∗ represents aconvolution integral. The transfer function between d i ( t )and d j ( t ) is obtained by deconvoluting the j -th componentwith the i -th component, and its impulse response g i,j ( τ )is calculated as follows: g i,j ( τ ) = Z ∞−∞ R ∞−∞ d j ( t ′ )e − j ωt ′ d t ′ R ∞−∞ d i ( t ′ )e − j ωt ′ d t ′ e j ωτ d ω. (17)If the above-mentioned conditions are satisfied, theimpulse response g i,j ( τ ) is approximated by g i,j ( τ ) ∝ δ ( τ + τ i − τ j ) . (18)Strictly speaking, the waveforms of d i ( t ) and d j ( t ) do notmatch perfectly, and, as a result, g i,j ( τ ) is not a deltafunction. To make g i,j ( τ ) as similar to a delta function aspossible, the fourth moment of g i,j ( τ ) is also included inour objective function as F ( W ) = Y ≤ i
0, which assumes a modelexpressed as g i,j ( τ ) ∝ δ ( τ + τ i − τ j ) and τ j > τ i . F ( W ) = Y ≤ i
Finally, the beam patterns (array factors) for weightvectors w , w , · · · preferably are orthogonal because eachweight vector is designed to receive a specific echo, whereasthe other weights are designed to reject the echo. This canbe formulated as F ( W ) = Y ≤ i A. Simulation Model and Parameters In this section, the performance of the conventionaland proposed methods are evaluated through numericalsimulations, which assumes a frequency of 79 GHz and anarray of M = 12 elements with element spacings of λ/ N = 2 is assumed; only two parts of thebody contribute to the scattering of echoes. Body parts 1and 2 are located at distances of 0.5 m and 0.3 m fromthe bottom of the antenna (see Fig. 1), i.e., x = − . x = 0 . − . ◦ and 14 . ◦ , respectively.Here, the power values of echoes from parts 1 and 2are assumed to be P = 0 dB and P = − P N = − 45 dB; thesignal-to-noise power ratio (S/N) for parts 1 and 2 are 45dB and 42 dB. The S/N values do not reflect the actualmeasurement of pulse waves because we are interestedonly in the time-varying components of the echoes thatcontain physiological information. The displacement dueto the pulse wave is d max = ± µ m, corresponding to aphase rotation of 19.0 ◦ . In this simulation, the power of thetime-varying component is lower than the time-invariantcomponent by 20 . 43 dB, resulting in effective S/N ratiosof 24.6 dB and 21.6 dB, respectively. We set γ = 10 andthe total measurement time T was set to T = 20 s.Fig. 2 shows the displacements d ( t ) and d ( t ) at bodyparts 1 and 2 assumed in the simulation, where thedisplacements are both triangular waves with a delay τ − τ = 300 ms. Since body part 1 is closer to the heartthan body part 2 and the distance between the parts is0.8 m, it is about 2.7 m/s when converted to PWV. Theleft and right figures in Fig. 3 show the I-Q plots of echoesfrom body parts 1 and 2 on the complex plane, both ofwhich draw an arc as approximated in Eq. (10). Fig. 1. System model with the radar and human in a prone position.Body parts 1 and 2 are targets in our simulation. Time (s) -100-50050100 D i s p l a c e m en t ( m ) Fig. 2. Simulated displacement waveforms d ( t ) (red) and d ( t )(black) of body parts 1 and 2. B. Performance Evaluation of Adaptive Array in Separat-ing Pulse Wave Echoes in a Simulation In this section, we apply the Capon method [17], [20]and MVDR [21] method to the simulated signal to evaluatetheir performance in separating radar echoes.Fig. 4 shows the I-Q plots of the echoes estimatedusing the MVDR method. The I-Q plots look differentfrom the ones in Fig. 3, which indicates that the signalseparation has not been performed correctly. Fig. 5 showsthe estimated displacement waveforms obtained from thesignals separated by the MVDR method. The waveformsalso seem different from the actual displacement in Fig. 2.Because the displacement waveforms at multiple bodyparts are similar and they are correlated, the MVDR -1 -0.5 0 0.5 1 Real -1-0.500.51 I m ag i na r y -1 -0.5 0 0.5 1 Real -1-0.500.51 I m ag i na r y Fig. 3. I-Q plots of simulated echoes s ( t ) (red) and s ( t ) (black). -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q Fig. 4. I-Q plots of complex-valued signals ˆ s ( t ) (left) and ˆ s ( t )(right) estimated using the MVDR method. Time (s) -1-0.500.51 D i s p l a c e m en t Fig. 5. Body displacement waveforms ˆ d ( t ) (red) and ˆ d ( t ) (black)estimated using the MVDR method. method does not work properly. The rms error ǫ i in thedisplacement waveform ˆ d i ( t ) is evaluated as ǫ i = s min η T Z T | d i ( t ) − η ˆ d i ( t ) | d t, (25)where the coefficient η is selected so that the error is min-imized because we are interested only in the displacementwaveform, not the scaling coefficient. For the estimationusing the MVDR method, the rms errors ǫ = 9 . µ m and ǫ = 10 . µ m were obtained. Fig. 6 shows impulse response g , ( τ ) calculated from the displacements ˆ d ( t ) and ˆ d ( t )obtained above. Although there are two peaks in g , ( τ ),the peak close to the actual PTT is located at 305.5 ms.Because the actual PTT assumed in the simulation is 300ms, the estimation error is 5 . C. Performance Evaluation of Independent ComponentAnalysis in Separating Pulse Wave Echoes in a Simulation In this section, we apply JADE-ICA to the simulatedsignals. To investigate its performance under an ideal con-dition, the actual number of targets is given. Fig. 7 showsthe I-Q plots of the estimated echoes. Compared with theestimates with the MVDR method (Fig. 4), the trajectoryseems closer to the actual one (Fig. 3). Fig. 8 showsthe displacement waveforms ˆ d ( t ) and ˆ d ( t ) estimated bythe JADE-ICA. Compared with the estimates using theMVDR method, the displacement waveforms in Fig. 8look closer to the actual ones. Despite this, both of the -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e Fig. 6. Impulse response g , ( τ ) estimated using the MVDR method. -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q Fig. 7. I-Q plots of complex-valued signals ˆ s ( t ) (left) and ˆ s ( t )(right) estimated using JADE-ICA. waveforms are distorted, with rms errors of ǫ = 5 . µ mand ǫ = 5 . µ m. Fig. 9 shows impulse response g , ( τ )calculated from the ˆ d ( t ) and ˆ d ( t ) obtained above. Thepeak of g , ( τ ) is located at 294.7 ms, and the estimationerror is 5 . D. Performance Evaluation of the Proposed Method inSeparating Pulse Wave Echoes in a Simulation Next, we investigate the performance of the proposedmethod in separating the simulated signals. The actualnumber of targets N = 2 is given in the same way as inthe evaluation of the MVDR method and JADE-ICA, andthus, the size of the unmixing matrix W is N × M = 2 × Time (s) -1-0.500.51 D i s p l a c e m en t Fig. 8. Skin displacement waveforms ˆ d ( t ) (red) and ˆ d ( t ) (black)estimated using the JADE-ICA algorithm. -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e Fig. 9. Impulse response g , ( τ ) estimated using the JADE-ICAalgorithm. The optimization problem W ∗ = arg max W ∈ C n × m F ( W ) (26)is solved to obtain an unmixing matrix W ∗ to estimate anecho ˆ s ( t ) = W ∗ x ( t ). The estimated echo ˆ s ( t ) is rotated sothat its principal component is directed in the real axis,and its real part is output as the displacement estimateˆ d ( t ).A genetic algorithm (GA) is used to maximize theobjective function in Eq. (26). Because the number ofcomplex unknowns is N M = 24 in W , a 24-dimensionalvector of complex numbers is used as a gene in the GA.The elements of unmixing matrix W are all initialized tobe 1’s. The number of individuals in each generation is100, and the value of the objective function is used as thefitness. Roulette selection with a probability proportionalto fitness is used for selection. In the crossover, the rowvectors in the selected two-generation individual matrixare exchanged to generate the two next-generation indi-viduals.Fig. 10 shows the I-Q plots of echoes ˆ s ( t ) (red line) andˆ s ( t ) (black line) estimated from the unmixing matrix W of the best individual in the 1st, 3rd, and 50th generationsof the GA maximizing F ( W ), where we see that the I-Qplots become flat by increasing F ( W ) as intended. Fig. 11shows body displacements ˆ d ( t ) (red line) and ˆ d ( t ) (blackline) estimated in the 1st, 3rd, and 50th generations of theGA. The estimated body displacement in the 50th gener-ation seems similar to the actual displacement shown inFig. 2, where the rms errors were ǫ = 2 . µ m and ǫ = 4 . µ m, indicating that the both errors are smaller than thosefor the MVDR method and JADE-ICA. Fig. 12 shows theimpulse responses g , ( τ ) obtained from these estimateddisplacements in the 1st, 3rd, and 50th generations of theGA-based optimization.Because ˆ d ( t ) and ˆ d ( t ) are the same in the first gener-ation, the impulse response g , ( τ ) has its peak at τ = 0,which reduces the values of F ( W ), lowering the valueof the objective function F ( W ). In the 50th generation,however, the impulse response g , ( τ ) has its value at τ > 0, satisfying the causality condition, increasing the -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q -1 -0.5 0 0.5 1 I -1-0.500.51 Q Fig. 10. I-Q plots of complex-valued signals ˆ s ( t ) (red, left) and ˆ s ( t )(black, right) estimated using the proposed method that maximizes F ( W ). Shown are the estimations from the 1st (top), 3rd (middle),and 50th (bottom) generations of the genetic algorithm. value of F ( W ). Specifically, the peak of g , ( τ ) in the50th generation is located at τ = 300 . . d and d using the proposed method, where 100 different randomseeds were used. The mean ± SD errors in estimating thedisplacements d and d were 2.9 µ m ± µ m and 4.5 µ m ± µ m, respectively, whereas using the MVDR andICA, the average errors were 9.6 µ m and 5.7 µ m. Fig. 14shows the rms errors in estimating the displacementsusing the ICA and the proposed method. For each S/N,100 random seeds were used, among which, errors largerthan 15 µ m were excluded in calculating the rms errors.The large errors were caused by insufficient numbers ofgenerations and population size in the GA. For S / N = 25and 30 dB, 2 cases out of 100 resulted in large errors;for S / N = 35 and 50 dB, 1 case out of 100 resulted inlarge errors. The figure shows that the proposed method Time (s) -1-0.500.51 D i s p l a c e m en t Time (s) -1-0.500.51 D i s p l a c e m en t Time (s) -1-0.500.51 D i s p l a c e m en t Fig. 11. Body displacements ˆ d ( t ) (red line) and ˆ d ( t ) (black line)estimated using the proposed method that maximizes F ( W ), for the1st (top), 3rd (middle), and 50th (bottom) generations of the geneticalgorithm. achieved a higher accuracy than did the ICA.Next, we investigate the performance of the 3 methodsin another simulation with x = − . x = 0 . 3, and P = P = 0 dB; the rest of the parameters were thesame. For the proposed method, simulations were run with100 random seeds. Among 100 random seeds, 3 seeds gavelarge PTT errors: 293.2, 293.0, and 294.0 ms, whereas theother random seeds gave small errors. Except for those3 cases, the mean error in estimating the PTT was 0.57ms ± µ m, 6.1 µ m, and 3.8 µ m ± µ m using the MVDR, ICA, and theproposed method, respectively. These results demonstratethe effectiveness of the proposed method. E. Application of the Proposed Method to ExperimentalData In this section, we apply the proposed method to exper-imental radar data from a participant. In the experiment,a 79-GHz millimeter-wave radar was used, as in the sim-ulation in the previous section. The radar is a frequency- -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e Fig. 12. Impulse response g , ( τ ) estimated using the proposedmethod that maximizes F ( W ). Shown are the estimations from the1st (top), 3rd (middle) and 50th (bottom) generations of the geneticalgorithm. RMS Error in Displacement ( m) C oun t Echo 1Echo 2 Fig. 13. Histogram of the rms errors in displacements using theproposed method. modulated CW system with a 99%-occupied bandwidth of3.5 GHz. The antenna is a MIMO array that uses threeelements for transmission and four elements for reception;a total of seven elements are on the same straight line,and the transmission and reception arrays have elementspacings of 2 λ and λ/ 2, respectively. Therefore, the vir- 25 30 35 40 45 50 55 60 S/N (dB) R M S E rr o r i n D i s p l a c e m en t ( m ) ProposedJADE-ICA Fig. 14. Average rms errors vs. S/N in estimating displacementsusing the ICA and the proposed method. tual array is a linear array with 12 elements with half-wavelength spacings, which is the same condition as thesimulation in the previous section.The origin of the xyz Cartesian coordinates is the centerof the array. The antenna array is in the xy plane, inwhich the array baseline is on the x -axis. The downwardvertical direction is the z axis. The array is installedfacing the floor at a height of 1.4 m from the top of theStyrofoam bed placed on the floor. The antenna elementshave a linear polarization in the y direction, and theelement beam widths in the E plane ( yz plane) and Hplane ( xz plane) are 8 ◦ and 70 ◦ , respectively. When thebeamforming method is used, the baseline length of thevirtual array is 11 λ/ 2, so the beam width of the arrayfactor in the H plane ( xz plane) is 9 . ◦ . The participantlies on the bed in a prone position with his body alignedalong the x direction. The distance between the upper partof the participant and the array center was about 1.25 m.The footprint of the element beam was 152 . × . x and y directions, which covered the majority ofthe person’s body, while suppressing echoes from the floor.Fig. 15 shows a photograph of the participant during themeasurement. Targets 1 and 2 correspond to the back andcalf, respectively. As a form of preprocessing, range gatingwas adopted; we removed signals from ranges other thanthat of the two body parts, to suppress clutter, as x ( t ) = Z r r x ′ ( r, t )d r, (27)where x ′ ( r, t )d r is a signal vector received for time t andrange r , and the two body parts are located within r ≤ r ≤ r , which was set manually.The proposed method was applied to the measuredsignal x ( t ), and used to investigate the method’s per-formance in separating radar echoes with experimentaldata. Fig. 16 shows the body displacements ˆ d ( t ) (redline) and ˆ d ( t ) (black line) estimated in the 1st, 3rd and50th generations of the GA optimization, where we seethat there is a delay between the displacement waveformsshown in red and black. The impulse responses g , ( τ )calculated from these waveforms are shown in Fig. 17. In Fig. 15. Photograph of a participant in the radar measurement. the 1st generation, a peak appears at τ = 0, and in the3rd generation, other peaks appear at τ = 0. Finally, inthe 50th generation, the main peak is shifted to τ > VI. Discussion A. Posture of Target Person In this study, we focused on the radar echoes from atarget person in a prone position only. The posture ofthe person determines the body parts contributing radarechoes, which affects the number of sources, the directionsof arrival, the necessary angular resolution and the skindisplacement waveforms that are used in our proposedmethod. It will be important to study how the postureof the target person affects the performance in separatingskin displacement signals and estimating PTT and PWV. B. Methods for Separating Correlated Signals Other than the proposed method, there are existingmethods for separating correlated signals. One exampleis the use of compressed sensing (CS). In the on-grid CSmethod, possible values for each parameter to estimate arerepresented by a set of discrete grid points, where most ofthe grid values are zeros, and this property is exploited toseparate signals. In recent works, in addition to the on-gird CS method, various off-grid and gridless CS methodshave been proposed and demonstrated to achieve bothhigh accuracy and robustness [24], [25]. The performanceof these methods, as applied to our problem setting, willbe investigated in our future work. VII. Conclusion In this study, we proposed a signal processing tech-nique for non-contact measurement of the pulse wavepropagating along the artery using an array radar system.The proposed technique allows us to measure the body Time (s) -1-0.500.51 D i s p l a c e m en t Time (s) -1-0.500.51 D i s p l a c e m en t Time (s) -1-0.500.51 D i s p l a c e m en t Fig. 16. Body displacements ˆ d ( t ) (red) and ˆ d ( t ) (black) estimatedusing the proposed method applied to the experimental data. Theestimations from the 1st (top), 3rd (middle) and 50th (bottom)generations of the genetic algorithm are shown. displacement at multiple parts of the human body simul-taneously. The echoes from multiple body parts are phase-modulated by a correlated displacement waveform, whichlowers the performance of conventional methods suchas the MVDR method and JADE-ICA. We proposed amathematical model of body displacements caused by thepulse waves. The model is based on a small-displacementapproximation and time-shift waveform approximation ofthe pulse waves. Based on this model, we formulated thesignal separation as an optimization problem. Then, aGA was adopted to maximize the objective function andobtain an unmixing matrix, to estimate the echoes frommultiple body parts. We also conducted an experimentwith a participant using a 79-GHz, 12-channel MIMOarray radar. Our study demonstrated that the accuraciesin estimating the body displacements and PTT were bothimproved by the proposed method, compared with theconventional MVDR and JADE-ICA methods, which indi-cates that the proposed method is a promising approach toachieve an accurate non-contact PWV measurement usingarray radar. -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e -1 -0.5 0 0.5 1 Delay (s) I m pu l s e R e s pon s e Fig. 17. Impulse response g , ( τ ) estimated using the proposedmethod applied to the experimental data. Shown are the estimationsfrom the 1st (top), 3rd (middle), and 50th (bottom) generations ofthe genetic algorithm. References [1] “Summary health statistics: National health in-terview survey, 2018, Table A-1b, A-1c” 2018.Accessed: Aug. 2, 2020. [Online]. Available:http://https://ftp.cdc.gov/pub/Health_Statistics/NCHS/NHIS/SHS/2018_SHS_Table_A-1.pdf .[2] D. Buxi, J.-M. Redouté, and M. R. Yuce, “Cuffless bloodpressure estimation from the carotid pulse arrival time usingcontinuous wave radar,” Proc. 2015 37th Annu. Int. Conf. IEEEEng. Med. Biol. Soc., Milan, Italy, pp. 5704–5707, doi:10.1109/EMBC.2015.7319687.[3] D. Buxi, J. Redouté, and M. R. Yuce, “Blood pressure estima-tion using pulse transit time from bioimpedance and continuouswave radar,” IEEE Trans. Biomed. Eng., vol. 64, no. 4, pp. 917–927, Apr. 2017, doi: 10.1109/TBME.2016.2582472.[4] D. Buxi, et al., “Systolic time interval estimation us-ing continuous wave radar with on-body antennas,” IEEEJ. Biomed. Health Inform., vol. 22, no. 1, pp. 129–139, Jan. 2018.[5] M. P. Ebrahim, et al., “Systolic blood pressure estimation usingwearable radar and photoplethysmogram signals,” Bari, Italy, pp. 3878–3882, doi:10.1109/SMC.2019.8914567.[6] M. P. Ebrahim, F. Heydari, T. Wu, K. Walker, K. Joe, J.-M. Redoute and M. R. Yuce, “Blood pressure estimation usingon-body continuous wave radar and photoplethysmogram invarious posture and exercise conditions,” Sci. Rep., vol. 9,no. 16346, 2019. doi: 10.1038/s41598-019-52710-8[7] M. Kuwahara, E. Yavari, and O. Boric-Lubecke, “Non-invasive,continuous, pulse pressure monitoring method,” Berlin, Germany,pp. 6574–6577, doi: 10.1109/EMBC.2019.8857439.[8] T. Lauteslager, M. Tømmer, T. S. Lande, and T. G. Constandi-nou, “Coherent UWB radar-on-chip for in-body measurementof cardiovascular dynamics,” IEEE Trans. Biomed. CircuitsSyst., vol. 13, no. 5, pp. 814–824, Oct. 2019, doi: 10.1109/TB-CAS.2019.2922775.[9] T. Tao, S. Hu, J. Peng, and S. Kuo, “An ultrawide-band radar based pulse sensor for arterial stiffnessmeasurement,” Proc. 2007 29th Annu. Int. Conf. IEEEEng. Med. Biol. Soc., Lyon, France, pp. 1679–1682, doi:10.1109/IEMBS.2007.4352631.[10] M. Tang, C. Liao, F. Wang, and T. Horng, “Noncon-tact pulse transit time measurement using a single-frequencycontinuous-wave radar,” Proc. 2018 IEEE/MTT-S Int. Mi-crow. Symp., Philadelphia, PA, USA, pp. 1409–1412, doi:10.1109/MWSYM.2018.8439326.[11] L. Lu, C. Li, and D. Y. C. Lie, “Experimental demonstrationof noncontact pulse wave velocity monitoring using multipleDoppler radar sensors,” Proc. 2010 Annu. Int. Conf. IEEEEng. Med. Biol., Buenos Aires, Argentina, pp. 5010–5013, doi:10.1109/IEMBS.2010.5627213.[12] R. Vasireddy, J. Goette, M. Jacomet, and A. Vogt, “Es-timation of arterial pulse wave velocity from Dopplerradar measurements: a feasibility study,” Proc. 2019 41stAnnu. Int. Conf. IEEE Eng. Med. Biol. Soc., Berlin, Germany,pp. 5460–5464, doi: 10.1109/EMBC.2019.8857644.[13] F. Michler et al., “Pulse wave velocity detection using a24 GHz six-port based Doppler radar,” Orlando, FL, USA, pp. 1–3, doi:10.1109/RWS.2019.8714521.[14] X. Wang, L. Wang, X. Li, and G. Bi, “Nuclear norm min-imization framework for DOA estimation in MIMO radar,” Signal Process., vol. 135, pp. 147–152, June 2017, doi:10.1016/j.sigpro.2016.12.031.[15] F. Wen, Z. Zhang, K. Wang, G. Sheng, and G. Zhang, “Angleestimation and mutual coupling self-calibration for ULA-basedbistatic MIMO radar,” Signal Process., vol. 144, pp. 61–67,Mar. 2018, doi: 10.1016/j.sigpro.2017.09.021.[16] L. Wan, X. Kong, and F. Xia, “Joint range-Doppler-angle esti-mation for intelligent tracking of moving aerial targets,” IEEEInternet Things J., vol. 5, no. 3, pp. 1625–1636, June 2018, doi:10.1109/JIOT.2017.2787785.[17] J. Capon, “High-resolution frequency-wavenumber spectrumanalysis,” Proc. IEEE , vol. 57, no. 8, pp. 1408–1418, Aug. 1969.[18] H. L. Van Trees, Detection, Estimation, and Modulation The-ory, Optimum Array Processing. Hoboken, NJ, USA: Wiley,2004.[19] J. Benesty, J. Chen, and Y. Huang, “A generalized MVDRspectrum,” IEEE Signal Process. Lett., vol. 12, no. 12, pp. 827–830, Dec. 2005.[20] J. Li, P. Stoica, and Z. Wang, “On robust Capon beamformingand diagonal loading,” IEEE Trans. Signal Process., vol. 51,no. 7, pp. 1702–1715, Jul. 2003.[21] K. Takao and N. Kikuma, “Tamed adaptive antenna array,” IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 388–394,Mar. 1986, doi: 10.1109/TAP.1986.1143821.[22] D. A. Linebarger and D. H. Johnson, “The effect of spatialaveraging on spatial correlation matrices in the presence ofcoherent signals,” IEEE Trans. Acoust. Speech Signal Process., vol. 38, no. 5, pp. 880–884, May 1990, doi: 10.1109/29.56037.[23] J. F. Cardoso and A. Souloumiac, “Blind beamforming for non-Gaussian signals,” IEE Proc.-F, vol. 140, no. 6, pp. 362–370,Dec. 1993, doi: 10.1049/ip-f-2.1993.0054.[24] J. Zheng, T. Yang, H. Liu, and T. Su, “Efficient data transmis-sion strategy for IIoTs with arbitrary geometrical array,” IEEETrans. Ind. Informat., doi: 10.1109/TII.2020.2993586.[25] J. Zheng, T. Yang, H. Liu, T. Su, and L. Wan, “Accurate de-tection and localization of UAV swarms-enabled MEC system,”