Deterministic method for non-contact vital signs monitoring via multiple radar signals and related monitoring system

Abstract

The present invention relates to the estimation of vital signs of living beings, human or animal, through radar signals of different types both in terms of the signals and the devices used for the transmission and reception of such signals. In particular, a new algorithm is developed for the real-time estimation of the heart rate optionally in combination with the respiratory rate and further capable of detecting apnea if desired. It employs a selection combining technique and a filtering method with memory which, in one embodiment, is inspired by the principles of Bayesian filtering applied to radar signals reflected from the body of the monitored living being. The invention also relates to the system and software implementing such method.

Description

"DETERMINISTIC METHOD FOR NON-CONTACT VI TAL S IGNS MONI TORING VIA MULTI PLE RADAR S IGNALS AND RELATED MONI TORING SYSTEM"DESCRIPTION

[0001] The present invention relates to the estimation of the vital signs of living beings, human or animal, through radar signals of various types, both in terms of the signals themselves and in terms of the devices used for the transmission and reception of such signals.

[0002] Among the innovative aspects, a new algorithm is developed for real-time estimation of the heart rate, optionally in combination with the respiratory rate, and further capable of detecting apnea if desired. It employs a selection combining technique and a memory-based filtering method, which in one variant is inspired by the principles of Bay esian filtering. Experimental results show that the proposed algorithm, tested on commercial radar devices, outperforms the well-known vital sign estimation methods available in the technical literature, at the cost of a limited increase in computational complexity.INTRODUCTION

[0003] In the past two decades, significant attention has been given to the development of accurate methods for contactless vital signs monitoring (VSM) of human beings and, in particular, for monitoring heart rate (HR) and breathing rate (BR) [1-3], This interest is motivated by the fact that, in various circumstances (for example, in the case of infected patients or individuals with mental health problems or suffering from extensive bums or wounds), the use of wearable sensors or traditional monitoring methods such as electrocardiography (ECG) or photoplethysmography (PPG) is neither feasible nor advisable [4], Furthermore, contactless methods for VSM can be very useful for monitoring the sleep of elderly people and children and to enable prolonged monitoring of patients at home.

[0004] A relevant technical option for contactless VSM is represented by modem integrated radar systems operating at millimeter waves. In fact, these devices, when positioned in front of a human chest wall, are potentially able to detect small physiological movements caused by both heartbeat and respiration. Compared to traditional methods such as ECG and PPG, radar systems make it possible to continuously and promptly monitor BR and HR without requiring additional workload for nurses. Furthermore, unlike light-based sensors (such as LIDAR and cameras), they preserve the privacy of monitored individuals and their measurements are not affected by clothing.

[0005] The same advantages described above also apply in the case of animal monitoring; in this case too, the contactless sensor allows for continuous measurement of vital parameters in order to assess the animal's well- beintystress level.

[0006] The technical literature on radar-based VSM has considered the use of different radar technologies, including continuous wave Doppler (CW Doppler) radars [5-9], impulse radio ultra-wideband (IR-UWB) radars

[0010] , and frequency-modulated continuous wave (FMCW) radars [11-14], Moreover, recently, several research efforts have been dedicated to evaluating the benefits of exploiting colocated multiple-input multiple-output (MIMO) radars, i.e., radar devices equipped with closely spaced arrays of transmitting (TX) and receiving (RX) antennas [10, 11, 15],

[0007] However, regardless of the technology adopted and the type of radar signals involved, a major challenge in the field of radar-based VSM research remains the development of reliable and accurate signal processingmethods for extracting information about vital signs from radar echoes. Most of the methods currently available in the technical literature are deterministic, as their derivation is based on specific mathematical models of the electromagnetic signal reflected by the human (or animal) chest; moreover, they extract vital signs through spectral analysis of the phase of the received signal.

[0008] Notable examples of this approach are represented by the arctangent demodulator (AD) [5] and the cumulative phase gradient (CPG) technique

[0016] ,

[0009] In particular, the AD method was originally proposed for detecting cardiopulmonary activity through a quadrature CW Doppler radar [5], The availability of a quadrature receiver allows for solving the so-called null point detection problem, i.e., the dependence of detection sensitivity on the target’ s position. However, the accuracy of the AD is limited by the presence of a direct current (DC) offset

[0017] and multiple high-order harmonics and intermodulation products.

[0010] The DC offset problem is due not only to the electronic components but also to reflections from stationary obj ects and other parts of the human (or animal) body different from the monitored subj ecf s chest wall.

[0011] The CPG method, on the other hand, was proposed for the simultaneous estimation of the BR of multiple human subjects

[0016] , In this case, relevant information on the phase of the received signal is acquired through a phase gradient (PG) operation, followed by spectral analysis. This method eliminates the DC offset and is capable of suppressing the contribution of respiratory higher-order harmonics (HOH) and intermodulation products caused by the presence of two subjects with (potentially) different BRs. However, its performance is influenced by the presence of noise (e.g., random body movements, RBMs) in the analyzed spectrum; furthermore, its use in HR estimation has not yet been studied.

[0012] It should also be noted that: a) as recently demonstrated in

[0016] , the AD is outperformed by the CPG method in the estimation of BR; b) as far as known, the use of these methods for apnea detection has not yet been studied.

[0013] This patent application proposes a new robust and accurate radar-based method for VSM of a single subject.

[0014] This method exploits the PG operation proposed in

[0016] for extracting relevant information from radar echoes. However, unlike

[0016] , the proposed method employs: a) a method for detecting and estimating dominant spectral components. In one implementation, the single frequency estimation and cancellation (SFEC) method recently proposed in

[0018] is used, although other state-of- the-art methods can also be employed in experimental trials, such as:1) the CLEAN algorithm (https: / / en. wikipedia. ortywiki / CLEAN falgorithm)), using samples from the discrete Fourier transform (DFT) of the input signal and operating via serial cancellation;2) the RELAX algorithm, an extension of the previous one including a refinement step; or3) the MUSIC algorithm (Multiple Signal Classification, hPps: / tyitykipedia.Qj>' wiki / MUSIC_(algorithm)), typically used for angle-of-arrival estimation in systems with multiple antennas and based on factorization of the observed noise covariance matrix affecting the received signal;b) an innovative technique, inspired by Bayesian filtering, to emphasize relevant spectral components and suppress spurious ones; c) a known number of generic radar signal emitting sources, arranged in a known configuration and configured to emit electromagnetic detection signals to be received, when reflected, by a plurality of receiving elements (antennas) connected to a processing unit.

[0015] A peculiar aspect is that such elements can be coordinated in time and do not require mutual phase synchronization, the latter being a condition necessary for known systems using more or less advanced beamforming techniques for detection via electromagnetic waves. In the present invention, a temporal base is used among the transceiver elements solely to allow coordinated transmission and / or reception among them, and the received signals are used to improve the quality of the acquired measurements through a specifically designed selection combining mechanism.

[0016] The transmitting and receiving elements can be organized in multiple single-input single-output (SISO) radars, one or more multiple-input single-output (MISO) radars, one or more single-input multiple-output (SIMO) radars, or one or more multiple-input multiple-output (MEMO) radars of commercially available types. It is essential to emphasize that the chosen configuration may also include different types of radar among those listed above (e.g., 2 SISO radars and 1 MEMO radar), but must necessarily provide a total number >1 of virtual antennas (VAs, see the section DETAILED DESCRIPTION for the definition of virtual antenna).

[0017] Notation. Throughout the present document, (•)* denotes the complex conjugate, Ja(x) denotes the Bessel function of the first kind of order a evaluated at x, 5R{x] e 3{x} denote the real and imaginary parts, respectively, of the complex variable x.DETAILED DESCRIPTION

[0018] The invention is detailed below also with reference to the attached figures, in which:Figure 1 : Overall architecture of the radar-based VSM system developed in the present work;Figure 2: Block diagram of the SIMO FMCW radar, transmission, reception and pre-processing parts;Figure 2b: Block diagram of the SIMO SFCW radar, transmission, reception and pre-processing parts;Figure 2c: Block diagram of the SIMO CW Doppler radar, transmission, reception and pre-processing parts;Figure 2d: Block diagram of the SIMO IR-UWB radar, transmission, reception and pre-processing parts;Figure 3 : Representation of the time evolution of the instantaneous frequency characteristic of each of the employed FMCW radars with a single TX antenna;Figure 3b: Representation of the time evolution of the instantaneous frequency characteristic of each of the employed SFCW radars with a single TX antenna;Figure 3c: Representation of the time evolution of the instantaneous frequency characteristic of each of the employed CW Doppler radars with a single TX antenna;Figure 3d: Representation of the time evolution of the instantaneous frequency characteristic of each of the employed IR-UWB radars with a single TX antenna;Figure 3e: Representation of the time evolution of the instantaneous frequency characteristic of each of theemployed FMCW radars: radar board 1 has a single TX antenna, while radar board 2 has 2 TX antennas;Figure 4: Representation of the time evolution of the instantaneous frequency of the signals irradiated during a superframe; two FMCW radars are considered;Figure 4b: Representation of the time evolution of the instantaneous frequency of the signals irradiated during a superframe; two SFCW radars are considered;Figure 4c: Representation of the time evolution of the instantaneous frequency of the signals irradiated during a superframe; two CW radars are considered;Figure 4d: Representation of the time evolution of the instantaneous frequency of the signals irradiated during a superframe; two IR-UWB radars are considered;Figure 4e: Representation of the time evolution of the instantaneous frequency of the signals irradiated during a superframe: a system consisting of one FMCW radar with 2 IX antennas, one IR-UWB radar with 1 TX antenna, and one SFCW radar with 1 TX antenna is considered;Figure 5 : Representation of a range profile for NF= 317 from a single virtual antenna;Figure 6: Representation of the generation of data blocks;Figure 7a: Block diagram of the algorithm (processing steps);Figure 7b: Block diagram of the algorithm (processing and pre-processing steps) for a plurality of virtual antennas and radar boards;Figure 7c: Detailed block diagram of the algorithm (processing and pre-processing steps);Figures 8a and 8b: Representation of the amplitude spectrum of the PG computed on a data block for one of the virtual antennas and corresponding generated GM;Figure 9: Representation of the adopted measurement configuration;Figure 10: Representation of the amplitude spectrum observed on a single virtual antenna of each radar device used by the method;Figure 11 : Representation of the amplitude spectra observed on the two virtual antennas of a single radar;Figure 12: Representation of the effects of apnea on the amplitude spectrum observed on a given virtual antenna;Figure 13: Representation of the amplitudes, versus <5b / M, of the multiple sinusoidal components of the sequences (a) {s[ ]}; (b) { / [ ]}; (c) {^[ ]}.

[0019] With reference to the figures, the overall architecture of the radar-based VSM system developed in the present work is illustrated in Figure 1. In this specific representation, three radars are mounted on a stand in front of the monitored person and are connected to the same personal computer; this latter device provides the timing reference that coordinates the activity of the radars and processes the acquired measurements to extract information about the vital signs.

[0020] In the example described herein, the system consists of a known number NDof radar boards (in the figure, / Vo=3) positioned in proximity of an observation area intended to accommodate a subject to be monitored, each equipped with / VTXTX antennas and NRXRX antennas. These devices, assumed for simplicity to be of the same type, are installed on a support and oriented so that the signals reach the subj ect to be monitored.

[0021] When a radar board is equipped with at least one antenna array (TX or RX), i.e., it can be classified as MISO, SIMO, or MIMO, each TX / RX antenna pair can be replaced with an equivalent virtual antenna (V A). Thi s VA is virtually positioned at abscissa xVA=(XTX + XRX) / 2 and ordinate yVA= (yTX+ yRX) / 2, where xTX(xRX) and yTX(yRX) are the abscissa and ordinate coordinates of the selected TX (RX) antenna (see hityvken wiki pedia:org / wiki / MIM()^ j^dar). Note that a necessary condition fbr a TX / RX antenna pair to be replaced by the virtual element with the position described above is the phase synchronization of the two said antennas. In principle, a TX / RX antenna pair with such elements physically located on two separate radar boards can be replaced by the corresponding virtual antenna if phase synchronization between them exists.

[0022] This concept of antenna "virtualization" also applies in the case of single antennas, i.e., in the case of two or more individually operating antennas, including at least one TX and one RX (SISO antennas), synchronized in phase and coordinated in time.

[0023] For each radar board considered, the number / VVAof available VAs is particularly important; they are defined as the / VVA= NTXNRXelements resulting from the combination of physical antennas (TX / RX antenna pairs).

[0024] The different radars are also temporally coordinated by a personal computer (PC), which operates them one at a time in a round-robin manner (i.e., they are turned on and off in sequence, so that, at any given moment, only one radar is active). As will be seen below, this coordination is aimed at and limited to avoiding the overlap of transmitted signals, whereas operations that require synchronization of received signals (such as beamforming in transmission or reception) are not strictly contemplated and thus not essential to the operation of the invention. In light of what was previously described, this particular technical solution does not allow for obtaining virtual antennas corresponding to TX / RX antenna pairs physically located on different radar boards.

[0025] For the technical solution presented in this invention, the total number of available VAs JVT0T= S i Jo NyAis of particular importance, where NyArepresents the number of VAs provided by the i -th radar board.

[0026] Furthermore, the aforementioned round-robin operation mode is to be considered as one of the possible implementations of a more general concept of coordinated transmission / reception operations among the various elements of the system; for example, the invention is fully compatible with management according to known pseudorandom activation sequences of the devices involved in the transmission / reception of radar detection signals.

[0027] The flow of measurements made available by each radar is sent to the aforementioned PC, on which software that implements one or more operational steps of the described method is executed. The remaining part of this section is dedicated to the description of the signal models relevant to the various types of radar signals of interest for achieving the objectives of the invention.

[0028] Generally, as better described below, the method includes a transmission step of one or more radar signals, a subsequent reception step of said signals, a pre-processing phase of the received signals, and a processing step of said pre-processed signals. In order to optimize the procedures and provide maximum flexibility in the application of the method, the pre-processing step has been divided in parts that are specific to the type of radar signal used (transmission, reception, and pre-processing); the processing step, instead, is applicable regardless of the type ofradar signal used are defined.

[0029] Consequently, the following sections present the different approaches according to the radar signal type, distinguishing between FMCW, SFCW, CW Doppler, and IR-UWB signals.Pre-processing of radar signals obtained by FMCW

[0030] As shown in Figure 2, the radio frequency (RF) signal radiated by the i -th radar (with i = 0, 1, ... , ND— 1) is generated by a voltage-controlled oscillator (VCO), characterized by a free-running frequency fc(assumed, for simplicity, to be independent of i) and powered by a periodic ramp generator. Therefore, the transmitted signal consists of a sequence of chirps, each of which has a durationT0= T + TR(1) where T and TRrepresent the ramp period and the reset time, respectively. In the considered radar system, the transmission of chirps is organized in frames; more precisely, each frame is composed of Ncconsecutive chirps and its duration isTv= NcT0+ 7 (2) where Tbcalled idle time, is the duration of the time interval separating each pair of consecutive frames.

[0031] The time evolution of the instantaneous frequency of the radiated signal is shown in Figure 3, where Nc= 3 is assumed. It is important to underline that, in the considered radar system, 7) is much greater than To. This allows allocating transmission and reception of each radar device in the idle time of the remaining ND— 1 devices, so as to avoid mutual interference.

[0032] All chirps transmitted (and received) by the NDdevices during the time interval between the beginning of the transmission of the first chirp within a frame related to the d-th radar device (d G {0,1, ... , ND— 1}) and the beginning of the first chirp in the next frame of the same radar device belong to the same superframe; a simple example of the time evolution of the instantaneous frequency of the signals radiated during a superframe is shown in Figure 4, referring to the case where ND= 2. Note that, in general, the duration of each superframe, denoted as TSF, equals the duration of the frame TF.

[0033] Hereinafter, the term frame will be used when referring to a single radar device. Conversely, when referring to multiple radar devices, the term superframe will be used.

[0034] The electromagnetic signal radiated by each radar is reflected by the chest of the monitored person. The derivation of a mathematical model for the radar echo requires the knowledge of the temporal evolution of the chest displacement, i.e., the time-varying distance between the radar and the chest itself. The chest displacement consists of [11,16]: a) relatively large movements due to the inspiration and expiration phases; b) small oscillations due to cardiac activity; c) RBMs of unpredictable amplitude.

[0035] If RBMs are ignored, i.e., it is assumed that the monitored subject does not perform RBMs or performs them only for time intervals limited with respect to the observation period, and a point-target model is adopted forthe human chest, the chest displacement at time t can be expressed as [6] :JR(t) A R(t) — Ro= <5b(t) + <5h(0 (3) where 6b(t) (<5h(0) indicate the contributions ofbreathing (heart) toand R(t) (Ro) represent the radar- to-chest distance at time t (distance measured in the absence of activity; note that, in general, this distance varies from radar to radar in the system shown in Figure 1). Accurate models for <5b(t) and <5h( t) were proposed in

[0019] , However, in the present work, the simplified point-target model (SPTM) illustrated in [20-22] is adopted. For this reason, it is assumed thatfor both the inspiration and expiration phases of the respiration period, and<5h(0 = bh,Mcos(mht), (5) for the heart displacement, wherehere, / b(Tb) and / h(Th) are respectively the breathing rate (breathing period) and the heart rate (heart period), and <5b / Mand <5h / Mrepresent the absolute maximum values of the breathing displacement and heart displacement, respectively.

[0036] It is worth recalling that: a) typical values of <5b / Mfor a human subject range from 8 to 24 mm, while those of <5h / Mrange from 0.4 to 0.5 mm

[0010] ; b) the human breathing (cardiac) frequency for an adult is usually between / b min= 0.1 Hz ( / h.min=0.883 Hz) and / b max= 0.5 Hz ( / tymax= 2.667 Hz) [11,16,23], Consequently, the periods are Tb / min= 10 seconds (Th min= 1.2 seconds) and Tb max= 2 seconds (Th max= 0.375 seconds).These values may reasonably differ if the monitored subject is an animal or even if the subject is a child whose standard parameters differ (typically greater than those of an adult). Therefore, the use of the present invention for monitoring vital signs of animals differs only in the setting of characteristic parameters, well known in the state of the art.

[0037] The embodiment of the following description refers to a single SIMO radar device, i.e., the active one.

[0038] As shown in Figure 1, where monitoring of a seated human subject is exemplified, the electromagnetic echo coming from the chest of the monitored person is captured by the NRXantennas forming the RX array. Furthermore, the NRXreceived signals are subjected to amplification, quadrature downconversion, and analog-to- digital conversion; these three tasks are performed by a low noise amplifier (LNA) and, for both in-phase and quadrature components, by a mixer and band-pass filter and by an analog-to-digital converter (ADC), respectively.

[0039] The resulting sample stream feeds a pre-processing stage, to which the remainder of this section is devoted. Note that, in the description of the mathematical models of the signals appearing in this section, their dependence on the radar index (i.e., on i) is neglected to facilitate notation, and the following assumptions are made on the chestdisplacement model (3): a) the displacement J R [n] is always much smaller than Ro(in practice, the chest displacement is on the order of a few millimeters, while the chest-to-radar distance is on the order of one meter); this derives from the working assumptions that the target model is point-like or simplified, excluding mathematical models working on distributed-type models. b) in the absence of RBMs, i.e., assuming that the monitored subject does not perform random movements or performs them only for time intervals limited compared to the observation period, the displacement undergoes small variations within each frame, since the frame period (or duration) TF(not exceeding 0.1 s in VSM) is significantly smaller than the heart period THR(whose minimum is about 0.4 s).

[0040] Consequently, the variations observed in the frequency of the downconverted received signals within each frame and, for a given frame, from antenna to antenna, as well as those observed in the phase of the same signals within each chirp can be neglected. Based on all this, the fc-th complex sample acquired through the v-th RX antenna during the s-th chirp interval of the n-th frame can be approximated as

[0024] :with v = 0,1, ... , lVRx — 1, k = 0,1, ... , 1V — 1, and s = 0,1, ... , NC— 1;

[0041] here, Tsis the sampling interval, N is the total number of samples acquired in a single chirp interval, [k, s, n] is the complex Gaussian noise sample affecting x® [k, s, n], and a® [s], fn, and t / Av)[s, n] indicate the amplitude, frequency, and phase, respectively, of the useful component of x^ [k, s, n] .

[0042] It can be shown that [11, 12, 25]: fn =j 2T (7) andwhere g = B / T is the chirp slope, B is the amplitude of the swept frequency interval (i.e., the radar bandwidth), A = is the radar wavelength, and R [n] is the distance of the v-th RX antenna from the chest during the n-th frame interval.Based on model (3), (8) can be rewritten as ip^[s, n] = ty0+ Jty00 [s, n], (9) where tyo = Y A «o (10) is a constant phase shift depending on the distance Ro, andAip(v)[s, n} = AR(v>[s, n} (11) represents the phase displacement due to the range variation R [s, n] .

[0043] The same considerations previously reported for a single SIMO radar can also be extended to the case of a single MISO radar.

[0044] In this case, the single RX antenna captures the electromagnetic echo coming from the chest of the monitored subject and generated by the electromagnetic radiation transmitted by the / VTXantennas forming the TX array. As shown in Figure 3e, where the time evolution of the instantaneous frequency characteristic of a system composed of a radar device with a single TX antenna and a radar device with two TX antennas is shown, the different TX antennas of the same radar device transmit chirps in time division multiplexing (TDM); this means that each frame consists of the echo received by the RX antenna and generated by the chirps transmitted by each single TX antenna. Practically, the / VTXtransmitted signals are used in the quadrature downconversion together with the received signal (after amplification). Subsequently, analog-to-digital conversion takes place.

[0045] The resulting sample stream feeds a pre-processing stage entirely analogous to that described previously for the SIMO radar. The phase displacement signal due to the range variation J R(v>[s, n] is expressed by (11); however, in this case, v = 0,1, ... , ATX— 1.

[0046] The same considerations previously reported for single SIMO and single MISO radars can also be extended to the case of a single MIMO radar.

[0047] In this case, the ARXantennas forming the RX array capture the electromagnetic echo coming from the chest of the monitored subj ect and generated by the signals transmitted by the / VTXantennas forming the TX array. The different TX antennas transmit chirps in TDM; this means that each frame consists of the echo received by each RX antenna and generated by the chirps transmitted by each single TX antenna. Practically, the / VTXtransmitted signals are used in the quadrature downconversion together with the ARXreceived signals (after amplification). Subsequently, analog-to-digital conversion takes place.

[0048] The resulting sample stream feeds a pre-processing stage entirely analogous to that described previously for the SIMO radar. The phase displacement signal due to the range variation J R [s, n] is expressed by (11); however, in this case, v = 0,1, ... , iVTXiVRX— 1.

[0049] The same considerations previously reported for single SIMO, single MISO, and single MIMO radars can also be extended to the case of a single SISO radar.

[0050] In this case, the RX antenna captures the electromagnetic echo coming from the chest of the monitored person and generated by the signal transmitted by the TX antenna. The received signals are subjected to amplification, quadrature downconversion, and analog-to-digital conversion.

[0051] The resulting sample stream feeds a pre-processing stage entirely analogous to that described previously for the SIMO radar. The phase displacement signal due to the range variation J R [s, n] is expressed by (l lbis)which is entirely equivalent to (11), but loses dependence on the VA index v .

[0052] It is important to underline that, although there exist radars equipped with real demodulators, i.e., devices that provide only the in-phase component of the demodulated signal, for the method described hereafter, the use of the I / Q demodulator is mandatory since both the in-phase and quadrature signals are used.

[0053] The estimation of vital signs is based on measurements acquired over NFframes; hereinafter, attention is focused, for the considered radar device, on the frames associated with n = 0,1, ... , NF— 1, and the processing performed on the set [x(v>[k, s, n] }, which collects the NCN complex measurements acquired on such frames, is described. This processing evolves through the three phases described below for the v-th VA, where the first phase is preferable but not strictly necessary to achieve the aims of the invention.1) Combination of chirps received during the same frame - The measurements acquired in the Ncchirp intervals of the n-th frame are combined to improve the signal-to-noise ratio (SNR). This step is based on the hypothesis that the chirp duration is much smaller than the minimum heart period (i.e., that To« fh min, with fh min= l / / h,max = 0.37 s), so that the approximationholds for s G {0,1, ... , 1VC— 1}.The combination is obtained throughfor k = 0,1, ... , N — 1 and n = 0,1, ... , NF— 1; this produces the new set {x(v>[k, n]} consisting of NNFcomplex measurements. A simple model for x [k, 7i] can be derived by substituting the right-hand side (RHS) of (6) into that of (13) and taking into account approximation (12); thus,where and

[0054] 2) Target detection -For each frame (i.e., for each n), the sequence [x(v>[k, n] } composed of N elements is subjected to zero-padding followed by a discrete Fourier transform (DFT) of order ZV0; here, No= M • N, where M indicates the selected oversampling factor. This yields the complex matrix of dimension NoX NFXMA [%(v)[Z, n]], (17) whereThen, the frequency bin index (i.e., range bin) associated with the target (i.e., the chest of the monitored subject) is estimated as f<v) [n] = argwhere the limits Zminand Zmaxon the search interval are chosen based on prior knowledge of the subj ecf s distance from the radar. Note that the use of a DFT of order Noin (18) divides the frequency interval of interest into Nodistinct frequency (range) bins.3) Extraction of relevant spectral information - Given Z^ [n] from (19), the relevant spectral information is extracted from matrix(17) by calculating the quantityforn = 0,1,

[0055] The sequences= 0,1, ... , NF— 1} feed the algorithm. The processing performed by this algorithm is based on the hypothesis that, for any v, the measurement s(v>[n] can be modeled as (the derivation of the following result is reported in Appendix A, Details on range compression and target selection)for n = 0,1, ... , Np — 1; here, b(v>represents the amplitude of the chest echo (depending on its position and reflectivity), [n] is the noise affecting the considered measurement and (see (9)-( 11))Pre-processing of radar signals obtained by SFCW

[0056] As shown in Figure 2b, the transmitter of an SFCW radar is similar to that of an FMCW radar, with the only difference being that the ramp generator of the latter is replaced by a waveform step generator. Therefore, the instantaneous frequency of the signal generated by the VCO used in an SFCW changes in steps within each frequency sweep radiated. The time evolution of the instantaneous frequency of the radiated signal is shown in Figure 3b, where / Vc= 3 is assumed. In this figure, TR, T, Toand 7) represent respectively the reset time, the frequency sweep duration, the total frequency sweep duration (To= TR+ T) and the idle time; instead, N and d / represent the total number and the width of each frequency step, respectively. Note that if Ncindicates the total number of frequency sweeps forming a single frame, each frame lasts Tp = NCTOseconds. It is important to underline that, in the considered radar system, 7) is much greater than To. This allows allocating transmission and reception of each radar device in the idle time of the remaining ND— 1 devices, so as to avoid mutual interference.

[0057] All frequency sweeps transmitted (and received) by the NDdevices during the time interval between the beginning of the transmission of the first frequency sweep within a frame related to the d -th radar device (d G {0, 1, ... , ND— 1}) and the beginning of the first frequency sweep in the next frame of the same radar device aresaid to belong to the same superframe; a simple example of the time evolution of the instantaneous frequency of the signals radiated during a superframe is shown in Figure 4b, referring to the case ZVD= 2. Note that, as in the FMCW radar case, the duration of each superframe, denoted as TSF, equals the duration of the frame TF.

[0058] If it is assumed that the sampling period Tsis equal to the duration of each frequency step (i.e., the sampling frequency fs= 1 / Tsis equal to J / ), a single complex sample is acquired at the RX side during each frequency step. Furthermore, in this case, the Zc -th complex sample available for the n-th frequency sweep can be expressed aswhere it can be shown that Tn= 2R0 / c. Note that znis a delay proportional to the distance between radar and monitored subject assuming no movement.

[0059] Observing the previous equation and (18), the only difference is represented by the sign of the argument of the complex exponential appearing on the RHS; for this reason, an SFCW radar can be seen as the dual of an FMCW radar. In light of this, the pre-processing steps described previously for the FMCW are indeed valid also for the SFCW with the only exception of equation (18), which is replaced bywith Z = 0,1, ... , N0— 1 and n = 0,1, ... , NF— 1.

[0060] The sequences= 0,1, ... , ZVF— 1} feed the algorithm. The processing performed by this algorithm is based on the hypothesis that, for any v, the measurement s(v>[n] can be modeled as (the derivation of the following result is reported in Appendix A, Details on range compression and target selection)for n = 0,1, ... , NF— 1; here, b(v>represents the amplitude of the chest echo (depending on its position and reflectivity), [n] is the noise affecting the considered measurement and (see (9)-( 11))Pre-processing of radar signals obtained by CW Doppler

[0061] As shown in Figure 2c, the transmitter of a CW Doppler radar comprises a local oscillator (LO) generating atone at frequency f0. The sinusoidal signal generated by the LO undergoes upconversion and amplification before transmission. The time evolution of the instantaneous frequency of the radiated signal is shown in Figure 3c. In this figure, TR, T, To, and 7) represent respectively the reset time, the duration of the continuous tone, the total duration of the continuous tone (To= TR+ T), and the idle time.

[0062] In this scenario, a frame consists of Nccontinuous tones. It is important to underline that, in the considered radar system, 7) is much greater than To. This allows allocating transmission and reception of each radar device in the idle time of the others (1VD— 1), thereby avoiding mutual interference. All continuous tones transmitted (and received) by the NDdevices during the time interval between the beginning of transmission of the first continuous tone within a frame related to the d-th radar device (d G {0,1, ... , ND— 1}) and the beginning of the first continuous tone in the next frame of the same radar device belong to the same superframe; a simple example of the time evolution of the instantaneous frequency of the signals radiated during a superframe is shown in Figure 4c, referring to the case ND= 2. Note that, as in the FMCW radar case, the duration of each superframe, denoted TSF, equals the duration of the frame TF.

[0063] At the RX side, the electromagnetic echo received and generated following the tones transmitted by the TX antennas is amplified by a low noise amplifier (LNA) and converted to baseband to extract the in-phase and quadrature components, expressed for the n-th ADC sample of the s-th tone of the frame bywhere, for the v-th virtual antenna (VA),represents the amplitude of the useful component of the signal, [s, 7i ] is the complex Gaussian noise sample affecting s'^ [s, TI], and (see (9)-(l 1))Following the steps previously described for the FMCW radar, the tones received during the same frame are summed asThe sequences {s^ [n]; n = 0,1, ... , NF— 1} feed the algorithm.Pre-processing of radar signals obtained by IR-UWB

[0064] The architecture of an IR-UWB radar is depicted in Figure 2d. At the TX side, a Gaussian pulse generator, triggered by a square wave generator, is used to generate the baseband signalwhere p (t — TIT0) represents the TI -th transmitted pulse and Tois the pulse repetition frequency.

[0065] The time evolution of s(t) is shown in Figure 3d. In this figure, TR, T, To, and 7) represent respectively theresettime, thepulse duration, thetotal pulseduration(T0= TR+ T), and the idle time. In this scenario, a frame consists of Ncpulses. It is important to underline that, in the considered radar system, 7} is much greater than To.This allows allocating transmission and reception of each radar device in the idle time of the others (1VD— 1), thereby avoiding mutual interference.

[0066] All pulses transmitted (and received) by the NDdevices during the time interval between the beginning of the transmission of the first pulse within a frame related to the d -th radar device (d G {0, 1, ... , ND— 1}) and the beginning of the first pulse in the next frame of the same radar device belong to the same superframe; a simple example of the time evolution of the amplitude of the signals radiated during a superframe is shown in Figure 4d, referring to the case / VD= 2.

[0067] Note that, as in the FMCW radar case, the duration of each superframe, denoted TSF, equals the duration of the frame TF. The signal s (t) is subjected to upconversion and amplification before transmission.

[0068] It is worth noting that the system previously described, comprising NDradar boards of the same type (i.e., for example, all FMCW type), is immediately extendable to an equivalent version consisting of NDradar boards of different types. An example of the time evolution of the instantaneous frequency characteristic in the case of a system composed of an FMCW radar (2 TX antennas), an IR.-UWB radar (1 TX antenna), and an SFCW radar (1 TX antenna) is shown in Figure 4e.

[0069] At the RX side, since the pulses of s(t) are not temporally overlapped, attention can be focused on the echo generated by the target (patient) in response to the s-th pulse. The received signal is thus preamplified and converted to baseband by a quadrature demodulator before undergoing analog-to-digital conversion at frequency fs= 1 / Ts, where Tsis the sampling period. The complex signal obtained for the fc-th sample during the s-th pulse of the n-th frame can be expressed aswhere aSv^ represents the amplitude of the useful signal component, p(v)[fc, s] ± p(kTs- 4V)- sT0), represents the s-th transmitted pulse, and n, k, ip^[s, n], w^[k,s, n] have the same meaning as the corresponding terms illustrated for (6). The pre-processing steps consist of a) matched filtering (see Figure 2d):b) subsequent combination of pulses belonging to the same frame byc) subsequent target detection and d) extraction of relevant spectral information as described for the FMCW radar case.It is important to observe that, unlike the FMCW case, the target detection step does not require calculation of aDFT, but only the calculation of the index associated with the target (via (19)). Indeed, the signal obtained after matched filtering and combination (x'^ [k, n]) provides information about the target range and can be easily expressed asThe sequences {s^ [n]; n = 0,1, ... , NF— 1} feed the algorithm.

[0070] The description of the proposed radar systems and adopted models deserves several comments, which are listed below.Multipath attenuation in VSM

[0071] The point target model illustrated above is commonly adopted in the technical literature on radar-based VSM. Unfortunately, this model does not take into account the fact that: a) the received signal is actually the superposition of multiple echoes, originating not only from the body of the monitored subj ect but also from any detectable target in the considered propagation environment (usually, an indoor environment); b) the human (or animal) chest, being quite close to the radar device (its typical distance is on the order of 1-2 meters), behaves as an extended target and consequently generates multiple echoes that may interfere constructively or destructively on each of the RX antennas. For these reasons, the communication channel characterizing the employed radar system is affected by multipath fading; this phenomenon can completely disturb the (weak) traces of respiration and heartbeat in radar-based VSM.

[0072] These considerations have allowed us to interpret experimental results; indeed, such results have shown that: a) The quality of signals acquired through the different virtual antennas (VAs) of a given radar device in the same observation interval can be substantially different

[0026] , In particular, as shown in the paragraph Impact of multipath on radar measurements, some of these signals are useless since their amplitude spectra do not exhibit significant peaks; consequently, the expected traces of respiration and heartbeat cannot be easily identified. b) The signal quality provided by a given virtual antenna varies over time. In other words, a virtual antenna that provides relevant spectral information in a certain observation interval may become useless in a different interval.

[0073] The impact of multipath fading can be mitigated by exploiting the spatial diversity provided by a single radar equipped with a reasonably large array; in this case, beamforming techniques can be employed [13, 27, 28], In the present work, a different strategy, simpler from the signal processing point of view and suitable for low-cost implementation, has been chosen. Indeed, diversity offered by multiple virtual antennas and optionally by multiple radars (equipped with any number of antennas) is exploited.

[0074] A solution using multiple radars was originally proposed by Rong et al.

[0029] to compensate human body movements in VSM through a pair of IR-UWB radars. In

[0030] , the same pair was employed for VSM of a human subject performing large movements, while in

[0031] , two IR-UWB radars with different central frequencies wereused for estimating vital signs of multiple humans.

[0075] In the present work, the NDradars may be equipped with a single TX and a single RX (in this case ND> 1) or with RX and TX arrays. The presence of multiple antennas (1VTOT> 1) allows adopting a selection combining strategy. For this reason, as illustrated in the following section, in each superframe the algorithm selects, among the ZVT0Tsequences {s^ [n]; n = 0,1, ... , NF— 1], the / VT0Tsequences that offer the spectrally richest information for vital signs detection.Exploitation of prior information on the position of the monitored subject

[0076] If the only target in the radar’ s field of view was the subj ecf s chest, Zmin= 0 and Zmax= No— 1 could be selected in (20). Unfortunately, in a real scenario, the presence of multiple targets gives rise to multiple spectral peaks, some of which may be stronger than that associated with the subj ecf s chest.When this happens, (19) leads to selecting the wrong bin, thereby losing vital signs information. This is illustrated in Figure 5, where the range profile from a single VA and for ZVF= 317 is shown. Here, the portions indicated by reference number 5001 refer to the strongest targets located within the radar’ s field of view (FOV). In this scenario, the range bin at 5.5 m (corresponding to the position of the wall) would be erroneously selected as the patient’s chest, which instead is located at about 1 m and is represented by the horizontal dashed line in the figure and indicated by reference number 5003. The range limits, represented by continuous horizontal lines indicated by reference number 5002, are set to prevent the wall (recognizable as the strongest target at 5.5 m) from being selected by ( 19) as the subj ecf s chest.

[0077] In the case of FMCW radar, if the position of the subject’s chest is known to belong to the interval [Z?min, Z?max], Zminand Zmax, i.e., the limits of the search interval, can be evaluated as and

[0078] respectively. Here again, = B / T is the chirp slope, Nois the total number of samples or intervals (order of the DFT applied to each frame), Tsis the sampling interval, and c is the speed of light.

[0079] In Figure 5, the introduction of range limitations at Z?min= 0.5 m and Rmax= 1.5 m, represented by continuous horizontal lines indicated by reference number 5002, allows correctly selecting the patient’s chest through ( 19)-(20).Data block collection

[0080] In fact, HR and BR are time-varying parameters and, in V SM, their fluctuations must be carefully tracked;for this purpose, it is necessary to exploit the correlation existing between observed values of these parameters in adjacent time intervals. For this reason, a buffering mechanism must be implemented to discard frames generated too far back in time.

[0081] A practical implementation of the aforementioned buffering consists of the following approach: the superframe stream provided by the NDradar devices is divided into partially overlapping blocks (called data blocks), each composed of NBconsecutive superframes (and consequently of duration NBTSFseconds), as shown in Figure 6; note that the start of each data block is shifted by Nssuperframe intervals relative to the start of the previous block, with Ns< NB. In this case, from the available superframe sequence, four partially overlapping blocks are extracted, each collecting measurements acquired over NB= 8 consecutive superframes, and the block shift / Vsis equal to 3.

[0082] New HR and BR estimates are evaluated on each data block and, consequently, are made available every NsTSFseconds; however, as illustrated in the following section, they are processed by a filtering algorithm to exploit the above correlation.

[0083] It is important to note that the duration of each data block must be selected so as to obtain sufficiently accurate BR estimates. Note that respiration, being a slower phenomenon than heart beating, requires a longer observation interval for its estimation.

[0084] Experiments conducted by the inventors have shown that the duration of each data block should cover at least two consecutive breathing periods; this means NBmust be selected so that the inequality NBTSF> 2 Tb minholds.DESCRIPTION OF THE PROCESSING ALGORITHM

[0085] The internal structure of the processing algorithm is described by the block diagram shown in Figure 7a. In this section, the signal processing performed by the blocks appearing in this scheme is analyzed. The description shows how the p-th data block is processed to estimate both HR and BR and how the estimates extracted from consecutive data blocks are filtered to improve overall accuracy.I) Phase gradient cateiilation (PG)

[0086] The first block of Figure 7a is fed by the setv = 0,1, ... , NT0T— 1} which collects all the measurements forming the p-th data block, whererepresents the contribution of the v-th virtual antenna; here

[0087] For any v, the (1VB— J ri) -dimensional vectorgapwhich conveys information on phase variations (i.e., on PG) observed on consecutive samples of the sequence {SpV')[n]}, is evaluated.

[0088] The n-th element of the vector gg (26) is evaluated aswhereis the n-th element of the so-called PG sequence

[0016] ; here, rin is a positive integer parameter called time lag. For this parameter, a small value is chosen to avoid the phase wrapping phenomenon occurring when the phase angle of a measured signal exceeds the range —n to n (or 0 to 2TT) and is thus "wrapped" within this interval. This can lead to misinterpretation of phase data, as phases exceeding n are brought back into the principal range, causing apparent discontinuities in the data.

[0089] Another implementation of the first block of Figure 7a consists in calculating the (1VB— 2) -dimensional vector gapwhere the n-th element, in this case, is evaluated aswhere(v) ( , O)\Pp = unwrap I ZSp 1 and the “unwrap” operator refers to the operation of adding or subtracting multiples of 2TT when the modulus of the difference between two consecutive samples exceeds TT (htps: / / en.vri&pefe^phase . and frequency). The above implementation allows obtaining a phase gradient through the second derivative of the phase. It is then possible to realize the phase gradient step via differentiation of any order and with any implementation variant.2) Spectral analysis and sparsifkation

[0090] In the second block of Figure 7a, the spectral content of each of the 1VTOTvectors {gg f} (26) is analyzed, and a sparse representation thereof is generated. These tasks can be performed by the SFEC algorithm illustrated in [32, Par. H.B]. This algorithm identifies, for any v, the dominant sinusoidal components of the sequence {g^p [n] } and estimates their parameters (i.e., the normalized frequency and the complex amplitude).

[0091] In practice, for any v, the vectoris used to generate the (1VB— J n) -dimensional vectorwherewith n = 0,1, ... , NB— An — 1 and q = 1,2. Then, for any v, the vectors {g^,zero-padded and subj ected to a discrete Fourier transform (DFT) of order No' \ here, Nf = MfNB— An) and M ' represents the selected oversampling factor. In this way, the No'-dimensional vectorwith v = 0,1, ...,NTOT— l and q = 0,1,2, is obtained; here,with I = 0,1, ... , No' — 1. Then, for any v, the set {G^; q = 0,1,2], which collects three different spectra, is processed by an algorithm that, in each of its iterations, identifies the dominant tone, estimates its parameters, and removes its contribution from the spectrum.

[0092] In this case, prior knowledge on the typical frequency ranges of HR and BR is exploited. For this reason, the spectral analysis and sparsification algorithm (for example, SFEC) is executed twice.

[0093] More precisely, in the first execution, it is employed to detect the Mbstrongest tones falling within the respiration frequency interval and to estimate their parameters. In this way, tone detection, which is based on the well-known periodogram method (see, e.g.,

[0033] ), is limited to the frequency range [ / b,min< / b.max] that includes typical BRs.

[0094] The second execution of the algorithm aims at detecting and estimating the Mhstrongest tones in a different frequency range, in particular in the interval [ / h.min > 2 / h,max] , which contains both the first harmonic (FH) and the second harmonic (SH) of the cardiac activity. This last choice is motivated by the fact that, as highlighted by experimental results of the inventors (see below), the SH of the heart can be significantly stronger than the associated FH; this phenomenon is unpredictable and clearly highlights the limits of the adopted model for chest displacement representation (see (3)-(5)).

[0095] The outputs of the spectral analysis and sparsification block are the Mx-dimensional vectors (see [18, Sec. m.B]) fr’W = pty [p]. tyFand[pi . 4i-i [pi] ■ (34) which collect the normalized frequencies and the corresponding complex amplitudes, respectively, of the Mxstrongest tones potentially related to respiration (if x = b) and to the heart (if x = h) detected on the v-thVA (with v = 0,1, ... , ATOT— 1). The normalized frequency F associated with frequency / is defined as F = / TSh / / Vo', where TSFand 7V0' represent, respectively, the superframe period and the DFT (or IDFT) order adopted by the spectral analysis and sparsification algorithm.

[0096] Thus, for any v, the vectorconsisting ofthe amplitudes ofthe same tones, is generated for x = b and ft; here, [p] | form =0,1, ... , MX- 1.

[0097] Note that the sets {Fx^ Cp]; v = 0,1, ... , NT0T— IJ andfa^ tp]; v = 0,1, ... , 1VTOT— 1} provide, for the v-th VA, a sparse representation of the spectral content potentially originating from respiration (if x = b) and from the heart (if x = h) during the p-th data block.3) Antenna selection

[0098] In the third block of Figure 7a, the quality of the spectral information provided by the v-th VA is evaluated based on the vector a^ tp] (35) (with v = 0,1, ... , 1VTOT— 1).Thus, only the best VAs are selected for subsequent processing. More precisely, the antenna selection strategy adopted in this exemplary embodiment consists of a) evaluating the quality index2 2Q'*” w = s" o-1(4? w) + V.V(4? w) ■ (36) for the v -th VA (with v = 0,1, ... , NT0T— 1); b) selecting the VAs associated with the 1VTOThighest values thereof.

[0099] Note that (36) corresponds to a possible embodiment of the quality index calculation. Furthermore, note that Q [p] (36) is proportional to the overall energy associated with the (Mb+ Mh) tones detected by the SFEC algorithm, and the choice of this index is motivated by our experimental results. Indeed, our measurements have shown that low-quality spectra (i.e., those in which spectral peaks related to cardiac and respiratory activities cannot be easily identified

[0034] ) are usually associated with low values of the adopted index. Further details are provided in Key insights derived from experimental results.

[0100] The spectral information available after antenna selection for the p-th data block are collected in the setstty, ... , t^TOT-i) (37) andSb[plwhich refer, respectively, to the heart and respiration; here, vkindicates the index of the fc-th selected VA.

[0101] Processing of the set <Sh[p] (37) is performed for the estimation of HR, and together with the set <Sb[p] (38) for apnea detection. For the estimation of BR, only the set 5b[p] (38) is used. These tasks are executed in the remaining three blocks of Figure 7a.4) Apnea detection

[0102] The sets [p] (37) and 5b[p] (38) are passed to the fourth block of Figure 7a. This block is used to detect apnea events, characterized by the disappearance of the breathing peak in the spectral region of interest. The output is a binary variable, denoted ( [p] , whose value is set to 1 (0) if an apnea event is detected (not detected) in the p-th data block. The strategy adopted in the present work for generating % [p] can be expressed as^[p]A(1 ifah,max [p] >ab ,max [p] (39)10 otherwise where <4,max[p] and ab,max[p] represent the largest elements of all vectors in the sets {a^ fp]; v = v0, 14,, respectively.

[0103] In other words, apnea event determination is based on the comparison between displacement information of the monitored subject’s body related to vectors representing displacements associated with cardiac pulsations [p] ) and displacements associated with thoracic expansions due to respiratory acts.

[0104] The underlying hypothesis of this criterion is that, when physiologically present, thoracic displacements due to respiratory acts exhibit peaks always greater than those related to cardiac activity.5) Heart rate estimation

[0105] The set <S"h[p] (37) feeds the fifth block of Figure 7a, whose task is to estimate heart rate.

[0106] In the development of an algorithm for this purpose, the following specific aspects, emerging from spectral analysis of the measurements, have been considered: a) Cardiac harmonic intensities are usually much weaker than those due to respiration. Furthermore, the strongest spectral peak observed in the typical cardiac frequency range does not necessarily represent the cardiac contribution, but may originate from multipath or noise contributions. b) The amplitude of the cardiac second harmonic (SH) can be significantly stronger than the associated first harmonic (FH). This phenomenon, first observed by Park et al.

[0035] , motivates the approach of exploiting both spectral contributions to obtain a reliable HR estimate. c) The intensities of cardiac FH and SH are random and vary over time.

[0107] These aspects make HR estimation a very difficult problem. In particular, referring to

[0016] , the technology described therein is not usable for HR estimation without deep modifications or improvements, which are also the subject of the present invention.The solution developed here, based on findings of

[0016] , exploits the correlation of multiple consecutively estimated HR values observed over consecutive data blocks and is based on a Bayesian filtering approach (see e.g.

[0036] ). In particular, one step of the method involves employing, for each data block defined in preprocessing, an iterative filter with memory on the parameter vector to generate a synthesized version of the parameter vector, at least partially cleaned from random body movements and multipath contributions.

[0108] Indeed, the estimation method of the present invention is structured according to two distinct phases: a time update and a measurement update, similarly to the general formulation of Bayesian filters. However, unlike known approaches (such as the Kalman filter and variants), which require known statistical models for the state transition function and measurement model, the proposed method does not rely on prior knowledge of such probabilistic models.The filter structure is maintained, but the implementation of the two phases is performed via a procedure based on state estimates at previous steps, which allows estimating the temporal evolution of HR based on observed data. Consequently, the method is applicable despite unknown statistical models, thus overcoming operational limitations of traditional filtering techniques.

[0109] In the present work, since both models were unavailable, the following heuristic approach was adopted: firstly, HR for each data block was represented as a discrete state variable, able to assume Nhdistinct integer values. Note that, since HR is discretized, the accuracy of this parameter estimate is limited by the selected discretization step, i.e., practically by the spectral resolution adopted.

[0110] Specifically, it is assumed that the values of this state variable belong to the setwhose elements represent the indices of the Nhfrequency bins resulting from the DFT processing of order No' (see (32) with q = 0) falling within the frequency interval [ h,min< 2 h,max] •

[0111] A simple algorithm was designed to update the probability that HR takes a specific value in c / Zh(i.e., falls within a specific frequency interval) in a given data block; this algorithm, representing the time update of the filtering method, can be summarized as follows.

[0112] Let a® [p — 1] be the probability that HR belongs to the n-th frequency interval during the (p — 1) -th data block, with n G c / Zh; then, the prediction[p] of the probability that HR falls within the m-th bin in the subsequent (i.e., the p-th) data block is calculated aswith m = n^n, n®n+ 1, . . . , n®n+ Ah— 1; here, denotes the probability of a state transition from theZ-th bin to the m-th bin, and ® is the set collecting indices of bins in the (p — l)-th data block from which it is possible to reach the m-th bin in the p-th data block.

[0113] In this work, the following specific choices were made:a) The setwas chosen,; here,andare the indexes of the minimum and maximum bins, respectively, among the bins in the (p — l)-th data block from which the m-th bin of the p-th data block can be reached. Lhis a positive integer defining the maximum size of the set J® (which, at most, consists of 2Lh+ 1 consecutive integers). b) The state transition probability / ?was evaluated aswhere JVtyx; p, a2) denotes the probability density function of a Gaussian random variable with mean p and vari •ance 2 .

[0114] On the other hand, the measurement update of this filtering method is expressed asfor n = n^n, n®n+ 1, . . . , n®n+ Nh— 1; here, y^h')[p] denotes the probability that, based on the SFEC output available for the p-th data block (i.e., the measurements for the considered update), HR falls within the n-th frequency interval. In practice, this probability is evaluated aswith n = n®n, n®n+ 1 n®n+ lVh- 1. Here,is the probability density function describing a Gaussian mixture (GM) with Mh= 2MhNT0Tcomponents, and p^ [p], and [p] denote, respectively, the weight (positive), mean, and standard deviation of its i-th component (note that the GM weights must satisfy the normalization condition= 1).

[0115] In the present algorithm, the parameters of fph>(x) (48) are evaluated as follows. First, the Mh-dimensional vector of absolute amplitudesand the Mh-dimensional vector of normalized frequenciestp],0.5 F<’“> [p],0.5 F<’*> [p] 0.5 F, / -'"’1’ [p]], (50) which collect respectively all normalized frequencies estimated by SFEC within the considered frequency interval and their amplitude values divided by two. Thus, the mean and variance of the i -th component of(x) (48) are evaluated asrespectively, with i = 1,2, . . . , Mh; here,[p] and a® [p] indicate the i-th element of Fh[p] (50) and ah[p] (49), respectively, a®x[p] — niax ah[p], d^ fp] = min ah[p], ahis a positive scale factor, and p is a positive parameter (called the interval control parameter) ensuring that the range of the denominator in the fraction ontheRHS of (52) is the interval [p, 1] (with 0 < p < 1 selected).

[0116] Once the measurement update is performed, the most probable state index (i.e., frequency bin) is calculated as n<h) [p] = argand the estimateof HR is evaluated based on this index.

[0117] The filtering algorithm is initialized by settingfor n6) Breathing rate estimation

[0118] The objective of the sixth block of Figure 7a is BR estimation. The developed algorithm is fed by the set5b[p] (38) and optionally by the apnea indicator ( [p] (39), a binary variable (0 or 1) indicating the presence or absence of apnea, here used to enable or disable steps of this algorithm. Its structure is similar to that described above for HR estimation; the main differences with respect to the previous algorithm are summarized below. It is worth noting that apnea determination is not fundamental to determining other parameters but can be advantageously used to avoid pointwise and real-time determination of the same when it is already known that there is an absence of respiratory acts.

[0119] Furthermore, apnea detection is faster than HR and BR determination and enables faster alerting of medical personnel about the patient's health risk. Main differences compared to the previous algorithm are: a) BR for each data block is represented as a discrete state variable that can assume Nbdistinct integer values, particularly any value in the set c / Zb= n®n+ 1, . . . , n®n+ Nb— 1 j; such values representing theindices of the Nbfrequency bins falling within the frequency interval [ / b,min > / b,max] • b) The vector of amplitudes ah[p] (49) and frequencies Fb[p] (50) are replaced by the corresponding MblVT0T- dimensional vector of absolute amplitudesand by the Mb / VT0T-dimensional vector of normalized frequenciesrespectively. Note that the FH of respiration is always stronger than the associated SH; for this reason, the spectral contribution of the latter harmonic is not considered. c) The set J® (42) and the state transition probability(45) are replaced byandrespectively; here,and Lbis a positive integer defining the maximum size of the set J® (which, at most, consists of 2Lb+ 1 consecutive integers). Note that, generally, the variance <Jband parameter Lbdiffer from their counterparts <Jband Lbin (43)-(44). This is due to the fact that, in general, the variability of breathing observed for agiven subject alwaysdiffers from cardiac variability. Moreover, HR is outside the subject’ s control, whereas BR is not. d) The probability y^ [p] that BR belongs to the n-th frequency interval is evaluated asTnb)[p] = fpb\n), (62) whereis the probability density function describing a GM with Mb= MbNT0Tcomponents, and, p® [p], and cr® [p] denote respectively the positive weight, mean, and standard deviation of its i -th component (note that GM weights satisfy the normalization= 1).The mean p ■b-)[p] and variance [p] of the i -th component of(x) are evaluated similarly to (51) and(52), respectively. e) The algorithm is initialized by settingfornf) The estimate of BR is forced to 0 Hz throughout an apnea event. In this case, the variable % [p] (39) is 1 and the BR estimation algorithm halts. When the apnea event ends, the BR estimation algorithm is reinitialized (see (64)).

[0120] This concludes the description of the processing algorithm, whose steps are summarized as follows: Input: Integer parameters NB, Ns, An, N0', Mb, Mb, real parameters ab, ab, T], and sets {s^; v = 0,1, , NTOT- 1} (see (1)), c / ty (2), c / fy, {mt(h)}, {mt(b)}, {£® }, and {^} (see (3) and (4), respectively). Initialization: Initialize[—1]} according to (5) and (6), respectively.Main: For p = 0 to (NF— NB~) / Nsdo a. Phase gradient calculation (PG): Calculate elements of each vector in the set {gg ?} (7) based on (8)-(9). b. Spectral analysis and sparsification: Execute the SFEC algorithm to generate sets of vectors {F^v)[p]} (10) and {A^ [p]} (11), with x = b and h. This requires calculating vectorsaccording to (14) and (15), respectively. c. Antenna selection: Generate set {Q(v>[p]} based on (16). d. Apnea detection: Calculate parameter % [p] according to (17).Ilfylp l = 1 thenSet fb[p] = 0 and compute} based on (6). else: e. Respiratory frequency filtering- Time update: Evaluate probabilities {a^ [p]} (see (18)).f. Respiratory frequency filtering - Measurement update: Calculate meansand variances[p]) } (see(19) and (20), respectively). Then calculate {y^b')[p]} (21) and {cr„b')[p]} (see (22) and (23), respectively). g. Respiratory rate estimation: Calculate the BR estimate fb[p] according to (24)-(25). end h. Cardiac frequency filtering - Time update: Calculate probabilities {a^ [p] } based on (18). i. Cardiac frequency filtering - Measurement update: Calculate means {p • [p] } and variances { [p] ) }based on (19) and (20), respectively. Then calculate {y^h')[p]} (26) and[p] } according to (22) and (23), respectively. j . Heart rate estimation: Calculate HR estimates fb[p] according to (24)-(25).Output: The estimates { / h[p]; p = 0,1, . . . , (NF— NB) / NS} and { / b[p]; p = 0,1, . . . , (NF— N^ / N^.COMMENTS ON THE PROCESSING ALGORITHM

[0121] The complete algorithm, including the preprocessing steps described in the DETAILED DESCRIPTION section and the processing steps described in the DESCRIPTION OF THE PROCESSING ALGORITHM section, and considered for 1VTOTvirtual antennas derived from the use of NDradar boards, is shown in Figure 7b. The block diagram shown in Figure 7c, instead, represents in greater detail the complete algorithm in the case of FMCW radar and 1VTOTvirtual antennas. This algorithm deserves several comments, which are listed below.1) Calculation of the phase gradient

[0122] As already mentioned previously, the sequence {p^ W } (^ (27)) conveys information about the phase gradient (PG) characterizing the sequence {s^ [n]} and referring to its elements spaced by An. Further details on the main properties of the PG are provided in Appendix B Analysis of the Phase Gradient, where it is demonstrated that: a) the calculation of the sequence {p^ } simplifies the identification of the first harmonic (FH) of the heart and respiration, enhancing their spectral contribution; b) the selection of the imaginary part of {f^ [n]} (27) aims to suppress both the undesired even-order harmonics (which appear in its real part) and the DC component (see also [16, Sec. in, eq. (12)]).2) Spectral analysis and sparsifkation

[0123] In the present work, the selected value for the parameter Mbis always greater than that chosen for Mb. This is due to the fact that the FH of respiration is usually characterized by a high signal-to-noise ratio (SNR).3) Antenna selection

[0124] The antenna selection criterion adopted (see (36)) is not influenced by the specific characteristics of the measurement scenario (for example, by the propagation environment, the monitored patient’s characteristics, and the radar positioning).

[0125] It should also be noted that the sets <Sh[p] (37) and 5b[p] (38) do not convey phase information. This is because the processing performed by the HR and BR estimators does not require coherently combining the contributions from different virtual antennas.4) Apnea detection

[0126] The apnea detection criterion is very simple and exploits the information provided by the 1VTOTselected virtual antennas in a non-coherent manner.5) Estimation of heart rate and respiratory rate

[0127] The development of algorithms for HR and BR estimation required several assumptions. It is important to emphasize that: a) The models (3^L(45) and (3^L(59) for the state transition probabilities were selected based on our experimental results; further details on this topic are provided in Paragraph 5.2. b) The sets ® (42) and ® (58) and the state transition probabilities (3^ (45) and (3^ (59) do not evolve over time. This assumption is reasonable when vital signs are estimated for a resting patient. c) The formula (52) adopted for the evaluation of the variance[p] of the i -th component of the Gaussian mixture (GM) (48) is motivated by the fact that a reduced spectral peak is associated with a low SNR and, consequently, with a large uncertainty; similar considerations apply to the evaluation of the variance [p] ■

[0128] The implications of this choice are exemplified by Figures 8a and 8b. In particular, the first figure shows the amplitude spectrum of the PG for a specific virtual antenna selected by the processing algorithm and the Mh= 5 peaks detected by the SFEC algorithm in the HR range, indicated by solid vertical lines and reference numeral 8001. The second figure, instead, shows the GM generated based on the processing algorithm output under the assumption that its components have equal weight; note that the weaker spectral peaks appearing in the first figure, indicated by reference numeral 8003, lead to broader Gaussian components, indicated by reference numeral 8004.

[0129] Further information on the considered example is provided in Table 1, listing the normalized frequency and amplitude of each spectral peak, as well as the mean and variance of the 2Mh= 10 components of the associated GM.

[0130] The following table reports the values of normalized frequency and amplitude characterizing the five spectral peaks highlighted in Figure 8a. Also listed are the means {p^} and variances { } of the ten GMcomponents generated based on these results and equations (51)-(52).NUMERICAL RESULTSThe objective of this section is threefold. First, the experimental apparatus developed for acquiring vital sign information from the radars and from the reference sensor is described. Subsequently, some important insights drawn from the analysis of the experimental results are commented on.1) Experimental configuration

[0131] Performances of the considered VSM algorithms were evaluated based on measurements acquired through a pair of Position2Go radars

[0041] on a human subject. These devices, produced by Infineon, are FMCW radars operating at the frequency / o = 24 GHz, and are equipped with a single TX antenna ( / VTX= 1) and a pair of RX antennas (IVRX = 2); therefore, in this experiment, VI OT= 4.

[0132] In the present work, the following values were selected for the parameters of each radar device:1) Chirp slope g = 0.78 MHz / / zs;2) Total number of samples per chirp N = 256;3) ADC sampling frequency fs= 1 MHz;4) Chirp duration To= 456 gs;5) Frame duration TF= 0.063 s;6) Number of chirps per frame Nc= 3.

[0133] Note that, in this case, the ramp time is equal to T = N / fs= 256 gs and the bandwidth of the radiated signal is B = 200 MHz.

[0134] The two radar devices were mounted on a support, in a non-planar configuration (the non-planarity of theradar boards used allows beter exploitation of spatial diversity, thus attenuating the impact of the multipath phenomenon), at about 1 m in front of the subject to be monitored, as illustrated in Figure 9. In the case of the proposed algorithm, these devices operate in round-robin mode. Indeed, the first chirp of each frame radiated by the second radar is activated after the transmission of the last chirp of the frame sent by the first radar has ended. In this way, it is ensured that the chirps of each radar are transmitted and received during the idle time (of duration T, ) of the other; this is exemplified in Figure 4.

[0135] The processing algorithm selects two of the four signals provided by all the virtual antennas (i.e., 1VTOT= 2). In the known art, instead, the measurements to be processed come from a single radar device (i.e., from two virtual antennas), which is selected at the beginning of data acquisition and is never changed; in that case, to combine the signal of the two virtual antennas, beamforming is used.

[0136] To acquire an accurate estimate of vital signs, a reference sensor was also used, the Shimmer3 ECG Unit

[0042] , This FDA-approved device, produced by Shimmer Research Ltd, uses a Bluetooth radio interface to transmit acquired measurements and is equipped with five electrodes, which must be placed at specific positions on the chest of the examined subject.

[0137] Figure 9 shows a representation of the adopted measurement configuration: two Position2Go radars, indicated by reference numeral 9001, installed on a support, are placed in front of the subject to be monitored. The Shimmer3 ECG unit, indicated by reference numeral 9002, is used as the reference sensor.

[0138] In the experimental acquisitions, the reference ECG signal was read from the voltage difference between the left leg electrode and the right mid-axillary electrode (LL-RA). The respiratory signal, instead, was generated by measuring the impedance between these two electrodes. The time-domain signals transmitted via the Bluetooth interface were subjected to FFT processing followed by peak selection for BR and HR estimation (unfortunately, apnea events cannot be detected).

[0139] Thanks to the availability of a dedicated API, the Shimmer3 measurements can be easily synchronized in time with those acquired through the radar devices, so that data acquisition can occur in real time. Indeed, the Shimmer3 device can share its time reference (i.e., the CPU timestamp) with the radar devices; this allowed us to synchronize the reference sensor and the radars with high precision.

[0140] The conducted experimental campaign made it possible to acquire the required radar signals and associated reference signals for multiple acquisitions and in two different scenarios: a) a single subject in conditions of absence of apnea; each acquisition lasts on average 70 seconds, for a total duration of about 610 seconds; and b) a single subject experiencing apnea events of arbitrary duration; each acquisition lasts on average 250 seconds, for a total duration of about 2000 seconds, and a single acquisition contains multiple apnea events.The chosen values for the parameters of the considered estimation algorithms are listed in the following Table 2.2) Main insights derived from the experiments! results

[0141] A thorough analysis of our experimental results highlighted the importance of various aspects to be carefully considered in radar-based vital sign estimation. This subsection focuses on four specific issues deemed particularly significant in VSM. a) On the importance of the second harmonic in heart rate estimation

[0142] As already mentioned in the DESCRIPTION OF THE PROCESSING ALGORITHM, spectral analysis of the acquisitions showed that the second harmonic (SH) of the heart rate is sometimes stronger than the associated FH. This heart rate doubling effect [5] is clearly visible in Figure 10, where, for a given acquisition, the amplitude spectrum observed on a single virtual antenna of each radar device is shown.

[0143] In this case, any strategy based on spectral peak detection inevitably leads to an erroneous HR estimate. The proposed algorithm, instead, exploits the presence of the second harmonic to improve the reliability of HR estimates.In particular, Figure 10 shows a representation of the amplitude spectrum observed on a single virtual antenna ofeach radar device used by the proposed algorithm. The SH of the heart displacement is identified by the vertical line indicated by reference 10001 corresponding to a frequency of 1.72 Hz; its strength is significantly greater than that of the associated FH, which is identified by the vertical line indicated by reference 10002 corresponding to a frequency of 0.86 Hz. b) Impact of multipath on radar measurements

[0144] Due to the multipath phenomenon and the short wavelength characterizing the radar devices used, significant variations are observed in the quality of the signals received on different antennas

[0026] , This is exemplified by Figure 11, where the amplitude spectra observed on the two VA of the same radar device are shown.

[0145] Figure 11 shows a representation of the amplitude spectra observed on the two VA of a single radar. The vertical lines indicated by reference numerals 11001 and 11002 identify, respectively, the BR and HR provided by the Shimmer3 device. The amplitude spectrum evaluated for the first virtual antenna (dashed curve) shows well- defined peaks centered at the heart and respiratory frequencies; the other amplitude spectrum (solid curve), instead, does not have this property and is deeply attenuated in the frequency range of interest.

[0146] In this case, in fact, the amplitude spectrum referring to the virtual antenna of the first radar is characterized by clearly visible peaks corresponding to the FH of respiration and heart, while that referring to the virtual antenna of the second radar is useless, as it shows no significant peaks in the frequency range of interest.

[0147] The same concept was exploited also in

[0034] to select the range bins containing relevant information on vital signs; in that work, the indicator used was the ratio between the in-band energy (i.e., the energy obtained by accumulating frequency responses within the respiration and heart rate intervals) and out-of-band energy.

[0148] It should also be noted that the qualify metric (36) Q(1)= 5.06 x 104obtained for the first radar is significantly greater than its counterpart= 6.56 x 10-6evaluated for the other radar. c) Apnea detection performance

[0149] As previously illustrated, the processing algorithm detects an apnea event when the most significant peak observed in the heart region is higher than its counterpart in the respiration region. This is exemplified by Figure 12, where the spectrum observed before an apnea event and that obtained during the event on the same virtual antenna are represented by the dashed and solid curves, respectively, indicated by references 12001 and 12002. Note that the position of the peak corresponding to the HR, indicated by reference numeral 12003, changes during this event (in particular, a physiological increase in HR is observed).

[0150] Experimental results also showed that: a) the missed detection probability (corresponding to the case where apnea is not detected during the entire event, nor immediately after its conclusion) for apnea events shorter than 15 seconds is about 0.15; otherwise, it is zero; b) false alarm events are very unlikely (in fact, no such events were observed in our experiments); c) a delay ranging between 3 s and 15 s (with an average of 10 s) is observed in apnea detection.

[0151] Figure 12 shows a representation of the effects of apnea on the amplitude spectrum observed on a given virtual antenna. The spectrum observed during apnea (solid curve), contrary to that referring to its absence (dashed curve), is characterized by the absence of the BR peak (whose position is identified by the vertical line indicated byreference numeral 12004). d) Variability of heart rate and respiratory rate in resting conditions

[0152] The HR and BR signals acquired by the reference sensor (Shimmer3) during the entire measurement campaign were used to estimate the values of the following parameters:1) cardiac state transition variance <Jb, used in evaluating the cardiac state transition probability(45);2) respiratory state transition variance <Jb, used in evaluating the respiratory state transition probability (59);3) maximum cardiac deviation Lh, i.e., the maximum size of the set ®4) maximum respiratory deviation Lb, i.e., the maximum size of the set

[0153] To estimate the state transition variances, the Matlab "Distrib" tool from the “Statistics and Machine Learning 12.4” toolbox was used, which provided <J = 1.75 x 10-3and <Jb= 2.60 x 10-3. Furthermore, our results showed that, for TSF= 0.063 s, Lh= 12 and Lb= 18 should be selected for the maximum size of the sets ® and ®, respectively. e) Sensitivity to the position of the monitored subject

[0154] The experimental results of the inventors have shown that the accuracy of the processing algorithm is marginally influenced by the position of the monitored subj ecf s torso, provided that the subj ect remains within the radar’ s FOV. It should also be noted that this is no longer true if beamforming is used in combination with known methods.

[0155] Throughout the description, similar reference numerals may be used to identify similar elements.

[0156] While various aspects of the embodiments are presented in the drawings, the drawings are not necessarily to scale unless specifically indicated.

[0157] It will be readily understood that the components of the embodiments as generally described herein and at least partially illustrated in the accompanying figures may be arranged and designed in a wide variety of different configurations.

[0158] Therefore, the following description, as also represented in the figures, is not intended to limit the scope of the present invention but is simply representative of some possible embodiments.

[0159] The present invention may be embodied in other specific forms without departing from its spirit or essential characteristics.

[0160] The embodiments described should be considered in all respects only as illustrative and not restrictive.

[0161] All changes that fall within the meaning and range of equivalency of the claims are to be encompassed withintheir scope.

[0162] Reference throughout this description to features, advantages, or similar language does not imply that all features and advantages that may be realized with the present invention should be or are present in any single embodiment of the invention.

[0163] Rather, the language referring to features and advantages is intended to mean that a feature, advantage, or particular characteristic described in connection with one embodiment is included in at least one embodiment of thepresent invention.

[0164] Consequently, discussions of features and advantages and similar language throughout this description may, but do not necessarily, refer to the same embodiment.

[0165] Furthermore, the features, advantages, and characteristics described of the invention may be combined in any suitable manner in one or more embodiments.

[0166] Those skilled in the art will recognize, in light of the present description, that the invention may be practiced without one or more of the specific features or advantages of a particular embodiment.

[0167] In other instances, additional features and advantages may be recognized in certain embodiments that may not be present in all embodiments of the invention.

[0168] Reference throughout this description to “an embodiment,” “a variant,” “an exemplary form,” or similar language means that a particular function, structure, or characteristic described in connection with the indicated embodiment is included in at least one embodiment of the present invention.

[0169] Consequently, references to “an embodiment,” “a variant,” and similar language throughout this description may, but do not necessarily, all refer to the same embodiment.

[0170] The components of the embodiments as generally described herein and illustrated in the accompanying figures may be utilized and designed in a wide variety of different configurations.

[0171] All changes falling within the meaning and scope of equivalence of the claims are to be encompassed withintheir scope.

[0172] Although the operations of the methods herein are shown and described in a particular order, the order of operations of each method may be altered such that some operations may be performed in a different order and / or such that some operations may be performed, at least partially, simultaneously with other operations.

[0173] Note that at least some of the operations for the methods may be implemented using software instructions stored on a computer-usable storage medium for execution by a computer.

[0174] As an example, one embodiment of a computer program product includes a computer-usable storage medium to store a computer-readable program that, when executed on a computer, causes the computer to perform operations as described herein.

[0175] Further, embodiments of at least portions of the invention may take the form of a computer program product accessible from or readable by a computer-usable or computer-readable medium providing program code for use by or in connection with a computer or any instruction execution system.

[0176] For the purposes of this description, a computer-usable or computer-readable medium may be any apparatus that can contain, store, communicate, propagate, or transport the program for use by or in connection with the system, apparatus, or device for instruction execution.

[0177] The computer-usable or computer-readable medium may be an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system (or apparatus or device), or a propagation medium.

[0178] Examples of a computer-readable medium include semiconductor memory or solid-state memory, magnetic tape, a removable computer diskette, random access memory (RAM), read-only memory (ROM),magnetic hard disk, and optical disk.

[0179] Current examples of optical disks include compact disk read-only memory (CD-ROM), compact disk read / write (CD-R / W), digital video disk (DVD), and Blu-ray disk.

[0180] It is, however, evident that the invention should not be considered limited to the particular arrangements illustrated above, which constitute only exemplary embodiments thereof, but that various modifications are possible, all within the reach of one skilled in the art, without thereby departing from the scope of protection of the invention itself, which is defined by the following claims.REFERENCES[1] M. Kebe, R Gadhafi, B. Mohammad, M. Sanduleanu, H. Saleh and M. Al-Qutayri, "Human vital signs detection methods and potential using radars: A review," Sensors, vol. 20, no. 5, p. 1454, 2020.[2] Y. Wu, H. Ni, C. Mao, J. Han and W. Xu, 'Non-intrusive human vital sign detection using mmWave sensing technologies: A review, "ACM Trans. Sen. Netw., vol. 20, n. 1, nov 2023.[3] M. Liebetruth, K. Kehe, D. Steinritz and S. Sammito, "Systematic literature review regarding heart rate and respiratory rate measurement by means of radar technology," Sensors, vol. 24, no. 3, 2024.[4] G. Patemiani, D. Sgrecda, A. Davoli, G. Guerzoni, P. Di Viesti, A. C. Valenti, M. Vitolo, G. M. Vitetta, and G. Boriani, "Radar-based monitoring of vital signs: A tutorial overview," Proceedings of the / EEC, vol. 111, no. 3, pp. 277-317, 2023.[5] B.-K. Park, O. Boric-Lubecke and V. M. Lubecke, "Arctangent demodulation with DC offset compensation in quadrature Doppler radar receiver systems," JEEE Transactions on Microwave Theory and Techniques, vol. 55, no. 5, pp. 1073-1079, 2007.[6] M. Nosrati and N. Tavassolian, "Accurate Doppler radar-based cardiopulmonary sensing using chest-wall acceleration," IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, vol. 3, n. 1, pp. 41-47, 2019.[7] C. Li and J. Lin, "Random body movement cancellation in Doppler radar vital sign detection," IEEE Transactions on Microwave Theory and Techniques, vol. 56, n. 12, pp. 3143-3152, 2008.[8] K. Han and S. Hong, "Differential-phase radar with amplitude-compensated complex signal demodulation for vital sign detection," in 2019 IEEE Asia-Pacific Microwave Conference (APMC), 2019, pp. 1390-1392.[9] C. Li, Y. Xiao and J. Lin, "Experiment and spectral analysis of a low-power Ka-band heartbeat detector measuring from four sides of a human body," IEEE Transactions on Microwave Theory and Techniques, vol. 54, n. 12, pp. 4464-4471, 2006.

[0010] I. Walterscheid and G. E. Smith, "Respiration and heartbeat monitoring using a distributed pulsed MIMO radar, " in 201739th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 2017, pp. 3449-3452.

[0011] Q. Wu, Z. Mei, Z. Lai, D. Li and D. Zhao, "A non-contact vital signs detection in a multi-channel 77GHz LFMCW radar system," IEEE Access, vol. 9, pp. 49614-49628, 2021.

[0012] M. Alizadeh, G. Shaker, J. C. M. D. Almeida, P. P. Morita and S. Safavi-Naeini, "Remote monitoring of human vital signs using mm-wave FMCW radar," IEEE Access, vol. 7, pp. 54958-54968, 2019.

[0013] A. Ahmad, J. C. Roh, D. Wang and A. Dubey, "Vital signs monitoring of multiple people using a FMCW millimeter-wave sensor," in 2018 IEEE Radar Conference (RadarConf 18), 2018, pp. 1450-1455.

[0014] L. Anitori, A. de Jong and F. Nennie, "FMCW radar for life-sign detection," in 2009 IEEE Radar Conference, 2009, pp. 1-6.

[0015] Q. Liu, H. Guo, J. Xu, H. Wang, A. Kageza, S. AIQami and S. Wu, "Non-contact non-invasive heart and respiration rates monitoring with MIMO radar sensing," in 2018 IEEE Global Communications Conference (GLOBECOM), 2018, pp. 1-6.

[0016] Q. Wu, X. Huang, Y. Chen, J. Li, and W. Zhu, "Multi-target respiration monitoring based on cumulative phase gradient approach," IEEE Geoscience andRemote Sensing Letters, pp. 1-1, 2023.

[0017] V. P. Tran, A. A. Al-Jumaily and S. M. S. Islam, "Doppler radar-based non-contact health monitoring for obstructive sleep apnea diagnosis: A comprehensive review, " Big Data and Cognitive Computing, vol. 3, n. 1, pag. 3, 2019.

[0018] P. Di Viesti, A. Davoli, G. Guerzoni and G. M. Vitetta, "Novel methods for approximate maximum likelihood estimation of multiple superimposed undamped tones and their application to radar systems," Techrxiv, 82021.

[0019] A. Albanese, L. Cheng, M. Ursino and N. W. Chbat, "An integrated mathematical model of the human cardiopulmonary system: Model development," American Journal of Physiology - Heart and Circulatory Physiology, vol. 310, n. 7, pp. H899-H921, 2016.

[0020] X. Yang, Y. Yu, H. Qian, X. Zhang and L. Zhang, "Body orientation and vital sign measurement with IR- UWB radar network," in 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), 2020, pp. 485-488.

[0021] D. Wang, S. Yoo and S. H. Cho, "Experimental comparison of IR-UWB radar and FMCW radar for vital signs," Sensors, vol. 20, n. 22, 2020.

[0022] M. Mercuri, Y.-H. Liu, I. Lorato, T. Torfs, A. Bourdoux and C. Van Hoof, "Frequency-tracking CW Doppler radar solving small-angle approximation and null point issues in non-contact vital signs monitoring," IEEE Transactions on Biomedical Circuits and Systems, vol. 11, no. 3, pp. 671-680, 2017.

[0023] M. Zhou, Y. Liu, S. Wu, C. Wang, Z. Chen and H. Li, "A novel scheme of high-precision heart rate detection with a mm-wave FMCW radar," IEEE Access, vol. 11, pp. 85 118-85 136, 2023.

[0024] S. M. Patole, M. Torlak, D. Wang and M. Ali, "Automotive radars: A review of signal processing techniques, " IEEE Signal Processing Magazine, vol. 34, n. 2, pp. 22-35, 2017.

[0025] A. Davoli, G. Guerzoni and G. M. Vitetta, "Machine learning and deep learning techniques for colocated MIMO radars: Atutorial overview," IEEE Access, vol. 9, pp. 33 704-33 755, 2021.

[0026] T. K. V. Dai, K. Oleksak, T. Kvelashvili, F. Foroughian, C. Bauder, P. Theilmann, A. E. Fathy and O. Kilic, "Enhancement of remote vital sign monitoring detection accuracy using multiple-input multiple-output 77 GHz FMCW radar," IEEE Journal of Electromagnetics, RF and Microwaves inMedicine and Biology, vol. 6, n. 1, pp.111-122, 2022.

[0027] M. Nosrati, S. Shahsavari, S. Lee, H. Wang and N. Tavassolian, "A concurrent dual-beam phased-array Doppler radar using MIMO beamforming techniques for short-range vital-signs monitoring," JEEE Transactions on Antennas and Propagation, vol. 67, no. 4, pp. 2390-2404, 2019.

[0028] I. Walterscheid, O. Biallawons and P. Berens, "Contactless respiration and heartbeat monitoring of multiple people using a 2-D imaging radar," in 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 2019, pp. 3720-3725.

[0029] Y. Rong, A. Dutta, A. Chiriyath, andD. W. Bliss, "Motion-tolerant non-contact heart-rate measurements from radar sensor fusion," Sensors, vol. 21, no. 5, 2021.

[0030] X. Zhang, X. Yang, Y. Ding, Y. Wang, J. Zhou and L. Zhang, "Contactless simultaneous breathing and heart rate detections in physical activity using IR-UWB radars," Sensors, vol. 21, n. 16, 2021.

[0031] Y. Rong, K. V. Mishra and D. W. Bliss, "Multiple moving targets heartbeat estimation and recovery using multi-frequency radars," in 2021 IEEE Radar Conference (RadarConf21), 2021, pp. 1-5.

[0032] P. Di Viesti, A. Davoli, G. Guerzoni and G. M. Vitetta, "Recursive algorithms for the estimation of multiple superimposed undamped tones and their application to radar systems," IEEE Transactions on Aerospace and Electronic Systems, vol. 59, n. 2, pp. 1834-1853, 2023.

[0033] J. Selva, "ML estimation and detection of multiple frequencies through periodogram estimate refinement," IEEE Signal Processing Letters, vol. 24, no. 3, pp. 249-253, 2017.

[0034] B. Zhang, B. Jiang, R Zheng, X. Zhang, J. Li and Q. Xu, "Pi-ViMo: Physiology-inspired robust vital sign monitoring using mmWave radars," ACM Trans. Internet Things, vol. 4, n. 2, maggio 2023.

[0035] B.-K. Park, S. Yamada, O. Boric-Lubecke and V. Lubecke, "Single-channel receiver limitations in Doppler radar measurements of periodic motion," in 2006 IEEE Radio and Wireless Symposium, 2006, pagg. 99-102.

[0036] S. Sarkka andL.Bayesian filtering and smoothing. Cambridge university press, 2023, vol. 17.

[0037] M. Pour Ebrahim, M. Sarvi and M. R Yuce, "Doppler radar system for sensing physiological parameters in walking and standing positions," Sensors, vol. 17, no. 3, p. 485, 2017.

[0038] Z. Chen, Y. Liu, C. Sui, M. Zhou and Y. Song, "A novel scheme for suppression of human motion effects in non-contact heart rate detection," IEEE Access, vol. 11, pp. 84241-84257, 2023.

[0039] M. Scherhaufl, F. Hammer, M. Pichler-Scheder, C. Kastl and A. Stelzer, "Radar distance measurement with Viterbi algorithm to resolve phase Ambiguity," IEEE Transactions onMicrowave Theory and Techniques, vol. 68, n. 9, pp. 3784-3793, 2020.

[0040] C. Li and J. Lin, "Complex signal demodulation and random body movement cancellation techniques for noncontact vital sign detection," m ' 2008 IEEEMTT-S InternationalMicrowave Symposiim Digest, 2008, pp. 567-570.

[0041] Infineon Technologies. (2024) Demo Position2Go. [Online], Available: https: / / www.infineon.com /

[0042] Shimmer Research. (2018) Shimmer3 ECG unit. [Online], Available: https: / / shimmersensing.com /

[0043] D. L. Colton and R Kress, Inverse acoustic and electromagnetic scattering theory. Springer, 1998, vol. 93.

[0044] A. Cuyt, V. Petersen, B. Verdonk, H. Waadeland and W. B. Jones, Handbook of continuedfractions for specialJunctions. Springer Science & Business Media, 2008.

[0045] G. N. Watson, A treatise on the theory of Bessel functions. Cambridge university press, 2012.

[0046] G. B. Arfken, H. J. Weber andF. E. Harris, Mathematical methods for physicists. Academic press, 2013.

[0047] G. H. Golub and C. F. Van Loan, Matrix computations. JHU press, 2013.APPENDIXA) Details rm interval compression and target selectionIn this appendix, model (21) is derived in the case of SIMO radar. To begin, the RHS of (14) is replaced with that of(l 8); this produces, for any virtual antenna (i.e., for any v E {0,1, ...,1^ — 1}), the two-dimensional complex sequence (X(v>[I, n]} wherewith I = 0,1, ... , No— 1 and n = 0,1, ... , NF— 1; here,represents the noise contribution. If the SNR is high, it is likely that the target detection strategy (20) selects, for a given n, the value of I (denoted lSv^ [n]) closest to ZV0 / nTs, such that

[0181] Note that the accuracy of the following approximation improves with increasing FFT order No. From (74), it easily follows thatwhereThis proves (21).B) Phase gradient analysis

[0182] In this appendix, the spectral content of the sequence {s^ [n] }, whose n-th element is defined by (21), is analyzed under the assumptions that: a) this sequence conveys displacement information from a single point target; and b) the target displacement is expressed by (3)-(5).For simplicity, in all following equations, the dependence on the block index p and antenna index v is omitted (e.g.,[ft] is denoted simply as s [n]) and noise contribution is neglected.

[0183] To start, substitute (4)-(5) into the RHS of (3) and then the resulting expression for target displacement into the RHS of (22) evaluated at t = nTF; this yields

[0184] Substituting the RHS of the last formula into that of (21) yieldswith n = 0,1, ... , 1VB— l; here,and b is defined in (15); note that parameter c is independent of the frame index n.

[0185] Based on the Jacobi-Anger expansion [43-46], (80) can be easily transformed into the formwhere s, [n] and sR[n] are defined respectively by (83) and (84),for any x andy.Figure 13 shows a representation of the amplitudes, versus <5b M, of the multiple sinusoidal components of the following sequences: (a) {s[n]} (82); (b) { / [ft]} (91); (c) {^[n]} (95). In all cases, parameters / c= 24 GHz, TF= 0.1 s, / b= 0.3 Hz, / h= 1 Hz, <5h M= 0.45 mm, and Jn = 1 are assumed.

[0186] From (82)-(85), it can be deduced that the amplitudes of the respiratory and cardiac harmonics, as well as those of the intermodulation components contained in the sequence {s [n]}, are proportional to c (81) and dependon the values of the coefficients {Jx(— 2n3b M / A)} and {Jy(47T<5hM / A)}.

[0187] Based on these mathematical results, assuming fc= 24 GHz (i.e., A = 12.5 mm), <5h M= 0.45 mm, and <5b MG [8,24] (see Section 2), it is possible to demonstrate that the FHs of the heart and respiration do not represent the dominant spectral components, since some of their HOHs have larger amplitude (see also

[0016] ).

[0188] This is exemplified in Figure 13a, where the amplitude of the multiple harmonics of heart and respiration is shown; note that coefficients C0,i and Ci,o we associated with the FHs of the heart and respiration, respectively.

[0189] These considerations lead to conclude that using the periodogram method to estimate BR and HR based on the amplitude spectrum of {s [n] } may lead to incorrect estimates of vital signs.

[0190] A recent work illustrated in

[0016] highlighted the advantages in estimating BR by replacing the complex sequence {s[n]; n = 0,1, ... , NB— 1} with the real sequence {g [n]; n = 0,1, ... , NB— An — 1], where#[n] A 3{ / [n]} (86) and f [n] = s[n] • (s[n + 4n])*, (87) represents the n-th element of the PG sequence (see (28)). If model (82) is adopted for s [n] , f [n] (87) can be easily transformed into the form / [n] = | c |2• yb[n] • yh[n], (88) whereandfor n = 0,1, ... , NB— An — 1; note that factors yb[n] and yb[n] represent the contributions due to respiration and heart, respectively. An alternative representation of f [n] (88) can be obtained by applying the Jacobi-Anger expansion to both yb[n] (89)andyb[n] (90); this yields f [n] (91), where / Jn] and / q[n] are defined respectively as (92) and (93):)for any x andy.

[0191] Equation (91) implies that the amplitudes of the respiratory and cardiac harmonics and those of the intermodulation components contained in the sequence { / [ft] } are proportional to the scaling factor | c |2(94) and to the values of the coefficients {Jx4n6b M / A. ■ sin(<nbJnTF / 2))} and {Jy2))}

[0192] Comparison between (94) and (85) highlights that the use of the PG operation (87) causes a reduction of the argument of the Bessel functions associated with respiration and heart by scaling factors sin / 2 andsin / 2, respectively. Generally, this leads to a reduction of the amplitude of HOHs and makes the FHsof the heart and respiration significantly stronger than the HOHs and intermodulation products.

[0193] This effect is exemplified in Figure 13b, where the amplitudes of the multiple harmonics of heart and respiration are shown for <5b MG [8,24] mm and <5h / M= 0.45 mm; note that the amplitude of the DC (due to the presence of term Do 0in (91)) is significant. Further advantages are provided by extracting the imaginary part of / [ft] , i.e., by evaluating g [n] (86). Indeed, substituting (91) into the RHS of (86) yields where / ] [ft] is defined as in (92).

[0194] This last result leads to the conclusion that the sequence {^ [ft]} contains no even-order harmonics, intermodulation products, or DC. This is exemplified in Figure 13 c, where the amplitudes of the multiple harmonics of heart and respiration are shown for <5b MG [8,24] mm and <5h / M= 0.4 mm.

[0195] These results demonstrate that: a) The FH of respiration is stronger than any other spectral component; b) The FH of the heart is stronger than any intermodulation product and any other HOH of heart and respiration.

[0196] These results hold for both sequences { / [ft] } and {g [ft]} if the presence of DC in the former sequence is neglected. Unfortunately, these mathematically obtained results do not fully agree with our experimental results(see the Section "Main insights derived from experimental results").

[0197] Indeed, our measurements showed that the spectrum of [g [n] } contains the second harmonic of the heart and that this component is often stronger than the FH. We believe the discrepancy between our mathematical results, and their experimental counterpart arises from the approximate model of chest displacement adopted in our analysis (see (3) and (5)).

Claims

CLAIMS1. Method for contactless monitoring of vital signs, in particular heart rate optionally in combination with breath rate, of at least one subject monitored by detecting and analyzing displacements in space of said at least one subject, said method comprising the steps of• Transmitting, by means of at least one transmitting element, one or more electromagnetic signals having known and predefined electromagnetic characteristics transmitted in proximity and at least partially towards said at least one subject;• Acquiring, by means of a number of receiving elements positioned in proximity of said at least one subj ect and such that there exist a plurality of virtual antennas, one or more electromagnetic signals reflected and at least partly coming from said subject by reflection of said electromagnetic signal;• Generating reflected digital signals by means of an analog-to-digital conversion of said reflected electromagnetic signals received by at least part of said receiving elements, said step comprising, for each virtual antenna, the steps of a. Performing a low-frequency and quadrature conversion of the acquired reflected analog signals; b. Generating digital samples representing the in-phase and quadrature components of the information content of said converted signals;• Performing a pre-processing phase of said reflected digital signals defined as a function of the known electromagnetic characteristics of said transmitted electromagnetic signal comprising the steps of a. Optionally, for each virtual antenna, combining the information contained in two or more of said reflected digital signals received at different and preferably consecutive time instants; b. Identifying in space said at least one subject, if not already identified; c. Collecting said combined reflected digital signals of the previous point a. by grouping them into data blocks acquired over a predefined time interval, where a data block comprises the information content from the reflected signals received in said predefined time interval.• Performing a processing phase of the pre-processed signals comprising the steps of a. For each virtual antenna and for each data block consequent to said virtual antenna, determine the phase gradient (PG) against the information relating to said pre-processed signals; b. For each virtual antenna and for each data block defined in the pre-processing phase, process the information content of said data block and the phase gradient as previously calculated, determining the dominant spectral components and describing said dominant spectral components in terms of significant parameters of such components, said parameters comprising, for each dominant spectral component, the frequency value and the related complex amplitude in terms of magnitude and phase;c. For each data block defined in the pre-processing phase, aggregate the parameters coming from said virtual antennas or from a selection thereof to form a parameter vector containing the information coming from said virtual elements or subset thereof; d. For each data block defined in the pre-processing phase, employ an iterative filter with memory on the parameter vector to generate a synthesized version of the parameter vector at least partially purified from the contribution of random movements of the body and of the multipath; e. For each filtered data block, select among one or more frequency values indicative of one or more vital signs of said at least one monitored subject.

2. Method according to claim 1, wherein the transmitted radar signal is of frequency modulated continuous wave (FMCW) type and said pre-processing phase comprises the steps of a. Optionally, for each virtual antenna, processing a combined signal, combining the signals received from said virtual antenna and obtained by reflection of the signal corresponding to the different chirps of a chirp train; b. Determining the frequency components of said combined signal, preferably by means of discrete Fourier transform (DFT) optionally preceded by zero-padding operation on the digital samples of said combined signal; c. Determining the position of the target by identifying the maximum of said frequency components; d. For each virtual antenna, grouping in blocks of NBdata said frequency components associated to said position of the target.

3. Method according to claim 1, wherein the transmitted radar signal is of continuous wave (CW Doppler) type, the pre-processing phase of said reflected digital signals comprising the steps of a. Optionally, for each virtual antenna, processing a combined signal, combining the reflected signals received from said virtual antenna and obtained by reflection of the signal corresponding to a sinusoid; b. For each virtual antenna, grouping in blocks of NBdata said samples of said combined signal.

4. Method according to claim 1, wherein the transmitted radar signal is of stepped-frequency continuous wave (SFCW) type, the pre-processing phase of said reflected digital signals comprising the steps of a. Optionally, for each virtual antenna, processing a combined signal, combining the reflected signals received from said virtual antenna and obtained by reflection of the signal corresponding to the different periods of a stepped waveform; b. Determining the time-domain signal of said combined signal, preferably by means of inverse discrete Fourier transform (IDFT) optionally preceded by zero-padding operation on the digital samples of said combined signal; c. Determining the position of the target by identifying the maximum of said time-domain signal;d. For each virtual element, grouping in blocks of NBdata said samples of the time-domain signal associated to said position of the target.

5. Method according to claim 1, wherein the transmitted radar signal is of impulse radio ultra-wideband (IR- UWB) type, the pre-processing phase of said reflected digital signals comprising the steps of a. Optionally, for each virtual antenna, processing a combined signal, combining the reflected signals received from said virtual antenna and obtained by reflection of the signal corresponding to the different pulses of a pulse train; b. Determining the position of the target by identifying the maximum of said reflected signal received by said virtual antenna and obtained by reflection of the signal corresponding to the different pulses of a pulse train; c. For each virtual antenna, grouping in blocks of lVBsaid samples of the time-domain signal associated to said position of the target.

6. Method according to one or more of the preceding claims wherein the virtual antennas of said plurality of virtual antennas are mutually time-coordinated and the acquisition of the reflected signals from the receiving elements occurs sequentially according to a mode providing the acquisition of a single element at a time and according to a sequence known a priori, such as, in round robin mode or according to a pseudo-random order.

7. Method according to one or more of the preceding claims wherein in said pre-processing and processing phases it is assumed that the movements of said subject are two or three orders of magnitude lower than the distance of said at least one subject from said transmitting and receiving elements.

8. Method according to one or more of the preceding claims wherein at least a second transmitting element is provided and the steps of the method are reiterated for each of the transmitting elements.

9. Method according to one or more of the preceding claims wherein in the pre-processing and / or processing phase the further steps are provided of- Predefining a range interval of distances between said receiving elements and the monitored subject, said range interval being between a minimum distance value and a maximum distance value;- Determining the approximate distance of objects in the environment as a function of the characteristics of said reflected digital signals;- Excluding the information relating to obj ects whose determined distance does not fall within said range interval.

10. Method according to one or more of the preceding claims wherein said data blocks are grouped such that the temporal duration of the information of a block (1VBTSF) is chosen so as to be greater than the duration of two consecutive breathing periods (2Th rnin).

11. Method according to one or more of the preceding claims wherein the iterative filtering with memory occurs through a filtering operation based on a Bayesian statistical approach, wherein:- the State Transition Model, i.e. the temporal evolution of the system state, is represented by the evolution of heart rate (HR), and where present breath rate (BR), as state variables;- the Measurement Model, i.e. the description of how observable measurements relate to the system state, is obtained from the measurements derived from radar observations;- the Time Update or Prediction, i.e. the evolution of the state probability distribution based on the state model, is represented by updating the a priori probabilities of heart rate and breath rate using the previous state distribution and the state transition model;- the Measurement Update or Correction, i.e. the update of the state probability distribution based on new observations, occurs by updating the a posteriori probabilities of heart rate and breath rate based on new radar measurements.

12. Method according to one or more of the preceding claims 1 to 11 wherein the vital sign of said at least one subject is heart rate and comprising the steps of- Setting a heart rate frequency range of interest, comprising a minimum heart rate frequency value ( / h,min) and a maximum heart rate frequency value ( / h,max)i- Limiting said determination of the dominant spectral components to a search range going from said minimum heart rate frequency value ( / h,min) to twice said maximum heart rate frequency value (2 / h rnax- Determining the heart rate in real time among the values within said range by discriminating between the first harmonic (FH) of the detected signal if within the range [ / h.min > / h,max]ar|d the second harmonic (SH) of the detected signal if within the range [2 / h min, 2 / h max] ..

13. Method according to one or more of the preceding claims 1 to 12 wherein the vital sign of said at least one subject comprises breath rate and comprising the steps of- Setting a breath rate frequency range of interest, comprised within a minimum breath rate frequency value ( / b,min)ar|d a maximum breath rate frequency value ( / tymax )i- Limiting said determination of the dominant spectral components to said breath rate frequency;- Determining the breath rate in real time among the values within said breath rate frequency.

14. Method according to one or more of the preceding claims wherein said step of selecting the data block related to at least one of said virtual antennas comprises the sub-steps of- Determining a quality index related to the signal associated to each virtual antenna, said quality index preferably calculated as a function of the intensity of the peaks of the dominant spectral components related to the vital signs of said at least one subject;- Selecting a subset of virtual antennas associated to the best quality index.

15. Method according to the previous claim, wherein the size of said subset of virtual antennas is predetermined as constant or is variable as a function of a predefined quality threshold parameter, said subset being composed of virtual antennas whose quality index value exceeds said threshold parameter.

16. Method according to one or more of the preceding claims wherein the breath rate determination is present and the processing phase further comprises the determination of the apnea event, calculated based on the disappearance of the breathing peak in the spectral region of interest, and optionally the breath rate determination occurs only if no apnea event is ongoing.

17. System for contactless monitoring of the vital signs of at least one monitored subject by detecting and analyzing displacements in space of said at least one subject, said system comprising: at least one transmitting element for transmitting radar signals; at least two receiving elements for receiving reflected radar signals; at least one control unit, said control unit being connected or connectable in data exchange with said transmitting and receiving elements and comprising: o a source of a reference timing signal for synchronization of said receiving and transmitting elements; o a memory for loading software encoding the steps of the method of one or more of the preceding method claims; o a processing unit capable of executing the steps of said software.

18. System according to claim 17 comprising one or more Single-Input Multiple-Output (SIMO) radar devices, distributed in space according to a known and predefined spatial distribution.

19. Computer program, such as firmware or software, comprising instructions in code form that allow the processing unit of the system according to claims 17 or 18, executing such program, to perform the steps of one or more of the preceding method claims.

Images