基于WiFi的CSI数据做呼吸频率检测 |
您所在的位置：网站首页 › psi数据和pos数据的区别 › 基于WiFi的CSI数据做呼吸频率检测 |
基于WiFi的CSI数据做呼吸频率检测
|
一、概述
本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、布设试验场景： 以上命令的含义，参考本大节第1步骤的配置记录博客。 此时设备会按每秒100个数据帧的速率持续发送数据，以上命令设置的发送数据量够发115天，要中断发送直接按ctrl+c即可。 4、选择另一台接收数据的设备，将本人修改的采数C代码log_to_file.c替换掉原先的log_to_file.c，先看修改后的log_to_file.c： /* * (c) 2008-2011 Daniel Halperin */ #include "iwl_connector.h" #include #include #include #include #include #include #include #include #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(&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(&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 *)&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 \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}; 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.sh7、在采数设备上执行读取数据并分析的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 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) 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个设备的连接线、垂直两个设备的连接线。除了躺在地上无法检测出呼吸速率显示为无人的误报外，其他情形都可检测出呼吸速率，当人走出房间，显示为无人。多人场景也可以检测出呼吸速率。 综上，本Demo检测出的呼吸速率可做参考，调整处理逻辑和参数可进一步改善结果，呼吸速率的误差不影响有人和无人的区分（除躺倒在地外），抗干扰能力较强，能适应不同环境。 |
3、选择一台发射数据的设备，输入如下发数据的命令：【本文地址】