基于WiFi的CSI数据做呼吸频率检测-python版(含代码和数据)

基于WiFi做研究的萌新可以系统学习大神总结的如下内容：
https://github.com/Mar s rocky/Awe some-WiFi-CSI-Sens in g

一、概述

本Demo无需机器学习模型，Demo 功能涉及的理论主要参考了硕士学位论文《基于WiFi的人体行为感知技术研究》，作者是南京邮电大学的朱XX，本人用python 复现了论文中呼吸频率检测的功能。Demo实现呼吸速率检测的主要过程为：

采数用的是C代码

1、通过 shell 脚本循环执行C代码进行csi 数据采集，形成一个个30秒的csi 数据文件（.d a t数据）；

解析和分析数据用python 代码

2、读取最新的.d a t数据文件，解析出csi数据；
3、计算 csi的振幅和相位，并对相位数据进行校准；
4、对振幅和相位数据进行中值滤波；
5、基于EMD 算法滤波；
6、基于FFT进行子载波筛选；
7、基于CA-CFAR 寻峰算法进行寻峰和呼吸速率统计；

二、操作 内容

1、配置好采数设备和代码运行环境，参考本人记录：
https://blog.csdn.net/Ac ecai01/a rt icle/d etails/129442761

2、布设试验场景：
在这里插入图片描述
3、选择一台发射数据的设备，输入如下发数据的命令：

xxx:~$: cd ~
xxx:~$: rfkill unblock wifi
xxx:~$: iwconfig
xxx:~$: sudo bash ./inject.sh wlan0 64 HT20
xxx:~$: echo 0x1c113 | sudo tee `sudo find /sys -name monitor_tx_rate`
xxx:~$: cd linux-80211n-csitool-supplementary/injection/
xxx:xxx$: sudo ./random_packets 1000000000 100 1 10000

以上命令的含义，参考本大节第1步骤的配置记录博客。
此时设备会按每秒100个数据帧的速率持续发送数据，以上命令设置的发送数据量够发115天，要中断发送直接按ctrl+c 即可。

4、选择另一台接收数据的设备，将本人修改的采数C代码 log_to_file.c 替换掉原先的log_to_file.c，先看修改后的log_to_file.c：

/*
 * (c) 2008-2011 Daniel Halperin <dhalperi@cs.washington.edu&gt;
 */
#include "iwl_connector.h"

#include <stdio.h&gt;
#include <stdlib.h&gt;
#include <string.h&gt;
#include <signal.h>
#include <unistd.h>
#include <arpa/inet.h>
#include <sys/socket.h>
#include <linux/netlink.h>

#define MAX_PAYLOAD 2048
#define SLOW_MSG_CNT 100

int sock_fd = -1; // the socket
FILE *out = NULL;

void check_usage(int argc, char **argv);

FILE *open_file(char *filename, char *spec);

void caught_signal(int sig);

void exit_program(int code);
void exit_program_err(int code, char *func);
void exit_program_with_alarm(int sig);

int main(int argc, char **argv)
{
	/* Local variables */
	struct sockaddr_nl proc_addr, kern_addr; // addrs for recv, send, bind
	struct cn_msg *cmsg;
	char buf[4096];
	int ret;
	unsigned short l, l2;
	int count = 0;

	/* Initialize signal*/
	signal(SIGALRM, exit_program_with_alarm);

	/* Make sure usage is correct */
	check_usage(argc, argv);

	/* Open and check log file */
	out = open_file(argv[1], "w");

	/* Setup the socket */
	sock_fd = socket(PF_NETLINK, SOCK_DGRAM, NETLINK_CONNECTOR);
	if (sock_fd == -1)
		exit_program_err(-1, "socket");

	/* Initialize the address structs */
	memset(&amp;proc_addr, 0, sizeof(struct sockaddr_nl));
	proc_addr.nl_family = AF_NETLINK;
	proc_addr.nl_pid = getpid(); // this process' PID
	proc_addr.nl_groups = CN_IDX_IWLAGN;
	memset(&amp;kern_addr, 0, sizeof(struct sockaddr_nl));
	kern_addr.nl_family = AF_NETLINK;
	kern_addr.nl_pid = 0; // kernel
	kern_addr.nl_groups = CN_IDX_IWLAGN;

	/* Now bind the socket */
	if (bind(sock_fd, (struct sockaddr *)&amp;proc_addr, sizeof(struct sockaddr_nl)) == -1)
		exit_program_err(-1, "bind");

	/* And subscribe to netlink group */
	{
		int on = proc_addr.nl_groups;
		ret = setsockopt(sock_fd, 270, NETLINK_ADD_MEMBERSHIP, &amp;on, sizeof(on));
		if (ret)
			exit_program_err(-1, "setsockopt");
	}

	/* Set up the "caught_signal" function as this program's sig handler */
	signal(SIGINT, caught_signal);

	/* Poll socket forever waiting for a message */
	while (1)
	{
		/* Receive from socket with infinite timeout */
		ret = recv(sock_fd, buf, sizeof(buf), 0);
		if (ret == -1)
			exit_program_err(-1, "recv");
		/* Pull out the message portion and print some stats */
		cmsg = NLMSG_DATA(buf);
		if (count % SLOW_MSG_CNT == 0)
			printf("received %d bytes: counts: %d id: %d val: %d seq: %d clen: %dn", cmsg->len, count, cmsg->id.idx, cmsg->id.val, cmsg->seq, cmsg->len);
		/* Log the data to file */
		l = (unsigned short)cmsg->len;
		l2 = htons(l);
		fwrite(&amp;l2, 1, sizeof(unsigned short), out);
		ret = fwrite(cmsg->data, 1, l, out);
		++count;
		if (count == 1)
		{
			/* Set alarm */
			/*alarm((*argv[2] - '0')); */
            alarm(atoi(argv[2]));
		}
		if (ret != l)
			exit_program_err(1, "fwrite");
	}

	exit_program(0);
	return 0;
}

void check_usage(int argc, char **argv)
{
	if (argc != 3)
	{
		fprintf(stderr, "Usage: %s <output_file> <time>n", argv[0]);
		exit_program(1);
	}
}

FILE *open_file(char *filename, char *spec)
{
	FILE *fp = fopen(filename, spec);
	if (!fp)
	{
		perror("fopen");
		exit_program(1);
	}
	return fp;
}

void caught_signal(int sig)
{
	fprintf(stderr, "Caught signal %dn", sig);
	exit_program(0);
}

void exit_program(int code)
{
	if (out)
	{
		fclose(out);
		out = NULL;
	}
	if (sock_fd != -1)
	{
		close(sock_fd);
		sock_fd = -1;
	}
	exit(code);
}

void exit_program_err(int code, char *func)
{
	perror(func);
	exit_program(code);
}

void exit_program_with_alarm(int sig)
{
	exit_program(0);
}

修改后的采数C代码可以实现自定义设定采数时长，时间参数单位为秒，可以设置10秒以上的数值。
该采数代码所在目录是:
~/linux-80211n-csitool-supplementary/netlink/
接着是编译该采数代码：

xxx:~$: cd ~
xxx:~$: cd linux-80211n-csitool-supplementary/netlink/
xxx:xxx$: make

编译后，当前目录会生成一个名为log_to_file的可执行文件，后面执行该文件（本文会用shell 脚本执行该文件）即可采数。

5、接着在采数设备上执行启动采数模式命令：

xxx:~$: cd ~
xxx:~$: sudo bash ./monitor.sh wlan0 64 HT20

执行上述命令后开始出现如下大片错误无需关注，最后会正常启动采数监听模式：

xxx@xxx:~$ sudo bash ./monitor.sh wlan0 64 HT20
[sudo] password for xxx: 
stop: Unknown instance: 
Bringing wlan0 down......
down: error fetching interface information: Device not found
wlan0: ERROR while getting interface flags: No such device
...
wlan0: ERROR while getting interface flags: No such device
Set wlan0 into monitor mode......
Bringing wlan0 up......
Set channel 64 HT20...
xxx@xxx:~$

6、在采数设备上执行循环采数的shell 脚本，shell 脚本make_data.sh 内容如下：

#!/bin/bash
# 存放数据的路径
org_p="/home/clife/csi_data/"
# 清空放数据的目录
dl=`rm -rf  ${org_p}*`

for i in {0..1000};
do
    echo "第${i}次采数"
    # 用整数命名数据文件
    fl_p="${org_p}${i}.dat"
    # 每采集30秒生成一个数据文件
    cais=`/home/clife/linux-80211n-csitool-supplementary/netlink/log_to_file $fl_p 30`
done

接着执行该shell 脚本启动采数：

xxx:xxx$: chmod +777 ./make_data.sh
xxx:xxx$: sudo ./make_data.sh

7、在采数设备上执行读取数据并分析的python 代码respiration_online.py，respiration_online.py 内容为：

# -*-coding:utf-8-*-
# -*-coding:utf-8-*-
import os
import time
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
# csi各种处理，参考宝藏工具：https://github.com/citysu/csiread
import csiread  # csiread/examples/csishow.py这里好多处理csi的基本操作，处理幅值和相位等等
import scipy.signal as signal
from PyEMD import EMD  #pip install EMD-signal
from scipy.fftpack import fft

# -----------------------------------------------求振幅和相位
# 参考：https://github.com/citysu/csiread 中utils.py和csishow.py
def scidx(bw, ng, standard='n'):
    """subcarriers index

    Args:
        bw: bandwitdh(20, 40, 80)
        ng: grouping(1, 2, 4)
        standard: 'n' - 802.11n， 'ac' - 802.11ac.
    Ref:
        1. 802.11n-2016: IEEE Standard for Information technology—Telecommunications
        and information exchange between systems Local and metropolitan area
        networks—Specific requirements - Part 11: Wireless LAN Medium Access
        Control (MAC) and Physical Layer (PHY) Specifications, in
        IEEE Std 802.11-2016 (Revision of IEEE Std 802.11-2012), vol., no.,
        pp.1-3534, 14 Dec. 2016, doi: 10.1109/IEEESTD.2016.7786995.
        2. 802.11ac-2013 Part 11: ["IEEE Standard for Information technology--
        Telecommunications and information exchange between systemsLocal and
        metropolitan area networks-- Specific requirements--Part 11: Wireless
        LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications
        --Amendment 4: Enhancements for Very High Throughput for Operation in
        Bands below 6 GHz.," in IEEE Std 802.11ac-2013 (Amendment to IEEE Std
        802.11-2012, as amended by IEEE Std 802.11ae-2012, IEEE Std 802.11aa-2012,
        and IEEE Std 802.11ad-2012) , vol., no., pp.1-425, 18 Dec. 2013,
        doi: 10.1109/IEEESTD.2013.6687187.](https://www.academia.edu/19690308/802_11ac_2013)
    """

    PILOT_AC = {
        20: [-21, -7, 7, 21],
        40: [-53, -25, -11, 11, 25, 53],
        80: [-103, -75, -39, -11, 11, 39, 75, 103],
        160: [-231, -203, -167, -139, -117, -89, -53, -25, 25, 53, 89, 117, 139, 167, 203, 231]
    }
    SKIP_AC_160 = {1: [-129, -128, -127, 127, 128, 129], 2: [-128, 128], 4: []}
    AB = {20: [28, 1], 40: [58, 2], 80: [122, 2], 160: [250, 6]}
    a, b = AB[bw]

    if standard == 'n':
        if bw not in [20, 40] or ng not in [1, 2, 4]:
            raise ValueError("bw should be [20, 40] and ng should be [1, 2, 4]")
        k = np.r_[-a:-b:ng, -b, b:a:ng, a]
    if standard == 'ac':
        if bw not in [20, 40, 80] or ng not in [1, 2, 4]:
            raise ValueError("bw should be [20, 40, 80] and ng should be [1, 2, 4]")

        g = np.r_[-a:-b:ng, -b]
        k = np.r_[g, -g[::-1]]

        if ng == 1:
            index = np.searchsorted(k, PILOT_AC[bw])
            k = np.delete(k, index)
        if bw == 160:
            index = np.searchsorted(k, SKIP_AC_160[ng])
            k = np.delete(k, index)
    return k

def calib(phase, k, axis=1):
    """Phase calibration

    Args:
        phase (ndarray): Unwrapped phase of CSI.
        k (ndarray): Subcarriers index
        axis (int): Axis along which is subcarrier. Default: 1

    Returns:
        ndarray: Phase calibrated

    ref:
        [Enabling Contactless Detection of Moving Humans with Dynamic Speeds Using CSI]
        (http://tns.thss.tsinghua.edu.cn/wifiradar/papers/QianKun-TECS2017.pdf)
    """
    p = np.asarray(phase)
    k = np.asarray(k)

    slice1 = [slice(None, None)] * p.ndim
    slice1[axis] = slice(-1, None)
    slice1 = tuple(slice1)
    slice2 = [slice(None, None)] * p.ndim
    slice2[axis] = slice(None, 1)
    slice2 = tuple(slice2)
    shape1 = [1] * p.ndim
    shape1[axis] = k.shape[0]
    shape1 = tuple(shape1)

    k_n, k_1 = k[-1], k[0]   # 这里本人做了修改，将k[1]改成k[0]了
    a = (p[slice1] - p[slice2]) / (k_n - k_1)
    b = p.mean(axis=axis, keepdims=True)
    k = k.reshape(shape1)

    phase_calib = p - a * k - b
    return phase_calib

# -----------------------------------------------EMD分解，去除高频噪声
# 参考：https://blog.csdn.net/fengzhuqiaoqiu/article/details/127779846
# 参考：基于WiFi的人体行为感知技术研究（南京邮电大学的一篇硕士论文）
def emd_and_rebuild(s):
    '''对信号s进行emd分解，去除前2个高频分量后，其余分量相加重建新的低频信号'''
    emd = EMD()
    imf_a = emd.emd(s)

    # 去掉前3个高频子信号，合成新低频信号
    new_s = np.zeros(s.shape[0])
    for n, imf in enumerate(imf_a):
        # 注意论文中是去除前2个，本人这里调整为去除前3个高频分量
        if n < 3:  
            continue
        new_s = new_s + imf
    return new_s

# -----------------------------------------------FFT变换筛选子载波
# 参考：https://blog.csdn.net/zhengyuyin/article/details/127499584
# 参考：基于WiFi的人体行为感知技术研究（南京邮电大学的一篇硕士论文）
def dft_amp(signal):
    '''求离散傅里叶变换的幅值'''
    # dft后，长度不变，是复数表示，想要频谱图需要取模
    dft = fft(signal)
    dft = np.abs(dft)
    return dft

def respiration_freq_amp_ratio(dft_s, st_ix, ed_ix):
    '''计算呼吸频率范围内的频率幅值之和,与全部频率幅值之和的比值
    dft_s: 快速傅里叶变换后的序列幅值
    st_ix: 呼吸频率下限的序号
    ed_ix: 呼吸频率上限的序号
    '''
    return np.sum(dft_s[st_ix:ed_ix])/np.sum(dft_s)

# ----------------------------------------------------------------------------- 均值恒虚警（CA-CFAR）
# 参考：https://github.com/msvorcan/FMCW-automotive-radar/blob/master/cfar.py
# 参考：基于WiFi的人体行为感知技术研究（南京邮电大学的一篇硕士论文）
def detect_peaks(x, num_train, num_guard, rate_fa):
    """
    Parameters
    ----------
    x : signal，numpy类型
    num_train : broj trening celija, 训练单元数
    num_guard : broj zastitnih celija，保护单元数
    rate_fa : ucestanost laznih detekcija，误报率

    Returns
    -------
    peak_idx : niz detektovanih meta
    """
    num_cells = len(x)
    num_train_half = round(num_train / 2)
    num_guard_half = round(num_guard / 2)
    num_side = num_train_half + num_guard_half

    alpha = 0.09 * num_train * (rate_fa ** (-1 / num_train) - 1)  # threshold factor

    peak_idx = []
    for i in range(num_side, num_cells - num_side):

        if i != i - num_side + np.argmax(x[i - num_side: i + num_side + 1]):
            continue

        sum1 = np.sum(x[i - num_side: i + num_side + 1])
        sum2 = np.sum(x[i - num_guard_half: i + num_guard_half + 1])
        p_noise = (sum1 - sum2) / num_train
        threshold = alpha * p_noise

        if x[i] > threshold and x[i] > -20:
            peak_idx.append(i)

    peak_idx = np.array(peak_idx, dtype=int)
    return peak_idx


if __name__ == '__main__':
    fs = 20  # 呼吸数据的采样率，设置为20Hz,数据包速率大于这个数的要进行下采样
    tx_num = 3
    rx_num = 3
    bpm_count_num = rx_num * tx_num * 2 * 10  # 理想情况下需要累加的呼吸速率个数

    is_sample = True   # 是否需要下采样
    sample_gap = 5     # 需要下采样则设置取数间隔

    # data_pt = 'E:/WiFi/test/data/csi_data/'
    data_pt = '/home/clife/csi_data/'

    while True:
        # 由于采数的shell脚本是不断产生30秒的数据文件的，为了不让数据文件撑爆硬盘，这里每次进入循环都要先删除多余的数据
        # 文件，留下最新的两个数据文件，因数据文件名是按整数来命名且依次递增的，文件名最大的两个文件是最新的文件。
        all_fl = sorted([int(item.split('.')[0]) for item in os.listdir(data_pt)])
        if len(all_fl)<2:
            time.sleep(2)
            continue
        for i in range(len(all_fl)-2):
            os.remove(data_pt+str(all_fl[i])+'.dat')
            
        # 取倒数第2个文件而不是最新的文件，可以确保拿到的文件已经采满30秒，而最新的数据文件可能正在写入数据。
        csifile = data_pt + str(all_fl[-2])+'.dat'  
        print('n', csifile)
        
        csidata = csiread.Intel(csifile, nrxnum=rx_num, ntxnum=tx_num, pl_size=10)
        csidata.read()
        csi = csidata.get_scaled_csi()
        print(csi.shape)

        # 等间隔抽样，为了将数据采样成20Hz，比如本人设置的发包率为100，那么sample_gap=5就可以降采样成20Hz
        if is_sample:
            csi = csi[0:-1:sample_gap,:,:,:]  
            print(csi.shape)


        # 振幅和相位计算
        csi_amplitude = np.abs(csi)                    # 求csi值的振幅
        csi_phase = np.unwrap(np.angle(csi), axis=1)   # 求csi值的相位
        csi_phase = calib(csi_phase, scidx(20, 2))     # 校准相位的值
        # print('csi_phase: ', csi_phase[:2, 1, 2, 1])

        # 中值滤波，去除异常点
        # 参考：https://blog.csdn.net/qq_38251616/article/details/115426742
        csi_amplitude_filter = np.apply_along_axis(signal.medfilt, 0, csi_amplitude.copy(), 3)   # 中值滤波,窗口必须为奇数，此处窗口为3
        csi_phase_filter = np.apply_along_axis(signal.medfilt, 0, csi_phase.copy(), 3)    # 中值滤波,窗口必须为奇数，此处窗口为3
        # print('csi_phase_filter: ', csi_phase_filter[:2, 1, 2, 1])

        # csi_amplitude_filter = csi_amplitude_filter[0:-1:5, :, :, :]
        # csi_phase_filter = csi_phase_filter[0:-1:5, :, :, :]
        # print(csi_phase_filter.shape)


        # emd分解信号-重建信号
        csi_amplitude_emd = np.apply_along_axis(emd_and_rebuild, 0, csi_amplitude_filter.copy())
        csi_phase_emd = np.apply_along_axis(emd_and_rebuild, 0, csi_phase_filter.copy())
        # print('csi_phase_emd: ', csi_phase_emd[:2, 1, 2, 1])


        # 基于振幅的fft变换筛选子载波，并针对挑选出的子载波进行寻峰和呼吸速率计算
        csi_dft_amp = np.apply_along_axis(dft_amp, 0, csi_amplitude_emd.copy())

        n = csi_dft_amp.shape[0]  # 采样点数
        # 0.15Hz对应dft中值的序号,呼吸频率下限
        l_ix = int(0.15*n/fs)
        # 0.5Hz对应dft中值的序号,呼吸频率上限
        u_ix = int(0.5*n/fs)+1
        # 计算呼吸频率值的占比
        csi_respiration_freq_ratio = np.apply_along_axis(respiration_freq_amp_ratio, 0, csi_dft_amp.copy(),l_ix, u_ix)
        # 针对1发1收对应的30个载波筛选出10个载波，进行呼吸频率计算
        sum_bpm = 0
        bpm_count = 0
        all_respiration_freq_ratio = 0
        for i in range(csi_respiration_freq_ratio.shape[1]):
            for j in range(csi_respiration_freq_ratio.shape[2]):
                temp = np.sort(csi_respiration_freq_ratio[:,i,j])
                for k in range(30):
                    if csi_respiration_freq_ratio[k,i,j] < temp[20]:  # 排名前10的才会进入下面的计算,如果temp[19]==temp[20]就会多出来一个
                        continue

                    amplitude_peak_idx = detect_peaks(csi_amplitude_emd[:, k, i, j].copy(), num_train=20, num_guard=8, rate_fa=1e-3)
                    phase_peak_idx = detect_peaks(csi_phase_emd[:, k, i, j].copy(), num_train=20, num_guard=8, rate_fa=1e-3)
                    amplitude_bpm = 0
                    phase_bpm = 0
                    try:
                        # 基于振幅计算的每秒呼吸次数
                        amplitude_bpm = (len(amplitude_peak_idx)-1)*fs/(amplitude_peak_idx[-1]-amplitude_peak_idx[0])
                        # 基于相位计算的每秒呼吸次数
                        phase_bpm = (len(phase_peak_idx)-1)*fs/(phase_peak_idx[-1]-phase_peak_idx[0])
                        # 呼吸心率必须大于1分钟9次
                        if ~pd.isna(amplitude_bpm) and ~pd.isna(phase_bpm) and amplitude_bpm > 0.15 and phase_bpm > 0.15:
                            sum_bpm = sum_bpm + amplitude_bpm + phase_bpm
                            bpm_count = bpm_count + 2
                            all_respiration_freq_ratio = all_respiration_freq_ratio + csi_respiration_freq_ratio[k,i,j]
                            # print(i, j, k, bpm_count)
                    except:
                        pass
                    # print(i, j, k, amplitude_bpm, phase_bpm)
        mean_respiration_freq_ratio = all_respiration_freq_ratio/bpm_count   # 呼吸频率范围的平均频率值
        print(bpm_count_num, bpm_count, round(mean_respiration_freq_ratio,4))
        # 下面的两个阈值需针对不同设备自行调整，本人自行采集了几个站位和静坐的数据以及几个无人情况下的数据，进行分析得出
        # 区分有人和无人的阈值
        if bpm_count/bpm_count_num > 0.7 and mean_respiration_freq_ratio > 0.03:  
            mean_bpm = sum_bpm / bpm_count
            print('rate ：', int(mean_bpm*60), '次/分钟')
        else:
            print('无人！')

以上代码各个功能模块都注释了参考出处，需要详细学习的可看参考链接或文献。

接着是执行该代码，进行持续的呼吸速率检测：

xxx:xxx$ /home/clife/anaconda3/bin/python37 respiration_online.py

以上命令中python37是本人设置的python.exe的软链接，知道是python 即可。

三、测试数据

链接：https://pan.baidu.com/s/1ZQIQT1bQot3-GOcnILS26g
提取码：1234

1、Demo计算本人的呼吸频率大致为21次/分钟，与标准成人的呼吸频率16~24次/分钟比较相符，如果你计算所得频率偏大，可以对数据进行进一步高频滤波（EMD分解后去掉更多高频分量）或者将FFT筛选子载波的频率范围缩小一些，使得最终用于CA-CFAR算法寻峰的载波曲线频率尽可能接近于呼吸信号的频率。
2、借助呼吸速率的计算，本Demo还可以在不同房间实现有人和无人的检测，有人时会给出呼吸速率值，无人则直接打印出无人结果，测试的case有：站立、坐椅子、坐地上、躺桌子上、躺地上，姿势方向有：平行2个设备的连接线、垂直两个设备的连接线。除了躺在地上无法检测出呼吸速率显示为无人的误报外，其他情形都可检测出呼吸速率，当人走出房间，显示为无人。多人场景也可以检测出呼吸速率。