Modelling, Analysis, and Simulation of the Micro-Doppler Effect in Wideband Indoor Channels with Confirmation Through Pendulum Experiments

This paper is about designing a 3D no n-stationary wideband indoor channel model for radio-frequency sensing. The proposed channel model allows for simulating the time-variant (TV) characteristics of the received signal of indoor channel in the presence of a moving object. The moving object is modelled by a point scatterer which travels along a trajectory. The trajectory is described by the object’s TV speed, TV horizontal angle of motion, and TV vertical angle of motion. An expression of the TV Doppler frequency caused by the moving scatterer is derived. Furthermore, an expression of the TV complex channel transfer function (CTF) of the received signal is provided, which accounts for the influence of a moving object and fixed objects, such as walls, ceiling, and furniture. An approximate analytical solution of the spectrogram of the CTF is derived. The proposed channel model is confirmed by measurements obtained from a pendulum experiment. In the pendulum experiment, the trajectory of the pendulum has been measured by using an inertial-measurement unit (IMU) and simultaneously collecting CSI data. For validation, we have compared the spectrogram of the proposed channel model fed with IMU data with the spectrogram characteristics of the measured CSI data. The proposed channel model paves the way towards designing simulation-based activity recognition systems.


Introduction
In wireless communications, the compound Doppler effect caused by the moving objects or bodies opened up opportunities for many applications. These applications track the scattered wave components by the moving bodies for drone detection [1], gesture recognition [2], human gait assessment for diagnosis and rehabilitation [3], and tracking human activities using no n-wearable radio-frequency-based (RF-based) elder-care [4]. Such waves contain the micro-Doppler effects corresponding to the moving bodies.
In channel modelling, the Doppler effect caused by moving scatterers has been modelled in two-dimensional (2D) stationary fixed-to-mobile radio in [5]. Then, this model has been extended for 2D no n-stationary fixed-to-fixed (F2F) indoor channels by considering the time-variant (TV) speed of the moving scatterer, angle of motion, angle of arrival, and angle of departure [6]. Later on, the TV Doppler frequency caused by the moving scatterer has been incorporated in three-dimensional (3D) channels by taking into account the TV azimuth angles of motion (AAOM), elevation angle of motion (EAOM), azimuth angle of departure (AAOD), elevation angle of departure (EAOD), azimuth angle integrated to obtain TV velocities and displacements (trajectories). The TV velocities and displacements suffer from linear and quadratic drifts, respectively. In [45,46], these drift problems were addressed by employing zero update algorithms.
In this paper, we design a 3D no n-stationary wideband channel model for activity recognition. We model a'moving object as a moving scatterer. The expressions for the TV speed, AAOM, EAOM, EAOD, AAOD, EAOA, and AAOA corresponding to the moving scatterer are presented. Then, the expression for the TV Doppler frequency caused by the moving scatterer is provided. Furthermore, the instantaneous channel phase of the moving scatterer and the complex channel transfer function (CTF) of the no n-stationary F2F channels are determined. Next, an approximate solution of the spectrogram of the complex CTF is presented to provide insight into the TV Doppler power characteristics of this model. We perform measurements of the calibrated CSI and IMU, simultaneously during a moving object experiment. The CSI data is calibrated by using a B2B connection to eliminate the TV phase distortions. Then, we feed the proposed channel model with the measured trajectory using the IMU. Finally, we show that the spectrogram of the calibrated measured CSI data and the channel model are matching. The TV mean Doppler shifts computed from both spectrograms are matching as well. The proposed channel model paves the way towards the design of software-RF-based human activity recognition and fall detection systems.
The rest of the paper is structured as follows. Section 2 exhibits the 3D multipath propagation scenario with moving and fixed objects. Section 3 presents the expressions of the TV channel parameters and the complex channel transfer function. The approximate analytical solution of the spectrogram of the complex channel function is provided in Section 4. Section 5 discusses the numerical and measurement results. The conclusion and future work are discussed in Section 6.

The 3D Geometrical Model
We consider the 3D geometrical model of a 3D multipath propagation channel shown in Figure 1. This figure shows a fixed Wi-Fi transmitter T x and a fixed Wi-Fi receiver R x which operate according to the IEEE 802.11n standard [19] with carrier frequency f 0 = 5.32 GHz and bandwidth B = 20 MHz. The positions of T x and R x are denoted by x T , y T , z T and x R , y R , z R , respectively. A moving object whose center of mass (CoM) is modelled for simplicity by a single moving (point) scatterer S M initially located at (x M , y M , z M ). The trajectory of the moving scatterer S M is described by a TV velocity vector v M (t) which can be expressed by the TV speed v M (t), the TV AAOM α v M (t), and the TV EAOM β v M (t). Each fixed object is modelled by a fixed scatterer S F m ( ) for m = 1, 2, . . . , M, where M denotes the number of fixed scatterers (objects). The TV parameters β T M (t), α T M (t), β R M (t), and α R M (t) designate the TV EAOD, TV AAOD, TV EAOA, and TV AAOA, respectively. We assume single-bounce scattering, i.e., each wave that is launched from T x is bounced by either a fixed scatterer S F m or a moving scatterer S M before arriving at R x .

The Channel Transfer Function
The TV velocity vector v M (t) of the moving scatterer S M is presented as where [·] T denotes the vector transpose operation. The velocities v M,x (t), v M,y (t), and v M,z (t) can be expressed in terms of the TV speed v M (t), where The function atan2(·) in (5) represents the four-quadrant inverse trigonometric tangent function that provides an azimuth angle α v M (t) ranging from −π to π, unlike the regular arctan(·) function that provides an angle ranging from −π/2 to π/2. Note that the elevation angle β v M (t) is within the range from −π/2 to π/2. By using the components of the TV velocity vector v M (t) in (2)-(4), one can compute the displacements x M (t), y M (t), and z M (t) of the moving scatterer S M as The superscript q in (17) represents the subcarrier index of OFDM communication systems that follows the IEEE 802.11n standard [19]. The parameter ∆ f (q) in (17) designates the subcarrier frequency which is given by for q ∈ {−28, −26, . . . , −2, −1, 1, 3, . . . , 27, 28}. In (20), the parameter ∆ represents the subcarrier frequencies difference which has a value of 312.5 kHz [19]. The function H M (t, ∆ f (q) ) designates the complex CTF of the moving scatterer S M and the parameter H F,m denotes the complex CTF corresponding the mth fixed scatterer S F m . The expression in (17) is similar to the one in [47] [Equation (21)]. The only difference is that the multipath effect associated with the fixed scatterers is taken into account by adding the second term in (17). The first term in (17) designates the TV part of the CTF corresponding to the moving scatterer S M with a fixed path gain c M and stochastic phase process θ M − 2π( f 0 + ∆ f (q) )τ M (t) associated with the qth subcarrier [see (18)]. The second term in (17) is time-invariant and represents the sum of the M received multipath components corresponding to the M fixed scatterers. Each component of the second term in (17) is characterized by a fixed path gain c m,F and a random phase variable θ m,F due to the interaction with the mth fixed scatterer S F m [see (19)]. It should be mentioned that the phases θ M and θ m,F are identically and independently distributed (i.i.d), each follows a uniform distribution over −π and π, i.e., θ M , θ m,F ∼ U (−π, π ]. The model presented in (17) which can be found in [47] [Equation (22)] as [48] f (q) where function f (q) max (t) designates the maximum Doppler shift caused by the moving scatterer S M which is given by From the expression in (21) and the relationship f , one can conclude that if the moving scatterer S M moves away from the T x and R x vicinity, the TV propagation delay τ M (t) and its slopeτ M (t) increase and the Doppler effect f M (t) has negative values, and vice versa. To obtain an approximate solution for the spectrogram of the CTF H(t, ∆ f (q) ) that will be discussed in Section 4, the Doppler frequency f for t l < t ≤ t l+1 and l = 0, 1, . . . , M,l denotes the slope of the approximated Doppler frequency f (q) The difference between two consecutive time instances t l+1 and t l , i.e., δ = t l+1 − t l is the same for all values of l = 0, 1, . . . , L − 1.
The TV mean Doppler shift B (1) f (q) (t) of the proposed 3D channel model can be computed by using (21) as [49] B (1) The expression in (25) denotes the squared path gain c 2 M multiplied by the Doppler frequency caused by the moving scatterer f (1) f (q) (t) will have values closer to those of the Doppler frequency of the moving scatterer f

Spectrogram Analysis
In this paper, we employ the spectrogram approach [50] to reveal the TV Doppler power spectrum of the proposed channel model. The spectrogram S H (q) ( f , t) of the CTF H(t, ∆ f (q) ) corresponding to the qth subcarrier index is computed in three steps. First, a sliding window w(t) is multiplied by the CTF H(t, ∆ f (q) ). In this paper, we choose a Gaussian window function [50] [Equation (2.3.1)] where the parameter σ w is called the Gaussian window spread. The window function w(t) is real, positive, and even. It has a no rmalized energy, i.e., ∞ −∞ w 2 (t) = 1. By multiplying the window function w(t) by the CTF where the variables t and t are the local time and the running time, respectively. The second step is to compute the short-time Fourier transform (STFT) By using the approximation of the TV Doppler shift provided in (23), the STFT X H (q) ( f , t) associated with the qth subcarrier is obtained as The expression in (29) is a complex Gaussian function with a TV mean f M,l (t) and a complex variance σ 2 x,M,l . Note that the complex variance σ 2 x,M,l in (30) is dependent on the slope k (23)]. The last step is to obtain the spectrogram S H (q) ( f , t) associated with the qth subcarrier by squaring the magnitude of the STFT X H (q) ( f , t);, i.e., are called the auto-term and the cross-term of the spectrogram S H (q) ( f , t), respectively. The auto-term is given by for t l < t ≤ t l+1 , where The auto-term S (a) (33) is an approximation that provides insight into the Doppler power spectrum of the proposed 3D no n-stationary channel model presented in Section 2. This term is real, positive, and consists of a sum of M + 1 weighted Gaussian functions. The first Gaussian function, which is due to the moving scatterer S M is weighted by the squared path gain c 2 M and centered on the approximated TV Doppler frequency f The cross-term S (c) (36) represents the undesired spectral interference term consisting of M(M + 1)/2 components which reduce the resolution of the spectrogram. This term is real but not necessarily positive. The operators {·} and (·) * denote the real value operator and the complex conjugate operator, respectively. The cross-term in (36) consists of two terms. The first term of (36) designates the sum of the components corresponding to the spectral interference caused by the fixed scatterers. The mth component of the second term denotes the spectral interference between the moving scatterer S M and the mth fixed scatterer S F m . The cross-term S (c) (36) is dependent on the random phases θ M and θ F,m unlike the auto-term S (a) t). Hence, the cross-term S (c) can be eliminated by taking the average over the random phases, i.e., E{S (c) The TV mean Doppler shift can be obtained by using the spectrogram as follows If the spectrogram S H (q) ( f , t) in (37) is replaced by the auto-term S (a)

Measurements and Numerical Results
In this section, we discuss and compare the TV Doppler power characteristics of our proposed channel model with those of measured CSI data. We will describe the processing of the measured trajectory during the measurements.

Measurement Scenario
To complement the TV Doppler power characteristics of the proposed channel model, measurements have been performed. The CSI data and the trajectory of a pendulum (moving object) have been measured simultaneously. Two laptops have been used for measuring the CSI as Wi-Fi T x and R x . An IMU sensor fusion has been used to measure the trajectory of the pendulum. Figure 2 illustrate the measurement scenario in xy and xz planes, respectively. The pendulum was a 3 kg medicine ball, covered with aluminum foil and attached to the ceiling by a rope, and was swinging in a horizontal direction perpendicular to the line-of-sight (LoS). The distance between the ceiling and the center of mass (CoM) of the ball L was 1.17 m and the horizontal distance between Wi-Fi T x antenna and the CoM of the ball was 1.5 m. The distance between Wi-Fi T x and R x antennas was 2 m and they had the same height value of 1.18 m. The initial location of the moving scatterer (ball) was the origin. The pendulum displacements x M (t) and z M (t) are computed as follows [36]: where g denotes the acceleration of gravity. The parameters x max and L in (38)

Motion Capturing Using IMU
A MetaMotionR sensor fusion (IMU) [51] was attached to the swinging ball. A smartphone was connected via Bluetooth to control the IMU and log the data files. The IMU was used to record quaternions and linear accelerations during the experiment. Euler angles were computed by using the recorded quaternions to rotate the measured linear accelerations. Next, the raw rotated linear accelerations were smoothed by using quadratic regression provided by the signal analysis toolbox in MATLAB 2019a. After that, the rotated linear accelerations were integrated and double integrated to obtain the velocities and the displacements (trajectories), respectively. Due to measurement errors of the IMU, the velocities and the displacements suffer from linear and quadratic drifts, respectively. To overcome this drift issue, zero-update (ZUPT) algorithms [45] were employed. Since the pendulum motion is periodic, its horizontal and vertical velocities reach zero when the horizontal and vertical accelerations approach their maximum or minimum values. Similarly, the values of horizontal and vertical displacements approach zero values when the velocities tend to their maximum or minimum values. Hence, by searching for the indices corresponding to the local maximum or minimum values of the accelerations, the velocity drift is removed between two consecutive indices. Also, by knowing the indices of the local maximum or minimum values of the drift-eliminated velocities, the drift of the displacement is removed. The source code of the ZUPT algorithm, where the sensors are placed on the toes of a walking person for position tracking, is available online [46]. This algorithm was repeated to also eliminate the drift of the displacement. Figure 3a depicts the TV drift-free horizontal displacements x M (t) of the captured data from the IMU and the mechanical model of the pendulum in (38) by using the pendulum parameters shown in Figure 2. The TV drift-free vertical displacements z M (t) of the captured data from the IMU and the mechanical model of the pendulum in (40) by using the pendulum parameters shown in Figure 2, are depicted in Figure 3b. A minimal error is noticed between the IMU data and the model in the order of centimeters during the whole interval of 15 s.

Capturing CSI Data
The CSI tool in [18] was installed to capture the CSI data (RF signals). Two HP Elitebook 6930p laptops equipped with Intel NIC5300 were used. An Ubuntu 14.04 LTS operating system was installed on both laptops. One laptop was the transmitter station in injector mode and the other laptop was the receiver operating in monitor mode. The carrier frequency f 0 was set to 5.32 GHz corresponding to channel 64 according to IEEE 802.11n standards [19]. The sampling frequency and the bandwidth were set to 1 KHz and 20 MHz, respectively. TV phase distortions exist due to carrier frequency offset [21][22][23], sampling frequency offset [24][25][26], and packet boundary delay [27,28]. These TV phase distortions were eliminated by using a B2B connection between the transmitter station and the receiver station as described in [35]. Since there was only one RF transmission port in the Wi-Fi T x , an RF power splitter ZFSC-2-10G+ from with two output ports was used. One of the output ports was used for the B2B connection and the other one was connected to the transmitting antenna. At the Wi-Fi receiver laptop, one of the ports was used for the B2B connection, and another port was connected to the receiver antenna. The port used for the B2B connection was connected to a 30 dB attenuator. RF cables 141-1MSM+ from Mini-Circuits ® were used as well. The processing of the captured CSI data was done by using MATLAB R2019a. Two matrices are stored in a file. One matrix contains the CSI data that corresponds to the captured signal with the fingerprint information associated with the motion of the pendulum and TV phase distortions. The other matrix corresponds to the B2B connection, i.e., it only contains the TV phase distortions. Then, the matrix that contains the fingerprint information and TV phase distortions is divided by the matrix corresponding to the B2B connection in elementwise form. The output matrix resulting from the elementwise division only contains the fingerprint information. After the elementwise division, a highpass filter has been used to reduce the power of zero-frequency components associated with the fixed scatterers and/or the line-of-sight. Regarding the channel model and its spectrogram, Figure 4a shows the block diagram of the proposed channel model discussed in Section 3 fed with IMU data as inputs and the computation of the spectrogram. Figure 4b shows the block diagram of the proposed channel model discussed in Section 3 fed with the mechanical model as inputs and the computation of the spectrogram. Note that the difference between Figure 4a and Figure 4b is how the trajectories are obtained to feed the channel model. If they are measured using IMU, then the preprocessing mentioned Section 5.2 should be considered before feeding them to the simulator. If they are computed using the expressions in (38)-(40), then they can be fed into the simulator directly. The channel model can be fed with the TV displacements from either the IMU (after applying ZUPT) or the mechanical model presented earlier in Section 5.2, as inputs. The carrier frequency of the simulator f 0 was set to 5.32 GHz for consistency with CSI measurement scenario. The number of the fixed scatterers M was chosen to be 6. The initial location of the moving scatterer S M and the locations of the Wi-Fi T x , Wi-Fi R x were set according to the experiment scenario as presented in Figure 2, i.e., they can be located anywhere, but the distances should be the same as those illustrated in Figure 2. Then, the TV displacements as presented in Figure 3 were added to the initial location of the moving scatterer S M .
After that, the TV Doppler frequency f (q) M (t) caused by the moving scatterer S M was computed according to (21). The path gains of the moving scatterer S M and each fixed scatterer S F m were computed by respectively. The parameter η is used to balance the contribution of the fixed and moving scatterers and was set to 0.8. The phases θ M and θ F,m were generated as realizations of random variables with uniform distribution from −π to π. Next, the STFT X H (q) ( f , t) for each subcarrier index q was computed according to (28). The window spread parameter σ w was set 31.1 ms. Finally, the spectrogram S( f , t) (orS( f , t) in case of using IMU data as inputs) was computed as the squared magnitude of the sum of the STFT over the subcarriers by the following expression For computing the spectrogram of the recorded CSI as exhibited in Figure 4c, the CTFĤ (q) (t, ∆ f (q) ) is recorded. Then, the STFTX H (q) ( f , t) was computed for each subcarrier q. After that the spectrogram S( f , t) is computed according to (42). Figure 5a-c exhibit the spectrograms ofS( f , t), S( f , t), andŜ( f , t) of the channel model with IMU data as inputs, the channel model fed with the mechanical model as inputs, and the recorded CSI data, respectively. It is shown that the TV Doppler power characteristics depicted by the spectrograms S( f , t), S( f , t), andŜ( f , t) in Figure 5a-c are fairly similar to each other, respectively. In Figure 5a-c, the Doppler frequency associated with the moving scatterer (pendulum) S M has negative values when the pendulum swings away from the Wi-Fi T x and Wi-Fi R x antennas and has positive values when it swings towards them. The Doppler frequency corresponding to the moving scatterer (pendulum) S M approaches zero values at the time instants in which the moving scatterer reach its local maximum and minimum displacement values [see Figure 3a,b]. Therefore, the speed of the pendulum v M (t) approaches zero values. Thus, the Doppler shift at these instants is zero according to (22). Figure 6 depicts the TV mean Doppler shiftsB (1) (t), B (1) (t), andB (1) (t) computed from the spectrogramsS( f , t), S( f , t), andŜ( f , t) using (37), respectively. There is a good match betweeñ B (1) (t), B (1) (t), andB (1) (t). The mean Doppler shifts have negative values at the time instants in which the pendulum (moving scatterer S M ) swings away from the Wi-Fi T x and Wi-Fi R x antennas.
They have positive values when the pendulum swings towards the Wi-Fi T x and Wi-Fi R x antennas. The TV mean Doppler shiftsB (1) (t), B (1) (t), andB (1) (t) approach zero values at the moments when the pendulum reaches its local maximum and minimum displacement values. There exists a slight drift in the values of the mean Doppler shiftB (1) (t) in between the time instants t ≈ 11.5 s and t ≈ 12.7 s due to the no ise of the measured CSI signal.   For quantitative evaluation, we collected CSI and IMU data for 20 experiments, i.e., K = 20. From the collected data measurement, we computed the normalized-mean-square-error (NMSE) γ k between the TV mean Doppler shiftB (1) k (t) of the proposed channel model fed with the IMU data as inputs and the TV mean Doppler shiftB (1) k (t) of the CSI data according to

Channel simulator
for k = 1, 2, . . . , K, where the parameter T obs denotes the observation interval which was set to 15 s, i.e., T obs = 15 s. Figure 7 depicts the NMSE γ k for each experiment. The maximum NMSE belongs to the first experiment and has a value of 0.1829, whereas the minimum NMSE, with a value of 0.0477, belongs to fourteenth experiment. The average NMSE equals 0.0932, and the variance of the NMSE is 0.0013.

Conclusions
In this paper, we proposed a no n-stationary wideband channel model and its TV Doppler power characteristics when there is a moving object in the 3D space. We derived the TV Doppler shift caused by the moving object in terms of the TV speed, AAOM, EAOM, AAOD, EAOD, AAOA, and EAOA. The TV Doppler characteristics of the proposed channel model were analyzed by using the spectrogram. Furthermore, we provided the approximate solution of the spectrogram of the channel model. We validated the proposed channel model by measuring the trajectory of the moving object using an IMU and calibrated CSI with B2B connection, simultaneously. Then, we fed the channel model with the trajectory data extracted form the IMU. The results showed a good agreement between the measured CSI and the channel model in terms of the spectrogram and the mean Doppler shift. We conclude that the proposed channel model can be used for designing simulation-based HAR systems. For the future, we aim to extend the proposed channel model for human activity recognition by modelling the moving human as multiple moving scatterers.