代码收藏家 基于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 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 <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};
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个设备的连接线、垂直两个设备的连接线。除了躺在地上无法检测出呼吸速率显示为无人的误报外，其他情形都可检测出呼吸速率，当人走出房间，显示为无人。多人场景也可以检测出呼吸速率。
综上，本Demo检测出的呼吸速率可做参考，调整处理逻辑和参数可进一步改善结果，呼吸速率的误差不影响有人和无人的区分（除躺倒在地外），抗干扰能力较强，能适应不同环境。
