Hands-on Wireless Sensing with Wi-Fi: A Tutorial

Abstract.

With the rapid development of wireless communication technology, wireless access points (AP) and internet of things (IoT) devices have been widely deployed in our surroundings. Various types of wireless signals (e.g., Wi-Fi, LoRa, LTE) are filling out our living and working spaces. Previous researches reveal the fact that radio waves are modulated by the spatial structure during the propagation process (e.g., reflection, diffraction, and scattering) and superimposed on the receiver. This observation allows us to reconstruct the surrounding environment based on received wireless signals, called “wireless sensing”. Wireless sensing is an emerging technology that enables a wide range of applications, such as gesture recognition for human-computer interaction, vital signs monitoring for health care, and intrusion detection for security management. Compared with other sensing paradigms, such as vision-based and IMU-based sensing, wireless sensing solutions have unique advantages such as high coverage, pervasiveness, low cost, and robustness under adverse light and texture scenarios. Besides, wireless sensing solutions are generally lightweight in terms of both computation overhead and device size. This tutorial takes Wi-Fi sensing as an example. It introduces both the theoretical principles and the code implementation of data collection, signal processing, features extraction, and model design. In addition, this tutorial highlights state-of-the-art deep learning models (e.g., CNN, RNN, and adversarial learning models) and their applications in wireless sensing systems. We hope this tutorial will help people in other research fields to break into wireless sensing research and learn more about its theories, designs, and implementation skills, promoting prosperity in the wireless sensing research field.

1. Wireless Sensing Background

1.1. What is Wireless Sensing

Various sensors and sensor networks have thoroughly extended human perception of the physical world. Nowadays, numerous sensors have been deployed to complete various sensing tasks, resulting in a significant deployment and maintenance overhead. This problem becomes increasingly troublesome when a large sensing scale is in demand. Taking indoor person tracking as an example, a specialized tracking system only covers a room-level area, which is too small compared with the moving region during a person’s daily life. Multiple tracking systems are needed to achieve practical sensing coverage in realistic living environments, e.g., houses, campuses, markets, airports, and offices, and the cost inevitably ramps up.

Given the cost limitations of the sensors, many pioneers tried to figure out an alternative solution during the past decade. Nowadays, various types of signals (e.g., Wi-Fi, LoRa, LTE) are filling out our living and working spaces for wireless communication, which can be leveraged to capture the environmental changes without causing extra overhead. According to the electromagnetism theory, the radio signals emitted by the transmitter (Tx) experience various physical phenomena such as reflection, diffraction, and scattering during the propagation process and form into multiple propagation paths. In this way, the superimposed multipath signals collected by the receiver carry spatial information about the signal propagation environment. Relying only on the ambient wireless signals and ubiquitous communication devices, wireless sensing emerges as a novel paradigm for environment sensing.

In recent years, wireless sensing technology has attracted many research interests to bring wireless sensing from the imagination into reality, by boosting sensing granularity, improving system robustness, and exploring application scenarios. Many of their works on wireless sensing have been published in flagship conferences and journals, such as ACM SIGCOMM, ACM MobiCom, ACM MobiSys, IEEE INFOCOM, USENIX NSDI, IEEE/ACM ToN, IEEE JSAC, and IEEE TMC. In addition, many famous companies are also exploring the productization of sensorless sensing, launching various IoT devices for human-computer interaction, security monitoring, and health care.

1.2. Comparison of Wireless Sensing and Computer Vision

Typical RF signals (300 kHz - 300 GHz) and visible light signals (380 THz - 750 THz) are essentially electromagnetic (EM) waves. When propagating in our physical world, the EM waves experience a variety of physical phenomena such as reflection, diffraction, and scattering. Multipath signals are eventually superimposed and received by the receiver. Therefore, the received superimposed signals carry the physical information of the signal propagation space.

Both the RF-based and the vision-based sensing algorithms share similar processes. They first analyze the received signals (radio signal at the antenna or visible light at the camera lens), from which the features reflecting the propagation space are extracted and finally resolved by algorithms to realize the sensing of the surrounding environment.

Compared with vision-based sensing, wireless sensing solutions have unique advantages such as high coverage (Chi et al., 2021), pervasiveness, low cost, and robustness under adverse light and texture scenarios (Chi et al., 2022).

Figure 1. Comparision of vision-based sensing and RF-based sensing processes.

1.3. Wireless Sensing Applications

Wireless sensing systems are capable of perceiving changes in surrounding environments, objects, and human bodies. In this subsection, we take passive human sensing applications as an example, which refers to a human-centered sensing application that doesn’t require the user to carry any device. Therefore, such an sensing application is also termed as device-free sensing or non-invasive sensing. Passive human sensing enables a wide range of applications, including smart homes, security surveillance, and health care.

In smart home applications, passive human sensing recognizes a person’s behavior or intention based on the user’s physical locations, gestures, and postures. Passive human sensing brings a better user experience without imposing restrictions on the user. For example, users can remotely control electrical devices, e.g., television, computer, or washer, by merely performing gestures in the air (Zheng et al., 2019; Abdelnasser et al., 2015). Likewise, when playing video games, users can interact with the computer by performing different postures (Jiang et al., 2020).

In security surveillance applications, traditional methods adopt infrared or RGB cameras to monitor illegal invasions, protect valuable properties, and deal with emergencies. However, cameras are constrained by the limited field of view or blockage of opaque or metallic objects, rendering these methods to fail when the target is not in the Line-of-Sight (LoS) area of the surveillance camera or hidden behind other objects. In contrast to visual surveillance, wireless sensing technology leverages radio signals, which provide omnidirectional coverage around the wireless devices and are less prone to blockages. For example, wireless signals can be used to detect illegal intrusions (Wu et al., 2015; Qian et al., 2014). Besides, they can also be used to detect if properties have been moved from their original places.

In health care applications, passive human sensing can be leveraged to detect vital signals such as human respiration, heartbeat, gait, and accidental fall. Specifically, some researchers have exploited Wi-Fi signals to detect human respiration (Wang et al., 2016) for sleep monitoring. Some other works (Zhang et al., 2020; Wu et al., 2020) have extracted gait patterns from Wi-Fi signals to recognize human identity. Recently, Wi-Fi signals have been further used to detect accidental falls to relieve the need for wearable sensors (Hu et al., 2020; Palipana et al., 2019).

2. Understanding CSI

Channel state information (CSI) lays the foundation of most wireless sensing techniques, including Wi-Fi sensing, LTE sensing, and so on. CSI provides physical channel measurements in subcarrier-level granularity, and it can be easily accessed from the commodity Wi-Fi network interface controller (NIC).

CSI describes the propagation process of the wireless signal and therefore contains geometric information of the propagation space. Thus, understanding the mapping relationship between CSI and spatial geometric parameters lays the foundation for feature extraction and sensing algorithm design.

This section focuses on two mainstream CSI models: the ray-tracing model and the scattering model. The two models are based on two perspectives of understanding the signal propagation process. Thus, they have unique advantages and apply to different scenarios.

2.1. Ray-tracing Model

In typical indoor environments, a signal sent by the transmitter arrives at the receiver via multiple paths due to the reflection of the radio wave. Along each path, the signal experiences a certain attenuation and phase shift. The received signal is the superimposition of multiple alias versions of the transmitted signal. Therefore, the complex baseband signal strength measured at the receiver at a specific time can be written as follows (Yang et al., 2013):

(1) $$V=\sum _{n=1}^N\parallel V_n\parallel e^{-j{\varphi }_n},$$

where $V_n$ and ${\varphi }_n$ are the amplitude and phase of the $n^{th}$ multipath component (note that the modulation scheme of the signal is implicitly considered), and $N$ is the total number of components. On this basis, the recieve signal strength indicator (RSSI) can be written as the received power in decibels (dB):

(2) $$\mathrm{RSSI}=10{\log }_2\left({\parallel V\parallel }^2\right).$$

As the superimposition of multipath components, RSSI not only varies rapidly with propagation distance changing at the order of the signal wavelength but also fluctuates over time, even for a static link. A slight change in specific multipath components may result in significant constructive or destructive multipath components, leading to considerable fluctuations in RSSI.

The essential drawback of RSSI is the failure to reflect the multipath effect. The wireless channel is modeled as a linear temporal filter to fully characterize individual paths, known as channel impulse response (CIR). Under the time-invariant assumption, CIR $h(t)$ is represented as:

(3) $$h(t)=\sum _{n=1}^N{\alpha }_n{e}^{-j{\varphi }_n}\delta (t-{\tau }_n),$$

where ${\alpha }_n$, ${\varphi }_n$, and ${\tau }_n$ are the complex antenuation, phase, and time delay of the $n^{th}$ path, respectively. $N$ is the total number of multipath and $\delta (\cdot )$ is the Dirac delta function. Each impulse represents a delayed multipath component, multiplied by the corresponding amplitude and phase.

In the frequency domain, the multipath causes frequency-selective fading, which is characterized by channel frequency response (CFR). CFR is essentially the Fourier transform of CIR. It consists of both the amplitude response and the phase response. Fig. 2 demonstrate a multipath scenario, the transmitted signal, the received signal, and the illustrative channel responses. Both CIR and CFR depict a small-scale multipath effect and are used for fine-grained channel measurement. Note that the complex amplitudes and antennuation are concerned in CIR and CFR, while another pair of parameters in terms of the signal power is Power Delay Profile (PDP) and Power Spectrum Density (PSD).

Figure 2. Multipath propagations, received signals, and channel responses.

CIR and CFR are measured by decoupling the transmitted signal from the received signal. Specifically, in the time domain, the received signal $r(t)$ is the convolution of transmitted signal $s(t)$ and channel impulse response $h(t)$:

(4) $$r(t)=s(t)\otimes h(t),$$

which indicates the recieved signal is generated from the transmit signal after it propagating from multipath channel.

Similarly, in the frequency domain, the received signal spectrum $R(f)$ is the multiplication of the transmitted signal spectrum $S(f)$ and the channel frequency response $H(f)$:

(5) $$R(f)=S(f)H(f).$$

Note that the $R(f)$ and $S(f)$ are the Fourier transform of the recieved signal $r(t)$ and the transmitted signal $s(t)$ respectively. In this way, Eqn. 4 and Eqn. 5 forms a beautiful “symmetric” relationship.

As demonstrated in Eqn. 4 and Eqn. 5, CIR can be derived from the deconvolution operation of received and transmitted signals, and CFR can be treated as the ratio of the received and the transmitted spectrums. Compared with multiplication, the convolution operation is generally time-consuming. Therefore, in most cases, the device focuses on calculating the CFR, and the CIR can be further derived from the CFR using the inverse Fourier transform (Patwari and Kasera, 2007):

(6) $$h(t)=\frac{1}{P_s}{𝔉}^{-1}\left\{S^{\ast }(f)R(f)\right\},$$

where ${𝔉}^{-1}$ denotes the inverse Fourier transform, $\ast$ is the conjugate operator, and $P_s$ approximates the transmitted signal power.

Although the derivation of CIR and CFR is independent of the modulation scheme, it could be more convenient to implement the process on commercial devices with specific modulation schemes. For the wireless standard where the orthogonal frequency division modulation (OFDM) is adopted (e.g., 802.11a/g/n/ac/ax), the amplitude and phase sampled on each subcarrier can be treated as a sampled version of the signal spectrum $S(f)$. On this basis, a sampled version of $H(f)$ can be easily get from the OFDM receivers.

Recent advances in the wireless community make it possible to get the sampled version of CFR from commercial-off-the-shelf (COTS) Wi-Fi NICs. The extracted CFR are often refered to as a series of complex numbers, which depicts the amplitude and phase of each subcarrier:

(7) $$H(f_j)=\Vert H(f_j)\Vert e^{\mathrm{\angle }H(f_j)},$$

where $H(f_j)$ is a sample at the $j^{th}$ subcarrier, with $\mathrm{\angle }H(f_k)$ denotes its phase. Most research paper treat the CFR at sampled at different subcarriers as the CSI data, which can be written as $𝑯=\left\{H(f_j)\right|j\in [1,J],j\in \mathbb{N}\}$, where $J$ denotes the total number of subcarriers.

The ray-tracing model establishes the relationship between geometric properties of the signal propagation space and the CSI data. Theoretically, by analyzing the multipath signal, various types of geometric information (e.g., the propagation distance, the reflection points) can be derived. Therefore, the ray-tracing model is widely adopted in wireless localization and tracking tasks. In addition, many wireless detection and sensing systems also extract geometric features based on the ray-tracing models, and further put them into machine learning models for regression or classification.

A significant drawback of the ray-tracing model is that it is based on a simple environmental assumption. Therefore, in a complex environment where the signal undergoes diffraction or is disturbed by various types of noise, it becomes difficult to accurately recover all the spatial information by only limited CSI data.

2.2. Scattering Model

The above-mentioned ray-tracing model characterizes signals with multiple propagation paths, which may not apply to rich-scattering environments. To enable wireless sensing in more complex scenarios, some works (Zhang et al., 2018; Wu et al., 2020) establish the scattering model for CSI measurements with Wi-Fi.

The CSI data reported by commodity Wi-Fi NIC can be modeled as:

(8) $$H(f,t)=\sum _{n=1}^N{\alpha }_n(f,t)e^{-j{\varphi }_n(f,t)},$$

where $N$ is the total number of reflection paths, ${\alpha }_n$ is the amplitude attenuation of path $n$, ${\varphi }_n$ is the corresponding phase.

As shown in Fig. 3, the scattering model treats all the objects in indoor environments as scatterers that diffuse the signals to all directions. The CSI observed by the receiver is added up with the portions contributed by the static (furniture, walls, etc.) and dynamic (arms, legs, etc.) scatterers. Intuitively, each scatter is deemed as a virtual Tx. Such modeling can be applied to typical indoor scenarios, where the rooms are crowded with furniture, and signals could propagate in almost all directions. On this basis, we can dismiss the specific signal propagation path and only statistically investigate the relationship between the observed CSI and the moving speed.

Figure 3. Wi-Fi signals in rich-scattering environment

We decompose the observed CSI into the portions contributed by individual scatterers:

(9) $$H(f,t)=\sum _{o\in {\mathrm{\Omega }}_s(t)}{H}_o(f,t)+\sum _{p\in {\mathrm{\Omega }}_d(t)}{H}_p(f,t),$$

where ${\mathrm{\Omega }}_s(t)$ and ${\mathrm{\Omega }}_d(t)$ are the sets of static and dynamic scatterers, $H_p(f,t)$ is the portion of observed CSI contributed by $p^{th}$ dynamic scatterer. For each scatterer, the diffused components undergo anonymous propagation processes and eventually add up at the receiver. Try to a three-dimensional coordinate with its origin at the $p^{th}$ scatterer and its z-axis align with the scatterer’s moving direction, the representation of $H_p(f,t)$ is further decomposed as follows (Zhang et al., 2018):

(10) $$H_p(f,t)={\int }_0^{2\pi }{\int }_0^{\pi }{h}_p(\alpha ,\beta ,f,t)\exp (-jk{v}_p\cos (\alpha )t)d\alpha d\beta ,$$

where $k=2\pi /\lambda$ and $\lambda$ is the signal wavelength, $v_p$ is the speed of the $p^{th}$ dynamic scatterer, $\alpha$ and $\beta$ are the azimuth and elevation angles, $\exp (-jk{v}_i\cos (\alpha )t)$ represents the phase shift of the signal on direction $(\alpha ,\beta )$. $h_p(\alpha ,\beta ,f,t)$ is the portion of observed CSI contributed by $p^{th}$ scatterer on direction $(\alpha ,\beta )$. Inherited from the physical properties of EM waves (Hill, 2009), $h_p(\alpha ,\beta ,f,t)$ can be viewed as a circularly-symmetric Gaussian random variable. For $\forall p,q\in {\mathrm{\Omega }}_d$ and $({\alpha }_1,{\beta }_1)\ne ({\alpha }_2,{\beta }_2)$, $h_p({\alpha }_1,{\beta }_1,f,t)$ is independent of $h_q({\alpha }_2,{\beta }_2,f,t)$.

Based on Eqn. 9 and Eqn. 10, we are able to derive the statistical property by analyzing the autocorrelation function (ACF) of $H(f,t)$ (Hill, 2009):

(11)

where $\mathrm{Cov}[\cdot ,\cdot ]$ is the covariance between two random variables, ${\sigma }_p^2(f)$ is the variance of $h_p(\alpha ,\beta ,f,t)$, $\mathrm{sinc}(kv\tau )=\frac{\sin (kv\tau )}{kv\tau }$, and $\tau$ is the time lag.

For some human activities like fall and walking, we can assume that the dynamic scatterers are mainly on the torso and have similar speeds $v$. Thus, Eqn. 11 is approximated with a much simpler form: ${\rho }_H(f,\tau )\approx \mathrm{sinc}(kv\tau )$. Using this formulation, we have quantitatively established the relationship between the huamn speed and the ACF of CSI. In practice, the speed $v$ is extracted by matching the first peak of the $\mathrm{sinc}(x)$ function and the first peak of ACF:

(12) $$v=\frac{x_0}{k{\tau }_0}=\frac{x_0\lambda }{2\pi {\tau }_0},$$

where $x_0$ is the constant value representing the location of the first peak of $\mathrm{sinc}(x)$ function, ${\tau }_0$ is the time lag corresponding to the first peak of the ACF.

As a straightforward example, we let a volunteer to walk and fall in a heavily furnished room and extract the speed with the CSI collected from a Wi-Fi link located in another room. As is shown in Fig. 4, even with non-line-of-sight (NLoS) occlusions and multipath pollution, the walking speed and fall speed are estimated consistently, demonstrating its robustness to environmental diversity.

Figure 4. The estimated speed and the ACF of CSI signals

The scattering model is generally applicable to various speed-oriented sensing The scattering model generally applies to various speed-oriented sensing tasks such as intrusion detection and fall detection. As a statistical model, it has unique advantages in complex scenarios.

However, the scattering model fails to construct a correspondence between CSI and geometric parameters directly, and thus is generally not applicable to precise localization and tracking tasks.

3. CSI Data Collection

CSI data collection is the first step toward the implementation of a practical wireless sensing system. This section introduces one of the most famous Wi-Fi sensing datasets, the Widar3 dataset, so that beginners can quickly get start with wireless sensing with minimal effort. This section also introduces three of the most famous CSI collection tools that researchers can use to try out collecting CSI data with COTS NICs.

3.1. Widar3 Dataset

Open datasets are essential to provide comprehensive knowledge for model training and a unified benchmark for model comparison. Open datasets are even more necessary in the wireless sensing field because RF signals are more sensitive to devices and deployment environments. However, the absence of high-quality and large-scale datasets has become the bottleneck that hindered the progress of wireless sensing technology. Existing wireless sensing datasets suffer from small scales and limited scenarios in 2019 when we started to build the Widar3 dataset. Widar3 is a wireless sensing dataset for human activity recognition. It is collected from commodity Wi-Fi NICs in the form of RSSI and CSI. It consists of 258,000 instances of hand gestures with a duration of 8,620 minutes and from 75 domains. Widar3 is so far the largest and most comprehensive dataset in this field and receives widespread attention from researchers all over the world. Widar3 dataset is publicly available at IEEE DataPort (official data repository) and continues evolving to contain more types of activities.

3.2. PicoScenes Platform

PicoScenes (Jiang et al., 2021) is a versatile and powerful middleware for CSI-based Wi-Fi sensing research. It is one of the very few tools that support the latest 802.11ac/ax protocols. It supports many prevalent commercial NICs, including Qualcomm Atheros AR9300 (QCA9300), Intel Wireless Link 5300 (IWL5300), Intel AX200 and Intel AX210. PicoScenes supports up to 27 NICs to work concurrently for packet injection and CSI measurement.

PicoScenes is architecturally versatile and flexible. It encapsulates all the low-level features into unified and hardware-independent APIs and exposes them to the upper-level plugin layer. As a result, users can quickly prototype their own measurement plugins.

The data reported by PicoScenes can be parsed in MATLAB as a struct, containing the CSI data of different packets, subcarriers, and antennas. The struct also includes other helpful information such as timestamps, RSSI, and the signal-to-noise ratio (SNR).

The homepage of this tool can be accessed for more detailed information.

3.3. Intel 5300 NIC CSI Tool

This CSI Tool (Halperin et al., 2011) is built upon the Intel WiFi Wireless Link 5300 802.11n MIMO radios, using a modified firmware and the open-source Linux wireless driver. It includes all the software and scripts required to collect, read, and parse CSI.

The IWL5300 provides 802.11n CSI of 30 subcarrier groups. Each group contains 2 adjacent subcarriers given 20 MHz bandwidth or 4 given 40 MHz bandwidth. Each CSI sample is a complex number, with a signed 8-bit resolution for both real and imaginary parts. One CSI record is a $A\times 30$ matrix, where $M$ is the number of pairs of transmitting and receiving antennas.

The homepage of this tool can be accessed for detailed information.

3.4. Atheros CSI Tool

Atheros CSI Tool (Xie et al., 2018) is an open-source 802.11n measurement and experimentation tool. It enables the extraction of detailed PHY wireless communication information from the Atheros WiFi NICs, including the Channel State Information (CSI), the received packet payload, and other information (the time stamp, the RSSI of each antenna, the data rate, etc.). Atheros-CSI-Tool is built on top of ath9k, an open-source Linux kernel driver supporting Atheros 802.11n PCI/PCI-E chips. Thus, this tool theoretically supports all types of Atheros 802.11n WiFi chipsets. We have tested it on Atheros AR9580, AR9590, AR9344, and QCA9558. Furthermore, Atheros CSI Tool is open source, and all functionalities are implemented in software without any modification to the firmware. Therefore, one can extend the functionalities of Atheros CSI Tool with their own codes under the GPL license.

Atheros-CSI-Tool works on various Linux distributions, e.g., Ubuntu, OpenWRT, Linino, etc.. Different Linux distribution works with different hardware. Ubuntu works for personal computers like laptops or desktops. OpenWRT works for embedded devices such as WiFi routers. Linino works for IoT devices, such as Arduino YUN. The official website provides the source code for the Ubuntu version and OpenWRT version of the Atheros CSI tool.

The homepage of this tool can be accessed for detailed information.

4. CSI Feature Extraction

The CSI features lay the fundation of wireless sensing. In particular, for different sensing tasks, choosing the most appropriate features can effectively improve the system performance. In addition, the quality of the extracted features determines the effectiveness of the sensing system.

For ease of illustration, code implementation used for feature extraction is provided below, which is a main function to call different function in the following subsections.

⬇

%{

CSI Feature Extraction for Wi-Fi sensing.

- Input: csi data used for calibration, and csi data that need to be sanitized.

- Output: sanitized csi data.

To use this script, you need to:

1. Make sure the csi data have been saved as .mat files.

2. Check the .mat file to make sure they are in the correct form.

3. Set the parameters.

Note that in this algorithm, the csi data should be a 4-D tensor with the size of [T S A L]:

- T indicates the number of csi packets;

- S indicates the number of subcarriers;

- A indicates the number of antennas (i.e., the STS number in a MIMO system);

- L indicates the number of HT-LTFs in a single PPDU;

Say if we collect a 10 seconds csi data steam at 1 kHz sample rate (T = 10 * 1000), from a 3-antenna AP (A = 3), with 802.11n standard (S = 57 subcarrier), without any extra spatial sounding (L = 1), the data size should be [10000 57 3 1].

%}

clear all;

addpath(genpath(pwd));

%% 0. Set parameters.

% Path of the raw CSI data.

src_file_name = ’./data/csi_src_test.mat’;

% Speed of light.

global c;

c = physconst(’LightSpeed’);

% Bandwidth.

global bw;

bw = 20e6;

% Subcarrier frequency.

global subcarrier_freq;

subcarrier_freq = linspace(5.8153e9, 5.8347e9, 57);

% Subcarrier wavelength.

global subcarrier_lambda;

subcarrier_lambda = c ./ subcarrier_freq;

% Antenna arrangement.

antenna_loc = [0, 0, 0; 0.0514665, 0, 0; 0, 0.0514665, 0]’;

% Set the linear range of the CSI phase, which varies with NIC types.

linear_interval = (20:38)’;

%% 1. Read the csi data for calibration and sanitization.

% Load the raw CSI data.

csi_src = load(src_file_name).csi; % Raw CSI.

%% 2. Perform various wireless sensing tasks.

% Test example 1: angle/direction estimation with imperfect CSI.

[packet_num, subcarrier_num, antenna_num, ~] = size(csi_src);

aoa_mat = naive_aoa(csi_src, antenna_loc, zeros(3, 1));

aoa_gt = [0; 0; 1];

error = mean(acos(aoa_gt’ * aoa_mat));

disp("Angle estimation error: " + num2str(error));

% Test example 2: distance estimation with CSI.

tof_mat = naive_tof(csi_src);

est_dist = mean(tof_mat * c, ’all’);

disp("The ground truth distance is: 10 m");

disp("The estimated distance is: " + num2str(est_dist) + " m");

4.1. Time of Flight

Figure 5. The relationship between ToF and CIR.

ToF is the time duration the signal propagates from the transmitter to the receiver along a specific path. Given the frequency $f$, the phase shift introduced by the ToF $\tau$ is:

(13) $${\varphi }_{\mathrm{ToF}}=-2\pi f\tau$$

As the superimposition of multipath signals, CSI can be represented based on the ray-tracing model:

(14) $$H(f)=\sum _{n=1}^N{\alpha }_n{e}^{-j2\pi f{\tau }_n}$$

where $N$ is the total number of multipath, and ${\alpha }_n$ and ${\tau }_n$ are the complex attenuation factor and time of flight (ToF) for the $n^{th}$ path, respectively. Theoretically, the ToF of all paths can be identified in CIR, which can be calculated by applying the inverse Fourier transform to CSI samples of all subcarriers. However, since the transmitter and the receiver lack synchronization, non-zero temporal shifts exist in CIR, and the absolute ToF is typically not accurate enough. The limited bandwidth also constrains the time resolution, causing meter-level ToF ambiguity The relationship between signal propagation path, ToF, and CIR is shown in Figure 5.

The following function naive_tof intends to extract the ToF of the strongest path (typically the shortest path) based on inverse Fourier transform.

⬇

function [tof_mat] = naive_tof(csi_data)

% naive_tof

% Input:

% - csi_data is the CSI used for ranging; [T S A E]

% - ifft_point and bw are the IFFT and bandwidth parameters;

% Output:

% - tof_mat is the rough time-of-flight estimation result; [T A]

global c, bw;

[pakcet_num, subcarrier_num, antenna_num, extra_num] = size(csi_data);

ifft_point = power(2, ceil(log2(subcarrier_num)));

% Get CIR from each packet and each antenna by ifft(CFR);

cir_sequence = zeros(packet_num, antenna_num, extra_num, ifft_point);

for p = 1:packet_num

for a = 1:antenna_num

for e = 1:extra_num

cir_sequence(p, a, e, :) = ifft(csi_data(p, :, a, e), ifft_point);

end

end

end

cir_sequence = squeeze(mean(cir_sequence, 4)); % [T ifft_point A]

half_point = ifft_point / 2;

half_sequence = cir_sequence(:, 1:half_point, :); % [T half_point A]

peak_indices = zeros(packet_num, antenna_num); % [T A]

for p = 1:pakcet_num

for a = 1:antenna_num

[~, peak_indices(p, a)] = max(half_sequence, [], 2);

end

end

% Calculate ToF;

tof_mat = peak_indices .* subcarrier_num ./ (ifft_point .* bw); % [T A]

end

4.2. Angle of Arrival and Angle of Departure

Figure 6. Angle of Arrival and Angle of Departure.

When a NIC equipes with multiple antennas, a local coordinate at the device can be created. As shown in Figure 6, for a transmitter, the angle of departure (AoD) $\phi$ represents the direction in the local coordinate along which the transmitted signal is emitted. For a receiver, the angle of arrival (AoA) $\theta$ represents the direction in the local coordinate along which the received signal is captured. Since the antennas are spatially separated, non-zero phase shifts between antennas are introduced. The phase shifts depend on the AoA/AoD. Specifically, suppose the relative location between two antennas is $𝚫𝒍=({\mathrm{\Delta }}_x,{\mathrm{\Delta }}_y)$ and the unit direciton vector of AoA is $𝒆=(cos\theta ,\sin \theta )$, the phase shift between the two antennas is:

(15) $${\varphi }_{\mathrm{AoA}}=\frac{2\pi }{\lambda }𝚫𝒍\cdot 𝒆$$

Then CSI can be modeled as:

(16) $$H(k)=\sum _{n=1}^N{\alpha }_n{e}^{-j\frac{2\pi }{\lambda }𝚫𝒍\cdot 𝒆}$$

where $k$ represents the $k^{th}$ antenna at the receiver. The same model applies to the AoD at the trasmitter side and the 3D space with azimuth and elevation angles.

In practice, algorithms such as Capon (Capon, 1969) and MUSIC (Schmidt, 1986) can be used to estimate the AoA/AoD of multiple paths from the CSI of the antenna array.

MUSIC analyses the incident signals on multiple antennas to find out the AoA of each signal. Specifically, suppose $D$ signals $F_1,\mathrm{\cdots },F_D$ arrive from directions ${\theta }_1,\mathrm{\cdots },{\theta }_D$ at $M>D$ antennas. The received signal at the $k^{th}$ antenna element, denoted as $X_k$, is a linear combination of the $D$ incident wavefronts and noise $W_k$:

or

(17) $$𝑿=𝑨𝑭+𝑾$$

where $𝐚({\theta }_k)$ is the array steering vector that characterizes added phase (relative to the first antenna) of each receiving component at the $k^{th}$ antenna. $𝑨$ is the matrix of steer vectors. As shown in Figure 6, for a linear antenna array with elements well synchronized,

(18) $$𝒂(\theta )=\left[\begin{array}{c}1\\ e^{-j\frac{2\pi }{\lambda }𝚫𝒍\mathbf{(}\mathrm{𝟏}\mathbf{)}\cdot 𝒆}\\ e^{-j\frac{2\pi }{\lambda }𝚫𝒍\mathbf{(}\mathrm{𝟐}\mathbf{)}\cdot 𝒆}\\ \mathrm{⋮}\\ e^{-j\frac{2\pi }{\lambda }𝚫𝒍\mathbf{(}𝑴\mathbf{-}\mathrm{𝟏}\mathbf{)}\cdot 𝒆}\end{array}\right]$$

Suppose $W_k\sim N(0,{\sigma }^2)$, and $F_k$ is a wide-sense stationary process with zero mean value, the $M\times M$ covariance matrix of the received signal vector $𝑿$ is:

(19) $𝑺$ $=\overline{𝑿{𝑿}^{\mathrm{H}}}$

$=𝑨\overline{𝑭{𝑭}^{\mathrm{H}}}{𝑨}^{\mathrm{H}}+\overline{𝑾{𝑾}^{\mathrm{H}}}$

$=𝑨𝑷{𝑨}^{\mathrm{H}}+{\sigma }^2𝑰$

where $𝑷$ is the covariance matrix of transmission vector $𝑭$. The notation ${(\cdot )}^{\mathrm{H}}$ represents conjugate transpose and $\overline{(\cdot )}$ represents expectation.

The covariance matrix $𝑺$ has $M$ eigenvalues ${\lambda }_1,\mathrm{\cdots },{\lambda }_M$ associated with $M$ eigenvectors ${𝒆}_1,{𝒆}_2,\mathrm{\cdots },{𝒆}_M$. Sorted in a non-descending order, the smallest $M-D$ eigenvalues correspond to the noise while the rest $D$ correspond to the $D$ incident signals. In other word, the $M$-dimension space can be divided into two orthogonal subspace, the noise subspace ${𝑬}_N$ expanded by eigenvectors ${𝒆}_1,\mathrm{\cdots },{𝒆}_{M-D}$, and the signal subspace ${𝑬}_S$ expanded by eigenvectors ${𝒆}_{M-D+1},\mathrm{\cdots },{𝒆}_M$ (or equivalently $D$ array steering vector $𝒂({\theta }_1),\mathrm{\cdots },𝒂({\theta }_D)$).

To solve for the array steering vectors (thus AoA), MUSIC plots the reciprocal of squared distance $Q(\theta )$ for points along the $\theta$ continue to the noise subspace as a function of $\theta$:

(20) $$Q(\theta )=\frac{1}{{𝒂}^{\mathrm{H}}(\theta ){𝑬}_N{𝑬}_N^{\mathrm{H}}𝒂(\theta )}$$

This yields peaks in $Q(\theta )$ at the bearing of incident signals. It is similar to apply MUSIC algorithm for AoD spectrum estimation.

The following function naive_aoa intends to estimate the 3D AoA based on the phase difference, which is similar to Eqn. 15. Note that the following algorithm only considers one path, and thus cannot be applied to mutlipath signals.

⬇

function [aoa_mat] = naive_aoa(csi_data, antenna_loc, est_rco)

% naive_aoa

% Input:

% - csi_data is the CSI used for angle estimation; [T S A E]

% - antenna_loc is the antenna location arrangement with the first antenna as a reference; [3 A]

% - est_rco is the estimated radio chain offset; [A 1]

% Output:

% - aoa_mat is the angle estimation result; [3 T]

global subcarrier_lambda;

[packet_num, subcarrier_num, antenna_num, extra_num] = size(csi_data);

csi_phase = unwrap(angle(csi_data), [], 2); % [T S A E]

% Get the antenna vector and its length.

ant_diff = antenna_loc(:, 2:end) - antenna_loc(:, 1); % [3 A-1]

ant_diff_length = vecnorm(ant_diff); % [1 A-1]

ant_diff_normalize = ant_diff ./ ant_diff_length; % [3 A-1]

% Calculate the phase difference.

phase_diff = csi_phase(:, :, 2:end, :) - csi_phase(:, :, 1, :) - permute(est_rco(2:end, :), [4 3 1 2]); % [T S A-1 E]

phase_diff = unwrap(phase_diff, [], 2);

phase_diff = mod(phase_diff + pi, 2 * pi) - pi;

% Broadcasting is performed, get the value of cos(theta) for each packet and each antenna pair.

cos_mat = subcarrier_lambda .* phase_diff ./ (2 .* pi .* permute(ant_diff_length, [3 1 2])); % [T S A-1 E]

cos_mat_mean = squeeze(mean(cos_mat, [2 4])); % [T A-1]

% Symbolic nonlinear optimization are performed.

syms x y

% aoa_sol = [x;y;(1-sqrt(x^2 + y^2))];

aoa_init = [sqrt(1/3);sqrt(1/3);sqrt(1/3)];

aoa_mat_sol = zeros(3, packet_num);

options = optimoptions(’lsqnonlin’, ’Algorithm’, ’levenberg-marquardt’, ’Display’, ’none’);

parfor p = 1:packet_num

cur_nonlinear_func = @(aoa_sol)ant_diff_normalize’ * aoa_sol - cos_mat_mean(p, :)’;

cur_aoa_sol = lsqnonlin(cur_nonlinear_func, aoa_init, [], [], options);

aoa_mat_sol(:, p) = cur_aoa_sol;

end

aoa_mat = aoa_mat_sol ./ vecnorm(aoa_mat_sol); % [3 T]

end

4.3. Phase Shift Spectrum

Figure 7. Phase shift spectrum (or Doppler spectrum) of three different moving path.

Non-zero phase shift $\mathrm{\Delta }\varphi$ across different packets is caused by the relative movement of the transmitter, receiver, or objects in the propagation path of the signal. It equals the changing rate of the path length of the signal. When multiple packets are received in sequence, the CSI corresponding to the $i^{th}$ received packet is:

(21) $$H(i)=\sum _{n=1}^N{\alpha }_n(i)e^{j{\varphi }_n(i)},$$

where ${\varphi }_n$ is the phase of the $n^{th}$ path (Chi et al., 2022). Extract the phase of the the $n^{th}$ path in packet $i$ and $i+1$ respectively, and calculate the phase shift as:

(22) $$\mathrm{\Delta }{\varphi }_n(i)={\varphi }_n(i+1)-{\varphi }_n(i).$$

Intuitively, the phase difference indicates the distance change of the $n^{th}$ path between two consecutive packets: $\mathrm{\Delta }{d}_n(i)=\frac{\mathrm{\Delta }{\varphi }_n(i)}{2\pi }\lambda$.

Take a step further, and apply the short-time Fourier transform (STFT) within a sliding window, we can get the spectrum as shown in Figure 7. The frequency axis reveals the change rate of consecutive CSI data, and implicitly contains the path length change rate. Figure 7 demonstrates the phase shift spectrum (or the Doppler Spectrum) for three different moving paths.

The following function naive_stft calculates the short-time Fourier transform of a series of CSI data. The generated spectrum can be used effectively for many wireless sensing tasks, like gesture recognition and fall detection.

⬇

function stft_mat = naive_spectrum(csi_data, sample_rate, visable)

% naive_spectrum

% Input:

% - csi_data is the CSI used for STFT spectrum generation; [T S A L]

% - sample_rate determines the resolution of the time-domain and

% frequency-domain;

% Output:

% - stft_mat is the generated STFT spectrum; [sample_rate/2 T]

% Conjugate multiplication.

csi_data = mean(csi_data .* conj(csi_data), [2 3 4]);

% Calculate the STFT and visualization.

stft_mat = stft(csi_data, sample_rate);

% Visualization (optional).

if visable

stft(csi_data, sample_rate);

end

end

4.4. Body-coordinate Velocity Profile

The limitation of the aforementioned spectrum is that, even the spectrum corresponding to the same activity will be different when the user moves at different locations or orientations relative to the Wi-Fi links. To resolve this problem, Widar3.0 (Zheng et al., 2019) proposes a domain-independent signal feature BVP (body-coordinate velocity profile) to characterize human activities.

Figure 8. Relationship between the BVP and Doppler spectrum. Each velocity component in BVP is projected onto the normal direction of a link, and contributes to the power of the corresponding radial velocity component in the Doppler spectrum.

The basic idea of BVP is shown in Fig. 8. A BVP $𝑽$ is quantized as a discrete matrix with dimension as velocity components decomposed along each axis of the body coordinates. For convenience, we establish the local body coordinates whose origin is the location of the person and positive x-axis aligns with the orientation of the person. The person’s location and orientation should be provided manually. Currently, it is assumed that the global location and orientation of the person are available. Then the known global locations of wireless transceivers can be transformed into the local body coordinates. Thus, for better clarity, all locations and orientations used in the following derivation are in the local body coordinates. Suppose the locations of the transmitter and the receiver of the $i^{th}$ link are ${\overrightarrow{l}}_t^{(i)}=(x_t^{(i)},y_t^{(i)})$, ${\overrightarrow{l}}_r^{(i)}=(x_r^{(i)},y_r^{(i)})$, respectively, then any velocity components $\overrightarrow{v}=(v_x,v_y)$ around the human body (i.e., the origin) will contribute its signal power to some frequency component, denoted as $f^{(i)}(\overrightarrow{v})$, in the Doppler spectrum of the $i^{th}$ link (Qian et al., 2017):

(23) $$f^{(i)}(\overrightarrow{v})=a_x^{(i)}{v}_x+a_y^{(i)}{v}_y.$$

$a_x^{(i)}$ and $a_y^{(i)}$ are coefficients determined by locations of the transmitter and the receiver:

(24) $a_x^{(i)}={\displaystyle \frac{1}{\lambda }}({\displaystyle \frac{x_t^{(i)}}{{\Vert {\overrightarrow{l}}_t^{(i)}\Vert }_2}}+{\displaystyle \frac{x_r^{(i)}}{{\Vert {\overrightarrow{l}}_r^{(i)}\Vert }_2}}),$

$a_y^{(i)}={\displaystyle \frac{1}{\lambda }}({\displaystyle \frac{y_t^{(i)}}{{\Vert {\overrightarrow{l}}_t^{(i)}\Vert }_2}}+{\displaystyle \frac{y_r^{(i)}}{{\Vert {\overrightarrow{l}}_r^{(i)}\Vert }_2}}),$

where $\lambda$ is the wavelength of Wi-Fi signal. As static components with zero Doppler spectrum (e.g., the line of sight signals and dominant reflections from static objects) are filtered out before the Doppler spectrum are calculated, only signals reflected by the person are retained. Besides, when the person is close to the Wi-Fi link, only signals with one-time reflection have prominent magnitudes (Qian et al., 2018). Thus, Eqn. 23 holds valid for the gesture recognition scenario. From the geometric view, Eqn. 23 means that the 2-D velocity vector $\overrightarrow{v}$ is projected on a line whose direction vector is ${𝒅}^{(i)}=(-a_y^{(i)},a_x^{(i)})$. Suppose the person is on an ellipse curve whose foci are the transmitter and the receiver of the $i^{th}$ link, then $d^{(i)}$ is indeed the average direction of the ellipse at the person’s location. Fig. 8 shows an example where the person generates three velocity components ${\overrightarrow{v}}_j,j=1,2,3$, and projection of the velocity components on the Doppler spectrum of three links.

Since coefficients $a_x^{(i)}$ and $a_y^{(i)}$ only depend on the location of the $i^{th}$ link, the relation of projection of the BVP on the $i^{th}$ link is fixed. Specifically, an assignment matrix ${𝑨}_{F\times N^2}^{(i)}$ can be defined:

(25)

where $f_j$ is the $j^{th}$ frequency sampling point in the Doppler spectrum, and ${\overrightarrow{v}}_k$ is velocity component corresponding to the $k^{th}$ element of the vectorized BVP $𝑽$. Thus, the relation between Doppler spectrum profile of the $i^{th}$ link and the BVP can be modeled as:

(26) $${𝑫}^{(i)}=c^{(i)}{𝑨}^{(i)}𝑽$$

where $c^{(i)}$ is the scaling factor due to propagation loss of the reflected signal.

Due to the sparsity of BVP, compressed sensing (Donoho, 2006) technique can be adopted to formulate the estimation of BVP as an $l_0$ optimization problem:

(27) $${\min }_V\sum _{i=1}^M|\mathrm{EMD}({𝑨}^{(i)}𝑽,{𝑫}^i)|+\eta {\Vert V\Vert }_0,$$

where $M$ is the number of Wi-Fi links. The sparsity of the number of the velocity components is coerced by the term $\eta {\Vert V\Vert }_0$, where $\eta$ represents the sparsity coefficients and $\parallel \cdot {\parallel }_0$ is the number of non-zero velocity components.

$\mathrm{EMD}(\cdot ,\cdot )$ is the Earth Mover’s Distance (Rubner and Tomasi, 2001) between two distributions. The selection of EMD rather than Euclidean distance is mainly due to two reasons. First, the quantization of BVP introduces approximation error, i.e., projection of velocity components to the Doppler spectrum bin might be adjacent to the true one. Such quantization error can be relieved by EMD, which takes the distance between bins into consideration. Second, there are unknown scaling factors between the BVP and Doppler spectrum, making the Euclidean distance inapplicable.

5. CSI Sanitization

The wireless sensing models and features described in previous sections are consistent with the EM propagation theory and geometry. However, they don’t consider the various types of noise caused by imperfect implementations of transceiver hardware (Xie et al., 2018; Xie et al., 2019). This section focuses on various CSI error sources and the corresponding error cancellation algorithms.

For ease of illustration, code implementation for sanitization is provided, which is a main function to call different error cancellation functions and test their performance.

⬇

%{

CSI Sanitization Algorithm for Wi-Fi sensing.

- Input: csi data used for calibration, and csi data that need to be sanitized.

- Output: sanitized csi data.

To use this script, you need to:

1. Make sure the csi data have been saved as .mat files.

2. Check the .mat file to make sure they are in the correct form.

3. Set the parameters.

Note that in this algorithm, the csi data should be a 4-D tensor with the size of [T S A L]:

- T indicates the number of csi packets;

- S indicates the number of subcarriers;

- A indicates the number of antennas (i.e., the STS number in a MIMO system);

- L indicates the number of HT-LTFs in a single PPDU;

Say if we collect a 10 seconds csi data steam at 1 kHz sample rate (T = 10 * 1000), from a 3-antenna AP (A = 3), with 802.11n standard (S = 57 subcarrier), without only one spatial stream (L = 1), the data size should be [10000 57 3 1].

%}

clear all;

addpath(genpath(pwd));

%% 0. Set parameters.

% Path of the calibration data;

calib_file_name = ’./data/csi_calib_test.mat’;

% Path for storing the generated calibration templated.

calib_template_name = ’./data/calib_template_test.mat’;

% Path of the raw CSI data.

src_file_name = ’./data/csi_src_test.mat’;

% Path for storing the sanitized CSI data.

dst_file_name = ’./data/csi_dst_test.mat’;

% Speed of light.

global c;

c = physconst(’LightSpeed’);

% Bandwidth.

global bw;

bw = 20e6;

% Subcarrier frequency.

global subcarrier_freq;

subcarrier_freq = linspace(5.8153e9, 5.8347e9, 57);

% Subcarrier wavelength.

global subcarrier_lambda;

subcarrier_lambda = c ./ subcarrier_freq;

% Antenna arrangement.

antenna_loc = [0, 0, 0; 0.0514665, 0, 0; 0, 0.0514665, 0]’;

% Set the linear range of the CSI phase, which varies with NIC types.

linear_interval = (20:38)’;

%% 1. Read the csi data for calibration and sanitization.

% Load the calibration data.

csi_calib = load(calib_file_name).csi; % CSI for calibration.

% Load the raw CSI data.

csi_src = load(src_file_name).csi; % Raw CSI.

%% 2. Choose different functions according to your task.

% Use cases:

% Make calibration template.

csi_calib_template = set_template(csi_calib, linear_interval, calib_template_name);

% Directly load the generated template.

csi_calib_template = load(calib_template_name).csi;

% Remove the nonlinear error.

csi_remove_nonlinear = nonlinear_calib(csi_src, csi_calib_template);

% Remove the STO (a.k.a SFO and PBD) by conjugate mulitplication.

csi_remove_sto = sto_calib_mul(csi_src);

% Remove the STO (a.k.a SFO and PBD) by conjugate division.

csi_remove_sto = sto_calib_div(csi_src);

% Estimate the CFO by frequency tracking.

est_cfo = cfo_calib(csi_src);

% Estimate the RCO.

est_rco = rco_calib(csi_calib);

%% 3. Save the sanitized data as needed.

csi = csi_remove_sto;

save(dst_file_name, ’csi’);

%% 4. Perform various wireless sensing tasks.

% Test example 1: angle/direction estimation with imperfect CSI.

[packet_num, subcarrier_num, antenna_num, ~] = size(csi_src);

est_rco = rco_calib(csi_calib);

zero_rco = zeros(antenna_num, 1);

aoa_mat_error = naive_aoa(csi_src, antenna_loc, zero_rco);

aoa_mat = naive_aoa(csi_src, antenna_loc, est_rco);

aoa_gt = [0; 0; 1];

error_1 = mean(acos(aoa_gt’ * aoa_mat_error));

error_2 = mean(acos(aoa_gt’ * aoa_mat));

disp("Angle estimation error (in deg) without RCO removal: " + num2str(error_1));

disp("Angle estimation error (in deg) with RCO removal: " + num2str(error_2));

% Test example 2: intrusion detection with CSI.

csi_sto_calib = sto_calib_div(csi_src);

intrusion_flag_raw = naive_intrusion(csi_src, 3);

intrusion_flag_sanitized = naive_intrusion(csi_sto_calib, 3);

disp("Intrusion detection result without SFO/PDD removal: " + num2str(intrusion_flag_raw));

disp("Intrusion detection result with SFO/PDD removal: " + num2str(intrusion_flag_sanitized));

5.1. Nonlinear Amplitude and Phase

The nonlinear amplitude and phase errors are caused by the imperfect analog domain filter implementation inside the hardware. Specifically, it causes the extracted CSI amplitude and phase to be equivalently processed by a nonlinear function. Let $f(\cdot )$ and $g(\cdot )$ be the nonlinear modes of CSI amplitude and phase, respectively, the errorous CSI can be written as:

(28) $$\tilde{H}(i,j,k)=\sum _{n=1}^{N}f({\alpha }_n)e^{-jg({\varphi }_n(i,j,k))}+N(i,j,k).$$

Specifically, during the OFDM modulation process, each subcarrier of should have the same gain. In other words, the amplitude-frequency characteristic of CSI should be a horizontal straight line when a coaxial cable is used to connect the transceiver ports. However, actual measurements show that even without the multipath radio channel, there is still a similar ”frequency selective fading” characteristic, i.e., the gain of each frequency band is different, showing an M-shaped amplitude-frequency characteristic curve. Similarly, when using a coaxial cable to connect the transceiver port, the CSI phase-frequency characteristics obtained by the NIC are not an ideal straight line with slope, but an S-shaped curve with certain nonlinearity.

After extensive research and experiments, we have observed two facts.

• First, for a specific type of NIC, the nonlinear amplitude/phase error of CSI is fixed. This means that the correction task can be accomplished if we use a known length coaxial cable connected to the transceiver port. Before performing the sensing task, measure a representative set of CSI amplitude and phase, record the nonlinear characteristics, and eliminate the nonlinearity in the subsequent steps.

• Second, we observe that the middle part of the subcarrier is free of nonlinear errors, and the nonlinear characteristics of both sides are also fixed.

Therefore, the function below the code performs the following steps to tackle the CSI nonlinearity:

1. Read in a set of CSI data measured using a coaxial cable.

2. Get its amplitude and phase.

3. Normalize its amplitude and record it as the amplitude template.

4. Unwrap the phases, then perform a linear fit to the middle part of its subcarriers. Subtract the linear fit result to obtain the nonlinear phase error template.

Finally, the nonlinear amplitude and nonlinear phase components are saved in the form of $\alpha e^{j\varphi }$, which indicates the nonlinear error template of a specific type of NIC.

⬇

function [csi_calib_template] = set_template(csi_calib, linear_interval, calib_template_name)

% set_template

% Input:

% - csi_calib is the reference csi at given distance and angle; [T S A L]

% - linear_interval is the linear range of the csi phase, which varies across different types of NICs;

% - calib_template_name is the saved path of the generated template;

% Output:

% - csi_calib_template is the generated template for csi calibration; [1 S A L]

[packet_num, subcarrier_num, antenna_num, extra_num] = size(csi_calib);

csi_amp = abs(csi_calib); % [T S A L]

csi_phase = unwrap(angle(csi_calib), [], 2); % [T S A L]

csi_amp_template = mean(csi_amp ./ mean(csi_amp, 2), 1); % [1 S A L]

nonlinear_phase_error = zeros(size(csi_calib)); % [T S A L]

for p = 1:packet_num

for a = 1:antenna_num

for e = 1:extra_num

linear_model = fit(linear_interval, squeeze(csi_phase(p, linear_interval, a, e))’, ’poly1’);

nonlinear_phase_error(p, :, a, e) = csi_phase(p, :, a, e) - linear_model(1:subcarrier_num)’;

end

end

end

csi_phase_template = mean(nonlinear_phase_error, 1); % [1 S A L]

csi_phase_template(1, linear_interval, :, :) = 0;

csi_calib_template = csi_amp_template .* exp(1i * csi_phase_template); % [1 S A L]

csi = csi_calib_template;

save(calib_template_name, ’csi’); % [1 S A L]

end

After getting the nonlinear error template, the sanitization process begins. For the raw CSI data csi_data collected in real time, we divide it by the error template csi_calib. This operation is equivalent to “dividing the original amplitude by the normalized nonlinear amplitude” and “subtracting the original phase from the nonlinear phase”, so that the returned CSI data csi_proc has sanitized amplitude and phase.

⬇

function [csi_remove_nonlinear] = nonlinear_calib(csi_calib, csi_calib_template)

% nonlinear_calib

% Input:

% - csi_src is the raw csi which needs to be calibrated; [T S A L]

% - csi_calib_template is the reference csi for calibration; [1 S A L]

% Output:

% - csi_remove_nonlinear is the csi data in which the nonlinear error has been eliminated; [T S A L]

csi_amp = abs(csi_calib); % [T S A L]

csi_phase = unwrap(angle(csi_calib), [], 2); % [T S A L]

csi_unwrap = csi_amp .* exp(1i * csi_phase); % [T S A L]

% Broadcasting is performed.

csi_remove_nonlinear = csi_unwrap ./ csi_calib_template;

end

5.2. Automatic Gain Control Uncertainty

Automatic gain control (AGC) induces a random gain $\beta$ in each received CSI packet.

(29) $$\tilde{H}(i,j,k)=\sum _{n=1}^N{\beta }_i{\alpha }_n{e}^{-j{\varphi }_n(i,j,k)}+N(i,j,k).$$

There are two ways to eliminate the AGC error: 1) Disable the AGC function of the wireless driver; 2) Compensate the amplitude of the measured CSI based on the reported AGC.

⬇

function [csi_remove_agc] = agc_calib(csi_src, csi_agc)

% rco_calib

% Input:

% - csi_src is the raw csi which needs to be calibrated; [T S A L]

% - csi_agc is the AGC amplitude reported by the NIC; [T, 1]

% Output:

% - csi_remove_agc is the csi data in which the AGC uncertainty has been eliminated; [T S A L]

% Broadcasting is performed.

csi_remove_agc = csi_src ./ csi_agc;

end

5.3. Radio Chain Offset

Radio chain offset (RCO) is the random phase variation ${ϵ}_{\varphi }$ introduced between different Tx/Rx chains (transceiver antenna pairs). The RCO is reset each time the NIC is powered up.

(30) $${\tilde{\varphi }}_n(i,j,k)=2\pi (f_c+\mathrm{\Delta }{f}_j+f_D){\tau }_n(i,j,k)+{ϵ}_{\varphi },$$

RCO induces a biast of the phase-frequency characteristic curve. It undermines the accuracy of features such as AoA or ToF. Fortunately, we found that this type of phase deviation is consistent between each successive packet sent and therefore doesn’t affect the performance of temporal tracking or sensing, and that this type of error can be eliminated by the following steps:

1. Once the NIC power-up, connect the transceiver port using a coaxial cable of known length and record the phase information calib_phase.

2. during subsequent measurements, subtracting this phase information from the measured phase csi_phase.

⬇

function [est_rco] = rco_calib(csi_calib)

% rco_calib

% Input:

% - csi_calib is the reference csi at given distance and angle; [T S A L]

% Output:

% - est_rco is the estimated RCO; [A 1]

antenna_num = size(csi_calib, 3);

csi_phase = unwrap(angle(csi_calib), [], 1); % [T S A L]

avg_phase = zeros(antenna_num, 1);

for a = 1:antenna_num

avg_phase(a, 1) = mean(csi_phase(:, :, a, 1), ’all’);

end

est_rco = avg_phase - avg_phase(1, 1);

end

5.4. Central Frequency Offset

Central frequency offset (CFO), which is caused by the frequency desynchronization on both sides of the transceiver, leads to random frequency shift ${ϵ}_f$ of each received CSI.

(31) $${\tilde{\varphi }}_n(i,j,k)=2\pi (f_c+\mathrm{\Delta }{f}_j+f_D+{ϵ}_f){\tau }_n(i,j,k)={\varphi }_n(i,j,k)+{ϵ}_f{\tau }_n(i,j,k).$$

The CFO induces an extra bias (i.e., an overall up and down shift) of the phase-frequency characteristic curve.

To eliminate CFO, we need to insert multiple HT-LTFs in the same PPDU (Wi-Fi data frame), and therefore obtaining multiple CSI measurements. Since the time interval between multiple HT-LTFs is strictly controlled to $4\mu s$ according to the 802.11 protocol, the phase difference between two HT-LTFs is induced by the CFO within $\mathrm{\Delta }t=4\mu s$. Thus, the approximate value of the CFO can be recovered.

⬇

function [est_cfo] = cfo_calib(csi_src)

% cfo_calib

% Input:

% - csi_src is the csi data with two HT-LTFs; [T S A L]

% Output:

% - est_cfo is the estimated frequency offset;

delta_time = 4e-6;

phase_1 = angle(csi_src(:, :, :, 1));

phase_2 = angle(csi_src(:, :, :, 2));

phase_diff = mean(phase_2 - phase_1, 3); % [T S A 1]

est_cfo = mean(phase_diff ./ delta_time, 2);

end

5.5. Sampling Frequency Offset and Packet Detection Delay

The sampling frequency offset (SFO), which appears to be an error in frequency domain, are generally considered as an equivalent time shift due to frequency asynchrony. The packet detection delay (PDD) is a time delay. Therefore, despite their distinct causes, they are often discussed together as a “time offset” together.

(32) $${\tilde{\varphi }}_n(i,j,k)=2\pi (f_c+\mathrm{\Delta }{f}_j+f_D)({\tau }_n(i,j,k)+{ϵ}_t)={\varphi }_n(i,j,k)+2\pi (f_c+\mathrm{\Delta }{f}_j+f_D){ϵ}_t.$$

This type of error is critical because the time delay can be confused with real ToF and thus affect the accuracy of the ranging accuarcy. Specifically, this deviation will be characterized in the phase-frequency characteristic curve as a change in slope, since this time deviation causes different phase changes with different sub-bands $\mathrm{\Delta }{f}_j$.

Currently, there is no “perfect algorithm” to solve this type of error. Conjugate multiplication and division are the only two methods to eliminate the SFO and PDD. The code of both of them are listed below. By appling the conjugate multiplication or division, the ${ϵ}_t$ is eliminated, at the cost of losing absolute ToF measurement.

⬇

function [csi_remove_sto] = sto_calib_mul(csi_src)

% sto_calib_mul

% Input:

% - csi_src is the csi data with sto; [T S A L]

% Output:

% - csi_remove_sto is the csi data without sto; [T S A L]

antenna_num = size(csi_src, 3);

csi_remove_sto = zeros(size(csi_src));

for a = 1:antenna_num

a_nxt = mod(a, antenna_num) + 1;

csi_remove_sto(:, :, a, :) = csi_src(:, :, a, :) .* conj(csi_src(:, :, a_nxt, :));

end

end

⬇

function [csi_remove_sto] = sto_calib_div(csi_src)

% sto_calib_div

% Input:

% - csi_src is the csi data with sto; [T S A L]

% Output:

% - csi_remove_sto is the csi data without sto; [T S A L]

antenna_num = size(csi_src, 3);

csi_remove_sto = zeros(size(csi_src));

for a = 1:antenna_num

a_nxt = mod(a, antenna_num) + 1;

csi_remove_sto(:, :, a, :) = csi_src(:, :, a, :) ./ csi_src(:, :, a_nxt, :);

end

end

To sum up, there are various forms of errors in Wi-Fi CSI measurements, including fixed bias and random errors. Each of them have different impacts on the localization, tracking, and sensing tasks. The erroneous CSI form can be finally written as:

(33) $$\tilde{H}(i,j,k)=\sum _{n=1}^N{\beta }_{i}f({\alpha }_n)e^{-jg({\varphi }_n(i,j,k))}+N(i,j,k),{\tilde{\varphi }}_n(i,j,k)=2\pi (f_c+\mathrm{\Delta }{f}_j+f_D+{ϵ}_f)({\tau }_n(i,j,k)+{ϵ}_t)+{ϵ}_{\varphi }.$$

6. Wireless Sensing with Deep Learning

This section introduces a series of learning algorithms, especially the prevalent deep neural network models such as CNN and RNN, and their applications in wireless sensing. This section also proposes a complex-valued neural network to accomplish learning and inference based on wireless features efficiently.

6.1. Convolutional Neural Network

Convolutional Neural Network (CNN) contributes to the recent advances in understanding images, videos, and audios. Some works (Zhao et al., 2018; Zheng et al., 2019; Jiang et al., 2020) have exploited CNN for wireless signal understanding in wireless sensing tasks and achieved promising performance. This section will present a working example to demonstrate how to apply CNN for wireless sensing. Specifically, we use commodity Wi-Fi to recognize six human gestures. The gestures are illustrated in Fig. 9. We deploy a Wi-Fi transmitter and six receivers in a typical classroom, and the device setup is sketched in Fig. 10. The users are asked to perform gestures at the five marked locations and to five orientations. The data samples can be found in our released dataset (Yang et al., 2020). We extract DFS from raw CSI signals and feed them into a CNN network. The network architecture is shown in Fig. 11.

Figure 9. Sketches of gestures evaluated in the experiment.

Figure 10. The setup of WiFi devices for gesture recognition task.

Figure 11. Convolutional Neural Network architecture.

We now introduce the implementation code in detail.

First, some necessary packages are imported. We use Keras (Chollet et al., 2015) API with TensorFlow as the backend to demonstrate how to implement the neural network.

⬇

import os,sys

import numpy as np

import scipy.io as scio

import tensorflow as tf

import keras

from keras.layers import Input, GRU, Dense, Flatten, Dropout, Conv2D, Conv3D, MaxPooling2D, MaxPooling3D, TimeDistributed, Bidirectional, Multiply, Permute, RepeatVector, Concatenate, Dot, Lambda

from keras.models import Model, load_model

import keras.backend as K

from sklearn.metrics import confusion_matrix

from keras.backend.tensorflow_backend import set_session

from sklearn.model_selection import train_test_split

Then we define some parameters, including the hyperparameters and the data path. The fraction of testing data is defined as 0.1. To simplify the problem, we only use six gesture types in the widar3.0 dataset.

⬇

# Parameters

fraction_for_test = 0.1

data_dir = ’widar30dataset/DFS/20181130/’

ALL_MOTION = [1,2,3,4,5,6]

N_MOTION = len(ALL_MOTION)

T_MAX = 0

n_epochs = 200

f_dropout_ratio = 0.5

n_gru_hidden_units = 64

n_batch_size = 32

f_learning_rate = 0.001

The program begins with loading data with the predefined function load_data. The loaded data are split into train and test by calling the API function train_test_split. The labels of the training data are encoded into the one-hot format with the predefined function onehot_encoding.

⬇

# Load data

data, label = load_data(data_dir)

print(’\nLoaded dataset of ’ + str(label.shape[0]) + ’ samples, each sized ’ + str(data[0,:,:,:,:].shape) + ’\n’)

# Split train and test

[data_train, data_test, label_train, label_test] = train_test_split(data, label, test_size=fraction_for_test)

print(’\nTrain on ’ + str(label_train.shape[0]) + ’ samples\n’ +\

’Test on ’ + str(label_test.shape[0]) + ’ samples\n’)

# One-hot encoding for train data

label_train = onehot_encoding(label_train, N_MOTION)

After loading and formatting the training and testing data, we defined the model with the predefined function build_model. After that, we train the model by calling the API function fit. The input data and label are specified in the parameters. The fraction of validation data is specified as 0.1.

⬇

# Train Model

model = build_model(input_shape=(T_MAX, 6, 121, 1), n_class=N_MOTION)

model.summary()

model.fit({’name_model_input’: data_train},{’name_model_output’: label_train},

batch_size=n_batch_size,

epochs=n_epochs,

verbose=1,

validation_split=0.1, shuffle=True)

After the training process, we evaluate the model with the test dataset. The predictions are converted from one-hot format to integers and are used to calculate the confusion matrix and accuracy.

⬇

# Testing...

print(’Testing...’)

label_test_pred = model.predict(data_test)

label_test_pred = np.argmax(label_test_pred, axis = -1) + 1

# Confusion Matrix

cm = confusion_matrix(label_test, label_test_pred)

print(cm)

cm = cm.astype(’float’)/cm.sum(axis=1)[:, np.newaxis]

cm = np.around(cm, decimals=2)

print(cm)

# Accuracy

test_accuracy = np.sum(label_test == label_test_pred) / (label_test.shape[0])

print(test_accuracy)

The predefined onehot_encoding function convert the label to one-hot format.

⬇

def onehot_encoding(label, num_class):

# label(ndarray)=>_label(ndarray): [N,]=>[N,num_class]

label = np.array(label).astype(’int32’)

label = np.squeeze(label)

_label = np.eye(num_class)[label-1]

return _label

The predefined load_data function is used to load all data samples and labels from a directory. Each file in the directory corresponds to a single data sample. Each data sample is normalized with the predefined normalize_data function. It is worth noting that the data samples have different time durations. We use a predefined zero_padding function to make their durations the same as the longest one.

⬇

def load_data(path_to_data):

global T_MAX

data = []

label = []

for data_root, data_dirs, data_files in os.walk(path_to_data):

for data_file_name in data_files:

file_path = os.path.join(data_root,data_file_name)

try:

data_1 = scio.loadmat(file_path)[’doppler_spectrum’] # [6,121,T]

label_1 = int(data_file_name.split(’-’)[1])

location = int(data_file_name.split(’-’)[2])

orientation = int(data_file_name.split(’-’)[3])

repetition = int(data_file_name.split(’-’)[4])

# Downsample

data_1 = data_1[:,:,0::10]

# Select Motion

if (label_1 not in ALL_MOTION):

continue

# Normalization

data_normed_1 = normalize_data(data_1)

# Update T_MAX

if T_MAX < np.array(data_1).shape[2]:

T_MAX = np.array(data_1).shape[2]

except Exception:

continue

# Save List

data.append(data_normed_1.tolist())

label.append(label_1)

# Zero-padding

data = zero_padding(data, T_MAX)

# Swap axes

data = np.swapaxes(np.swapaxes(data, 1, 3), 2, 3) # [N,6,121,T_MAX]=>[N,T_MAX,6,121]

data = np.expand_dims(data, axis=-1) # [N,T_MAX,6,121]=>[N,T_MAX,6,121,1]

# Convert label to ndarray

label = np.array(label)

# data(ndarray): [N,T_MAX,6,121,1], label(ndarray): [N,]

return data, label

The normalize_data function is used to normalize the loaded data samples. Each data sample has a dimension of $[6,121,T]$, in which the number ”6” represents the number of Wi-Fi receivers, the number ”121” represents the frequency bins, and the ”T” represents the time durations. To normalize a sample, we scale the data to be in the range of $[0,1]$ for each time snapshot.

⬇

def normalize_data(data_1):

# data(ndarray)=>data_norm(ndarray): [6,121,T]=>[6,121,T]

data_1_max = np.amax(data_1,(0,1),keepdims=True) # [6,121,T]=>[1,1,T]

data_1_min = np.amin(data_1,(0,1),keepdims=True) # [6,121,T]=>[1,1,T]

data_1_max_rep = np.tile(data_1_max,(data_1.shape[0],data_1.shape[1],1)) # [1,1,T]=>[6,121,T]

data_1_min_rep = np.tile(data_1_min,(data_1.shape[0],data_1.shape[1],1)) # [1,1,T]=>[6,121,T]

data_1_norm = (data_1 - data_1_min_rep) / (data_1_max_rep - data_1_min_rep + sys.float_info.min)

return data_1_norm

The zero_padding function is used to align all the data samples to have the same duration. The padded length is specified by the parameter T_MAX.

⬇

def zero_padding(data, T_MAX):

# data(list)=>data_pad(ndarray): [6,121,T1/T2/...]=>[6,121,T_MAX]

data_pad = []

for i in range(len(data)):

t = np.array(data[i]).shape[2]

data_pad.append(np.pad(data[i], ((0,0),(0,0),(T_MAX - t,0)), ’constant’, constant_values = 0).tolist())

return np.array(data_pad)

In this function, we define the network structure. The input layer is specified with the API function Input, which has the parameters to define the input shape, data type, and the layer name. Following the input layer, we use a three-dimensional convolutional layer, a max-pooling layer, and two fully connected layers. The output layer is specified with the API function Output, which has the parameters to define the activation function, the dimension, and the name. At last, we finalize the model with the API function Model and compile. The optimizer is specified to RMSprop, and the loss is specified to categorical_crossentropy.

⬇

def build_model(input_shape, n_class):

model_input = Input(shape=input_shape, dtype=’float32’, name=’name_model_input’) # (@,T_MAX,6,121,1)

# CNN

x = Conv3D(16,kernel_size=(5,3,6),activation=’relu’,data_format=’channels_last’,\

input_shape=input_shape)(model_input) # (@,T_MAX-4,6,121,1)=>(@,T_MAX-4,4,116,16)

x = MaxPooling3D(pool_size=(2,2,2))(x) # (@,T_MAX-4,4,116,16)=>(@,T_MAX-4,2,58,16)

x = Flatten()(x) # (@,T_MAX-4,2,58,16)=>(@,(T_MAX-4)*2*58*16)

x = Dense(256,activation=’relu’)(x) # (@,(T_MAX-4)*2*58*16)=>(@,256)

x = Dropout(f_dropout_ratio)(x)

x = Dense(128,activation=’relu’)(x) # (@,256)=>(@,128)

x = Dropout(f_dropout_ratio)(x)

model_output = Dense(n_class, activation=’softmax’, name=’name_model_output’)(x) # (@,128)=>(@,n_class)

# Model compiling

model = Model(inputs=model_input, outputs=model_output)

model.compile(optimizer=keras.optimizers.RMSprop(lr=f_learning_rate),

loss=’categorical_crossentropy’,

metrics=[’accuracy’]

)

return model

6.2. Recurrent Neural Network

Recurrent Neural Network (RNN) is designed for modeling temporal dynamics of sequences and is commonly used for time series data analysis like speech recognition and natural language processing. Wireless signals are highly correlated over time and can be processed with RNN. Some works (Zheng et al., 2019; Jiang et al., 2020) have demonstrated the potential of RNN for wireless sensing tasks. In this section, we will present a working example of combining CNN and RNN to perform gesture recognition with Wi-Fi. The experimental settings are the same as in Sec. 6.1. We also extract DFS from the raw CSI as the input feature of the network. The network architecture is shown in Fig. 12.

We now introduce the implementation code in detail.

Figure 12. Recurrent Neural Network architecture.

Most of the code is the same as in Sec. 6.1 except for the model definition. To define the model, we use two-dimensional convolutional layer and max-pooling layer on the dimensional except for the time dimension of the data. We adopt the GRU layer as the recurrent layer.

⬇

def build_model(input_shape, n_class):

model_input = Input(shape=input_shape, dtype=’float32’, name=’name_model_input’) # (@,T_MAX,6,121,1)

# CNN+RNN

x = TimeDistributed(Conv2D(16,kernel_size=(3,6),activation=’relu’,data_format=’channels_last’,\

input_shape=input_shape))(model_input) # (@,T_MAX,6,121,1)=>(@,T_MAX,4,116,16)

x = TimeDistributed(MaxPooling2D(pool_size=(2,2)))(x) # (@,T_MAX,4,116,16)=>(@,T_MAX,2,58,16)

x = TimeDistributed(Flatten())(x) # (@,T_MAX,2,58,16)=>(@,T_MAX,2*58*16)

x = TimeDistributed(Dense(128,activation=’relu’))(x) # (@,T_MAX,2*58*16)=>(@,T_MAX,128)

x = TimeDistributed(Dropout(f_dropout_ratio))(x)

x = TimeDistributed(Dense(64,activation=’relu’))(x) # (@,T_MAX,128)=>(@,T_MAX,64)

x = GRU(n_gru_hidden_units,return_sequences=False)(x) # (@,T_MAX,64)=>(@,64)

x = Dropout(f_dropout_ratio)(x)

model_output = Dense(n_class, activation=’softmax’, name=’name_model_output’)(x) # (@,64)=>(@,n_class)

# Model compiling

model = Model(inputs=model_input, outputs=model_output)

model.compile(optimizer=keras.optimizers.RMSprop(lr=f_learning_rate),

loss=’categorical_crossentropy’,

metrics=[’accuracy’]

)

return model

6.3. Adversarial Learning

Except for the basic neural network components, some high-level network architectures also play essential roles in wireless sensing. Similar to computer vision tasks, wireless sensing also suffer from domain misalignment problem. Wireless signals can be reflected by the surrounding objects during propagation and will be flooded with target-irrelevant signal components. The sensing system trained in one deployment environment can hardly be applied directly in other settings without adaptation. Some works (Jiang et al., 2018) try to adopt adversarial learning techniques to tackle this problem and achieve promising performance. This section will give an example of how to apply this technique in wireless sensing tasks. Specifically, we build a gesture recognition system with Wi-Fi, similar to that in Sec. 6.1. We try to achieve consistent performance across different human locations and orientations. The network architecture is shown in Fig. 13.

We now introduce the implementation code in detail.

Figure 13. Adversarial learning network architecture.

In the Widar3.0 dataset (Yang et al., 2020), we collect gesture data when the users stand at different locations. As discussed in Sec. 4.3, human locations have significant impact on the DFS measurements. To mitigate this impact, we treat human locations as different domains and build an adversarial learning network to recognize gestures irrespective of domains. In the program, we first load data, labels, and domains from the dataset and split them into train and test. Both label and domain are encoded into the one-hot format.

⬇

# Load data

data, label, domain = load_data(data_dir)

print(’\nLoaded dataset of ’ + str(label.shape[0]) + ’ samples, each sized ’ + str(data[0,:,:,:,:].shape) + ’\n’)

# Split train and test

[data_train, data_test, label_train, label_test, domain_train, domain_test] = train_test_split(data, label, domain, test_size=fraction_for_test)

print(’\nTrain on ’ + str(label_train.shape[0]) + ’ samples\n’ +\

’Test on ’ + str(label_test.shape[0]) + ’ samples\n’)

# One-hot encoding for train data

label_train = onehot_encoding(label_train, N_MOTION)

domain_train = onehot_encoding(domain_train, N_LOCATION)

After loading and formating data, we built the network and trained it from scratch. The training data, label, and domain are passed to the API function fit for training.

⬇

# Train Model

model = build_model(input_shape=(T_MAX, 6, 121, 1), n_class=N_MOTION, n_domain=N_LOCATION)

model.summary()

model.fit({’name_model_input’: data_train},{’name_model_output_label’: label_train, ’name_model_output_domain’: domain_train},

batch_size=n_batch_size,

epochs=n_epochs,

verbose=1,

validation_split=0.1, shuffle=True)

After the training process finishes, we evaluate the network with the test samples. Note that the adversarial network has both label and domain prediction outputs. We only use the label output for accuracy evaluation.

⬇

# Testing...

print(’Testing...’)

[label_test_pred,_] = model.predict(data_test)

label_test_pred = np.argmax(label_test_pred, axis = -1) + 1

# Confusion Matrix

cm = confusion_matrix(label_test, label_test_pred)

print(cm)

cm = cm.astype(’float’)/cm.sum(axis=1)[:, np.newaxis]

cm = np.around(cm, decimals=2)

print(cm)

# Accuracy

test_accuracy = np.sum(label_test == label_test_pred) / (label_test.shape[0])

print(test_accuracy)

Different from that in Sec. 6.1, we load data, label, and domain in the load_data function. The domain is defined as the location of the human, which is embedded in the file name.

⬇

def load_data(path_to_data):

global T_MAX

data = []

label = []

domain = []

for data_root, data_dirs, data_files in os.walk(path_to_data):

for data_file_name in data_files:

file_path = os.path.join(data_root,data_file_name)

try:

data_1 = scio.loadmat(file_path)[’doppler_spectrum’] # [6,121,T]

label_1 = int(data_file_name.split(’-’)[1])

location_1 = int(data_file_name.split(’-’)[2])

orientation_1 = int(data_file_name.split(’-’)[3])

repetition_1 = int(data_file_name.split(’-’)[4])

# Downsample

data_1 = data_1[:,:,0::10]

# Select Motion

if (label_1 not in ALL_MOTION):

continue

# Normalization

data_normed_1 = normalize_data(data_1)

# Update T_MAX

if T_MAX < np.array(data_1).shape[2]:

T_MAX = np.array(data_1).shape[2]

except Exception:

continue

# Save List

data.append(data_normed_1.tolist())

label.append(label_1)

domain.append(location_1)

# Zero-padding

data = zero_padding(data, T_MAX)

# Swap axes

data = np.swapaxes(np.swapaxes(data, 1, 3), 2, 3) # [N,6,121,T_MAX]=>[N,T_MAX,6,121]

data = np.expand_dims(data, axis=-1) # [N,T_MAX,6,121]=>[N,T_MAX,6,121,1]

# Convert label and domain to ndarray

label = np.array(label)

domain = np.array(domain)

# data(ndarray): [N,T_MAX,6,121,1], label(ndarray): [N,], domain(ndarray): [N,]

return data, label, domain

To define the network, we use a CNN layer and an RNN layer as the feature extractor, which is similar to that in Sec. 6.2. In the gesture recognizer and domain discriminator, we use two fully-connected layers and an output layer activated by softmax function, respectively. We use categorical cross-entropy loss for both label prediction and domain prediction outputs. The domain prediction loss is weighted with loss_weight_domain and subtracted from the label prediction loss.

⬇

def build_model(input_shape, n_class, n_domain):

model_input = Input(shape=input_shape, dtype=’float32’, name=’name_model_input’) # (@,T_MAX,6,121,1)

# CNN+RNN+Adversarial

x = TimeDistributed(Conv2D(16,kernel_size=(3,6),activation=’relu’,data_format=’channels_last’,\

input_shape=input_shape))(model_input) # (@,T_MAX,6,121,1)=>(@,T_MAX,4,116,16)

x = TimeDistributed(MaxPooling2D(pool_size=(2,2)))(x) # (@,T_MAX,4,116,16)=>(@,T_MAX,2,58,16)

x = TimeDistributed(Flatten())(x) # (@,T_MAX,2,58,16)=>(@,T_MAX,2*58*16)

x = TimeDistributed(Dense(128,activation=’relu’))(x) # (@,T_MAX,2*58*16)=>(@,T_MAX,128)

x = TimeDistributed(Dropout(f_dropout_ratio))(x)

x = TimeDistributed(Dense(64,activation=’relu’))(x) # (@,T_MAX,128)=>(@,T_MAX,64)

x = GRU(n_gru_hidden_units,return_sequences=False)(x) # (@,T_MAX,64)=>(@,64)

x_feat = Dropout(f_dropout_ratio)(x)

# Label prediction part

x_1 = Dense(64, activation=’relu’)(x_feat) # (@,64)=>(@,64)

x_1 = Dense(32, activation=’relu’)(x_1) # (@,64)=>(@,32)

model_output_label = Dense(n_class, activation=’softmax’, name=’name_model_output_label’)(x_1) # (@,32)=>(@,n_class)

# Domain prediction part

x_2 = Dense(64, activation=’relu’)(x_feat) # (@,64)=>(@,64)

x_2 = Dense(32, activation=’relu’)(x_2) # (@,64)=>(@,32)

model_output_domain = Dense(n_domain, activation=’softmax’, name=’name_model_output_domain’)(x_2) # (@,32)=>(@,n_domain)

model = Model(inputs=model_input, outputs=[model_output_label, model_output_domain])

model.compile(optimizer=keras.optimizers.RMSprop(lr=f_learning_rate),

loss = {’name_model_output_label’:custom_loss_label(), ’name_model_output_domain’:custom_loss_domain()},

loss_weights={’name_model_output_label’:1, ’name_model_output_domain’:-1*loss_weight_domain},

metrics={’name_model_output_label’:’accuracy’, ’name_model_output_domain’:’accuracy’}

)

return model

The pre-defined custom_loss_label and custom_loss_domain are categorical crossentropy losses for both label prediction and domain prediction.

⬇

def custom_loss_label():

def lossfn(y_true, y_pred):

myloss_batch = -1 * K.sum(y_true*K.log(y_pred+K.epsilon()), axis=-1, keepdims=False)

myloss = K.mean(myloss_batch, axis=-1, keepdims=False)

return myloss

return lossfn

⬇

def custom_loss_domain():

def lossfn(y_true, y_pred):

myloss_batch = -1 * K.sum(y_true*K.log(y_pred+K.epsilon()), axis=-1, keepdims=False)

myloss = K.mean(myloss_batch, axis=-1, keepdims=False)

return myloss

return lossfn

6.4. Complex-valued Neural Network

In this section, we will present a more complicated wireless sensing task with deep learning. Many wireless sensing approaches employ Fast Fourier Transform (FFT) on a time series of RF data to obtain time-frequency spectrograms of human activities. FFT suffers from errors due to an effect known as leakage, when the block of data is not periodic (the most common case in practice), which results in a smeared spectrum of the original signal and further leads to misleading data representation for learning-based sensing. Classical approaches reduce leakage by windowing, which cannot eliminate leakage entirely. Considering the significant fitting capability of deep neural networks, we can design a signal processing network to learn an optimal function to minimize or nearly eliminate the leakage and enhance the spectrums, which we call the Signal Enhancement Network (SEN).

The signal processing network takes as input a spectrogram transformed from wireless signals via STFT, removes the spectral leakage in the spectrogram, and recovers the underlying actual frequency components. Fig. 16 shows the training process of the network. As shown in the upper part of Fig. 16, we randomly generate ideal spectrums with 1 to 5 frequency components, whose amplitudes, phases, and frequencies are uniformly drawn from their ranges of interest. Then, the ideal spectrum is converted to the leaked spectrum following the process in the following equation to simulate the windowing effect and random complex noises:

(34) $$\widehat{𝐬}=\mathrm{𝐀𝐬}+𝐧,$$

where $𝐬$ and $\widehat{𝐬}$ are the ideal and estimated frequency spectrum, respectively, $𝐧$ represents the additive Gaussian noise vector, and $𝐀$ is the convolution matrix of the windowing function in the frequency domain. The $i^{th}$ column of $𝐀$ is:

(35) $${𝐀}_{(:,𝐢)}=\mathrm{FFT}(\varpi )\ast \delta (i).$$

where $\varpi$ represents the windowing function of FFT in time domain.

The amplitude of the noise follows a Gaussian distribution and its phase follows a uniform distribution in $[0,2\pi ]$. The network takes the leaked spectrum as input and outputs the enhanced spectrum close to the ideal one. Thus, we minimize the $L_2$ loss $L={\Vert \mathrm{SEN}(\widehat{𝐬})-𝐬\Vert }_2$ during training. During inference, the spectrums measured from real-world scenarios are normalized to $[0,1]$ and fed into the network to obtain the enhanced spectrum.

We now present the implementation code in detail.

(a)

(b)

Figure 14. Comparison between (a) real-valued and (b) complex-valued neurons.

Figure 15. Complex-valued neural network architecture.

Figure 16. The training process of the signal processing network.

Different from previous sections, we use the PyTorch platform to implement the network. PyTorch provides the interface to implement custom layers, which makes the implementation of the CVNN much more convenient. We first import some necessary packages as follows.

⬇

import os,sys,math,scipy,imp

import numpy as np

import scipy.io as scio

import torch, torchvision

import torch.nn as nn

from scipy import signal

from math import sqrt,log,pi

from torch.fft import fft,ifft

from torch.nn.functional import relu, softmax, cross_entropy

from torch import sigmoid,tanh

from torch.nn import MSELoss as MSE

Some parameters are defined in this part.

⬇

# Definition

wind_len = 125

wind_type = ’gaussian’

n_max_freq_component = 3

AWGN_amp = 0.01

str_modelname_prefix = ’./SEN_Results/SEN_’ + wind_type + ’_W’ + str(wind_len)

str_model_name_pretrained = str_modelname_prefix + ’_E1000.pt’

feature_len = 121

padded_len = 1000

crop_len = feature_len

blur_matrix_left = []

# Hyperparameters

n_begin_epoch = 1

n_epoch = 10000

n_itr_per_epoch = 500

n_batch_size = 64

n_test_size = 200

f_learning_rate = 0.001

The program begins with the following code. We first generate the convolution matrix of the windowing function (Eqn. 35) with the pre-defined function generate_blur_matrix_complex, which will be introduced shortly. Then we define the SEN model and move it to the GPU processor with the API function cuda. After the model definition, we train the model with synthetic spectrograms, during which process we save the trained model every 500 epochs. The trained and saved model can be directly loaded and used to enhance spectrograms.

⬇

if __name__ == "__main__":

# Generate blur matrix

blur_matrix_right = generate_blur_matrix_complex(wind_type=wind_type, wind_len=wind_len, padded_len=padded_len, crop_len=crop_len)

# Define model

print(’Model building...’)

model = SEN(feature_len=feature_len)

model.cuda()

# Train model

print(’Model training...’)

train(model=model, blur_matrix_right=blur_matrix_right, feature_len=feature_len, n_epoch=n_epoch, n_itr_per_epoch=n_itr_per_epoch, n_batch_size=n_batch_size, optimizer=torch.optim.RMSprop(model.parameters(), lr=f_learning_rate))

This generate_blur_matrix_complex function is used to generate the convolution matrix of the windowing function. The core idea behind this function is to enumerate all the frequencies, apply the window function on the sinusoid signal, and generate the corresponding spectrums with FFT. After obtaining the convolution matrix, we can bridge the gap between the ideal and the leaked spectrograms with Eqn. 34. In other words, we can directly get the leaked spectrograms by multiplying the idea spectrograms with the convolution matrix.

⬇

def generate_blur_matrix_complex(wind_type, wind_len=251, padded_len=1000, crop_len=121):

# Generate matrix used to introduce spec leakage in complex domain

# ret: (ndarray.complex128) [crop_len, crop_len](row first)

# Row first: each row represents the spectrum of one single carrier

# Steps: carrier/windowing/pad/fft/crop/unwrap/norm

# Parameters offloading

fs = 1000

n_f_bins = crop_len

f_high = int(n_f_bins/2)

f_low = -1 * f_high

init_phase = 0

# Carrier

t_ = np.arange(0,wind_len).reshape(1,wind_len)/fs # [1,wind_len] (0~wind_len/fs seconds)

freq = np.arange(f_low,f_high+1,1).reshape(n_f_bins,1) # [n_f_bins,1] (f_low~f_high Hz)

phase = 2 * pi * freq * t_ + init_phase # [n_f_bins,wind_len]

signal = np.exp(1j*phase) # [n_f_bins,wind_len]~[121,251]

# Windowing

if wind_type == ’gaussian’:

window = scipy.signal.windows.gaussian(wind_len, (wind_len-1)/sqrt(8*log(200)), sym=True) # [wind_len,]

else:

window = scipy.signal.get_window(wind_type, wind_len)

sig_wind = signal * window # [n_f_bins,wind_len]*[wind_len,]=[n_f_bins,wind_len]~[121,251]

# Pad/FFT

sig_wind_pad = np.concatenate((sig_wind, np.zeros((n_f_bins,padded_len-wind_len))),axis=1) # [n_f_bins,wind_len]=>[n_f_bins,padded_len]

sig_wind_pad_fft = np.fft.fft(sig_wind_pad, axis=-1) # [n_f_bins,padded_len]~[121,1000]

# Crop

n_freq_pos = f_high + 1

n_freq_neg = abs(f_low)

sig_wind_pad_fft_crop = np.concatenate((sig_wind_pad_fft[:,:n_freq_pos],\

sig_wind_pad_fft[:,-1*n_freq_neg:]), axis=1) # [n_f_bins,crop_len]~[121,121]

# Unwrap

n_shift = n_freq_neg

sig_wind_pad_fft_crop_unwrap = np.roll(sig_wind_pad_fft_crop, shift=n_shift, axis=1) # [n_f_bins,crop_len]~[121,121]

# Norm (amp_max=1)

_sig_amp = np.abs(sig_wind_pad_fft_crop_unwrap)

_sig_ang = np.angle(sig_wind_pad_fft_crop_unwrap)

_max = np.tile(_sig_amp.max(axis=1,keepdims=True), (1,crop_len))

_min = np.tile(_sig_amp.min(axis=1,keepdims=True), (1,crop_len))

_sig_amp_norm = _sig_amp / _max

sig_wind_pad_fft_crop_unwrap_norm = _sig_amp_norm * np.exp(1j*_sig_ang)

# Return

ret = sig_wind_pad_fft_crop_unwrap_norm

return ret

This function is to generate one batch of spectrograms in both the leaked form and the idea form. The process is the same as that illustrated in Eqn. 34.

⬇

def syn_one_batch_complex(blur_matrix_right, feature_len, n_batch):

# Syn. HiFi, blurred and AWGN signal in complex domain

# ret: (ndarray.complex128) [@,feature_len]

# blur_matrix_right: Row first (each row represents the spectrum of one single carrier)

# Syn. x [@,feature_len]

x = np.zeros((n_batch, feature_len))*np.exp(1j*0)

for i in range(n_batch):

num_carrier = int(np.random.randint(0,n_max_freq_component,1))

idx_carrier = np.random.permutation(feature_len)[:num_carrier]

x[i,idx_carrier] = np.random.rand(1,num_carrier) * np.exp(1j*( 2*pi*np.random.rand(1,num_carrier) - pi ))

# Syn. x_blur [@,feature_len]

x_blur = x @ blur_matrix_right

# Syn. x_tilde [@,feature_len]

x_tilde = x_blur + 2*AWGN_amp*(np.random.random(x_blur.shape)-0.5) *\

np.exp(1j*( 2*pi*np.random.random(x_blur.shape) - pi ))

return x, x_blur, x_tilde

This part demonstrates the code on how to implement the SEN network. According to the API of PyTorch, both the __init__ and the forward interfaces should be implemented with customized algorithms. In the __init__ function, we defined five complex-valued fully-connected layers. In the forward function, we defined the network structure by concatenating the five FC layers and specifying the input and output layers.

⬇

class SEN(nn.Module):

def __init__(self, feature_len):

super(SEN, self).__init__()

self.feature_len = feature_len

self.fc_1 = m_Linear(feature_len, feature_len)

self.fc_2 = m_Linear(feature_len, feature_len)

self.fc_3 = m_Linear(feature_len, feature_len)

self.fc_4 = m_Linear(feature_len, feature_len)

self.fc_out = m_Linear(feature_len, feature_len)

def forward(self, x):

h = x # (@,*,2,H)

h = tanh(self.fc_1(h)) # (@,*,2,H)=>(@,*,2,H)

h = tanh(self.fc_2(h)) # (@,*,2,H)=>(@,*,2,H)

h = tanh(self.fc_3(h)) # (@,*,2,H)=>(@,*,2,H)

h = tanh(self.fc_4(h)) # (@,*,2,H)=>(@,*,2,H)

output = tanh(self.fc_out(h)) # (@,*,2,H)=>(@,*,2,H)

return output

This m_Linear class leverages the interface of PyTorch to define the customized complex-valued fully-connected layer, which is the implementation of the network structure illustrated in Fig. 14.

⬇

class m_Linear(nn.Module):

def __init__(self, size_in, size_out):

super().__init__()

self.size_in, self.size_out = size_in, size_out

# Creation

self.weights_real = nn.Parameter(torch.randn(size_in, size_out, dtype=torch.float32))

self.weights_imag = nn.Parameter(torch.randn(size_in, size_out, dtype=torch.float32))

self.bias = nn.Parameter(torch.randn(2, size_out, dtype=torch.float32))

# Initialization

nn.init.xavier_uniform_(self.weights_real, gain=1)

nn.init.xavier_uniform_(self.weights_imag, gain=1)

nn.init.zeros_(self.bias)

def swap_real_imag(self, x):

# [@,*,2,Hout]

# [real, imag] => [-1*imag, real]

h = x # [@,*,2,Hout]

h = h.flip(dims=[-2]) # [@,*,2,Hout] [real, imag]=>[imag, real]

h = h.transpose(-2,-1) # [@,*,Hout,2]

h = h * torch.tensor([-1,1]).cuda() # [@,*,Hout,2] [imag, real]=>[-1*imag, real]

h = h.transpose(-2,-1) # [@,*,2,Hout]

return h

def forward(self, x):

# x: [@,*,2,Hin]

h = x # [@,*,2,Hin]

h1 = torch.matmul(h, self.weights_real) # [@,*,2,Hout]

h2 = torch.matmul(h, self.weights_imag) # [@,*,2,Hout]

h2 = self.swap_real_imag(h2) # [@,*,2,Hout]

h = h1 + h2 # [@,*,2,Hout]

h = torch.add(h, self.bias) # [@,*,2,Hout]+[2,Hout]=>[@,*,2,Hout]

return h

This loss_function defines the loss function of the SEN network. The loss is defined as the Euclidean distance between the idea spectrum and the netwrok predicted spectrum. Only the amplitude of the spectrums are considered.

⬇

def loss_function(x, y):

# x,y: [@,*,2,H]

x = torch.linalg.norm(x,dim=-2) # [@,*,2,H]=>[@,*,H]

y = torch.linalg.norm(y,dim=-2) # [@,*,2,H]=>[@,*,H]

# MSE loss for Amp

loss_recon = MSE(reduction=’mean’)(x, y)

return loss_recon

In this train function, we implement the training process of SEN as that described in Fig. 16. In each epoch, we generate and train the network with multiple iterations. For each iteration, we generate a batch of synthetic spectrums in leaked format. Each leaked spectrum have an idea spectrum as the label.

⬇

def train(model, blur_matrix_right, feature_len, n_epoch, n_itr_per_epoch, n_batch_size, optimizer):

for i_epoch in range(n_begin_epoch, n_epoch+1):

model.train()

total_loss_this_epoch = 0

for i_itr in range(n_itr_per_epoch):

x, _, x_tilde = syn_one_batch_complex(blur_matrix_right=blur_matrix_right, feature_len=feature_len, n_batch=n_batch_size)

x = complex_array_to_bichannel_float_tensor(x)

x_tilde = complex_array_to_bichannel_float_tensor(x_tilde)

x = x.cuda()

x_tilde = x_tilde.cuda()

optimizer.zero_grad()

y = model(x_tilde)

loss = loss_function(x, y)

loss.backward()

optimizer.step()

total_loss_this_epoch += loss.item()

if i_itr % 10 == 0:

print(’--------> Epoch: {}/{} loss: {:.4f} [itr: {}/{}]’.format(

i_epoch+1, n_epoch, loss.item() / n_batch_size, i_itr+1, n_itr_per_epoch), end=’\r’)

# Validate

model.eval()

x, _, x_tilde = syn_one_batch_complex(blur_matrix_right=blur_matrix_right, feature_len=feature_len, n_batch=n_batch_size)

x = complex_array_to_bichannel_float_tensor(x)

x_tilde = complex_array_to_bichannel_float_tensor(x_tilde)

x = x.cuda()

x_tilde = x_tilde.cuda()

y = model(x_tilde)

total_valid_loss = loss_function(x, y)

print(’========> Epoch: {}/{} Loss: {:.4f}’.format(i_epoch+1, n_epoch, total_valid_loss) + ’ ’ + wind_type + ’_’ + str(wind_len) + ’ ’*20)

if i_epoch % 500 == 0:

torch.save(model, str_modelname_prefix+’_E’+str(i_epoch)+’.pt’)

This function converts the complex-valued arrays to double channel real-valued arrays. This is because the GPU only supports the calculations of real numbers. We use a little trick to implement the complex-valued network by separating the real and imaginary parts into two real arrays.

⬇

def complex_array_to_bichannel_float_tensor(x):

# x: (ndarray.complex128) [@,*,H]

# ret: (tensor.float32) [@,*,2,H]

x = x.astype(’complex64’)

x_real = x.real # [@,*,H]

x_imag = x.imag # [@,*,H]

ret = np.stack((x_real,x_imag), axis=-2) # [@,*,H]=>[@,*,2,H]

ret = torch.tensor(ret)

return ret

This function converts the double channel real-valued arrays into complex-valued arrays.

⬇

def bichannel_float_tensor_to_complex_array(x):

# x: (tensor.float32) [@,*,2,H]

# ret: (ndarray.complex64) [@,*,H]

x = x.numpy()

x = np.moveaxis(x,-2,0) # [@,*,2,H]=>[2,@,*,H]

x_real = x[0,:]

x_imag = x[1,:]

ret = x_real + 1j*x_imag

return ret

After training the SEN network with sufficient epochs, we test the performance with spectrograms collected from Wi-Fi.

First, we define some parameters. The STFT window width is set to 125, and the window type is set to “gaussian”. The path to the pre-trained model and the CSI file is selected.

⬇

W = 125

wind_type = ’gaussian’

str_model_name = ’./SEN_Results/SEN_’ + wind_type + ’_W’ + str(W) + ’_E500.pt’

file_path_csi = ’Widar3_data/20181130/user1/user1-1-1-1-1.mat’

The program begins with loading the pre-trained model. Then, the CSI data is loaded and transformed to spectrograms with the predefined function csi_to_spec. After that, the complex-valued spectrograms are transformed to double real-valued channel tensors and processed with the SEN model. The results and the raw spectrograms are stored in files.

⬇

if __name__ == "__main__":

# Load trained model

print(’Loading model...’)

model = torch.load(str_model_name)

print(’Testing model...’)

model.eval()

with torch.no_grad():

# Import raw spectrogram

data_1 = csi_to_spec()

# Enhance spectrogram

x_tilde = complex_array_to_bichannel_float_tensor(data_1) # [6,121,T]=>[6,121,2,T]

x_tilde = x_tilde.permute(0,3,2,1) # [6,121,2,T]=>[6,T,2,121]

y = model(x_tilde.cuda()).cpu() # [6,T,2,121]

y = bichannel_float_tensor_to_complex_array(y) # [6,T,121]

y = np.transpose(y,(0,2,1)) # [6,T,121]=>[6,121,T]

scio.savemat(’SEN_test_x_tilde_complex_W’ + str(W) + ’.mat’, {’x_tilde’:data_1})

scio.savemat(’SEN_test_y_complex_W’ + str(W) + ’.mat’, {’y’:y})

This function first transforms the CSI data into spectrograms and crops the concerned frequency range between $[-60,60]$ Hz. Then, it unwraps the frequency bins and performs normalizations.

⬇

def csi_to_spec():

global file_path_csi

global W

signal = scio.loadmat(file_path_csi)[’csi_mat’].transpose() # [6,T] complex

# STFT

_, spec = STFT(signal, fs=1000, stride=1, wind_wid=W, dft_wid=1000, window_type=’gaussian’) # [6,1000,T]j

# Crop

spec_crop = np.concatenate((spec[:,:61], spec[:,-60:]), axis=1) # [1,1000,T]j=>[1,121,T]j

# Unwrap

spec_crop_unwrap = np.roll(spec_crop, shift=60, axis=1) # [1,121,T]j

# Normalize

spec_crop_unwrap_norm = normalize_data(spec_crop_unwrap) # [6,121,T] complex

if np.sum(np.isnan(spec_crop_unwrap_norm)):

print(’>>>>>>>>> NaN detected!’)

ret = spec_crop_unwrap_norm

return ret

This function transforms the time domain CSI data into frequency domain spectrograms with the API function scipy.signal.stft.

⬇

def STFT(signal, fs=1, stride=1, wind_wid=5, dft_wid=5, window_type=’gaussian’):

assert dft_wid >= wind_wid and wind_wid > 0 and stride <= wind_wid and stride > 0\

and isinstance(stride, int) and isinstance(wind_wid, int) and isinstance(dft_wid, int)\

and isinstance(fs, int) and fs > 0

if window_type == ’gaussian’:

window = scipy.signal.windows.gaussian(wind_wid, (wind_wid-1)/sqrt(8*log(200)), sym=True)

elif window_type == ’rect’:

window = np.ones((wind_wid,))

else:

window = scipy.signal.get_window(window_type, wind_wid)

f_bins, t_bins, stft_spectrum = scipy.signal.stft(x=signal, fs=fs, window=window, nperseg=wind_wid, noverlap=wind_wid-stride, nfft=dft_wid,\

axis=-1, detrend=False, return_onesided=False, boundary=’zeros’, padded=True)

return f_bins, stft_spectrum

This function scales the spectrograms to normalize the values into $[0,1]$.

⬇

def normalize_data(data_1):

# max=1

# data(ndarray.complex)=>data_norm(ndarray.complex): [6,121,T]=>[6,121,T]

data_1_abs = abs(data_1)

data_1_max = data_1_abs.max(axis=(1,2),keepdims=True) # [6,121,T]=>[6,1,1]

data_1_max_rep = np.tile(data_1_max,(1,data_1_abs.shape[1],data_1_abs.shape[2])) # [6,1,1]=>[6,121,T]

data_1_norm = data_1 / data_1_max_rep

return data_1_norm

Fig. 17 demonstrates the raw and enhanced spectrograms of pushing and pulling gestures.

(a)

(b)

Figure 17. Illustration of the spectrogram of a pushing and pulling gesture. (a) The measured spectrogram and (b) the enhanced spectrogram from the SEN.