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

一、概述

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

采数用的是C代码

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

解析和分析数据用python代码

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

二、操作内容

1、配置好采数设备和代码运行环境，参考本人记录：
https://blog.csdn.net/Acecai01/article/details/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>*/
#include "iwl_connector.h"#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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 100int 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, bindstruct 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(&proc_addr, 0, sizeof(struct sockaddr_nl));proc_addr.nl_family = AF_NETLINK;proc_addr.nl_pid = getpid(); // this process' PIDproc_addr.nl_groups = CN_IDX_IWLAGN;memset(&kern_addr, 0, sizeof(struct sockaddr_nl));kern_addr.nl_family = AF_NETLINK;kern_addr.nl_pid = 0; // kernelkern_addr.nl_groups = CN_IDX_IWLAGN;/* Now bind the socket */if (bind(sock_fd, (struct sockaddr *)&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, &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: %d\n", 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(&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 %d\n", 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};
doecho "第${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 indexArgs: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—Telecommunicationsand information exchange between systems Local and metropolitan areanetworks—Specific requirements - Part 11: Wireless LAN Medium AccessControl (MAC) and Physical Layer (PHY) Specifications, inIEEE 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 andmetropolitan area networks-- Specific requirements--Part 11: WirelessLAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications--Amendment 4: Enhancements for Very High Throughput for Operation inBands below 6 GHz.," in IEEE Std 802.11ac-2013 (Amendment to IEEE Std802.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 kdef calib(phase, k, axis=1):"""Phase calibrationArgs:phase (ndarray): Unwrapped phase of CSI.k (ndarray): Subcarriers indexaxis (int): Axis along which is subcarrier. Default: 1Returns:ndarray: Phase calibratedref:[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.ndimslice1[axis] = slice(-1, None)slice1 = tuple(slice1)slice2 = [slice(None, None)] * p.ndimslice2[axis] = slice(None, 1)slice2 = tuple(slice2)shape1 = [1] * p.ndimshape1[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 - breturn 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:  continuenew_s = new_s + imfreturn 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 dftdef 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_halfalpha = 0.09 * num_train * (rate_fa ** (-1 / num_train) - 1)  # threshold factorpeak_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]):continuesum1 = 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_trainthreshold = alpha * p_noiseif x[i] > threshold and x[i] > -20:peak_idx.append(i)peak_idx = np.array(peak_idx, dtype=int)return peak_idxif __name__ == '__main__':fs = 20  # 呼吸数据的采样率，设置为20Hz,数据包速率大于这个数的要进行下采样tx_num = 3rx_num = 3bpm_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)continuefor 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就可以降采样成20Hzif 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/115426742csi_amplitude_filter = np.apply_along_axis(signal.medfilt, 0, csi_amplitude.copy(), 3)   # 中值滤波,窗口必须为奇数，此处窗口为3csi_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 = 0bpm_count = 0all_respiration_freq_ratio = 0for 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]就会多出来一个continueamplitude_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 = 0phase_bpm = 0try:# 基于振幅计算的每秒呼吸次数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_bpmbpm_count = bpm_count + 2all_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_countprint('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个设备的连接线、垂直两个设备的连接线。除了躺在地上无法检测出呼吸速率显示为无人的误报外，其他情形都可检测出呼吸速率，当人走出房间，显示为无人。多人场景也可以检测出呼吸速率。

综上，本Demo检测出的呼吸速率可做参考，调整处理逻辑和参数可进一步改善结果，呼吸速率的误差不影响有人和无人的区分（除躺倒在地外），抗干扰能力较强，能适应不同环境。