一种基于Hilbert-Huang变换时频分析的呼吸频率检测方法

本发明属于利用商用wifi设备实现感知的，具体为一种基于hilbert-huang变换时频分析的呼吸频率检测方法，用于通过处理信道状态信息(channel stateinformation,csi)信号实现非接触的呼吸频率检测。

背景技术：

1、随着生活水平的提高，日常的呼吸频率的监测引起了人们的关注。在医学领域，呼吸频率是衡量一个人健康状况的重要生理指标，呼吸相关疾病包括慢性阻塞性肺疾病、睡眠呼吸暂停综合征、夜间低通气综合征等。如果能较早地发现人体的呼吸频率的异常变化，就能为人们提供早期的医疗建议，从而提高健康水平。现有的研究方法分为可穿戴式的检测系统和非接触式的检测系统两类。然而，常用的可穿戴设备，如睡眠多导仪、呼吸带等，这类产品携带不便，不适用日常。因此，基于非接触式的方案是一个可行的选择。非接触性方案包括基于视频和基于无线技术。基于视频的方案由于可能存在潜在的隐私问题，易受光照影响，使用环境受到很大限制。而基于无线技术可以很好地保护隐私，对光线不敏感。因此，基于无线技术的呼吸频率检测得到了研究者的青睐。

2、在无线技术领域中，基于雷达的方案通常使用多普勒雷达测量由呼吸心跳等产生的周期性运动，但缺点是需要部署专用的硬件。基于(radio frequency identification,rfid)的方案通常是利用标签反射的电波信息对呼吸进行检测，但是rfid的标签信息容易在未经许可的情况下被非法读取，甚至被篡改。由于商用wifi设备造价低，并且目前wifi设备非常普及。因此，基于wifi信号进行呼吸频率检测已然成为研究的热门。并且基于wifi的方案是不需要额外的硬件设备部署的，因为它们是建立在现有的室内wifi设备的基础上。最初研究人员们使用wifi的接收信号强度指示(received signal strength indicator,rssi)信息对呼吸频率进行检测。受限于接收信号强度指数本身的粗粒度和信号波动的不确定性，早期基于wifi rssi的传感精度和应用范围都受到很大限制。除此之外，wifi设备上还有另一个可利用的信息是csi。与rssi不同，csi能够在物理层提供更细粒度的无线信道信息，因此被认为是精确传感的替代解决方案。csi包含不同子载波上的信道幅度和相位信息，能够区分多路径特性。这些信息揭示了载波随着传输距离的变化而发生的信号散射、反射和功率衰减现象，通过分析和研究csi的变化，推测导致通道状态变化的物理环境的变化，从而实现非接触式呼吸感知。

3、然而，现有的方法需要近距离检测。主要原因是csi感知依赖于微弱的反射信号，呼吸引起的细微信号变化容易被噪声所掩盖。另外，wifi收发器没有严格的时钟同步，会导致随机相位偏移。csi数据由于受到传输速率、传输功率以及wi-fi网卡内部csi参考电平波动的影响，会引入异常值和高频噪声。此外，现有的工作也存在盲点的问题，即当人体靠近感知装置时，在某些位置也无法有效地检测到呼吸。感知距离短和“盲点”约束极大地限制了现有方法的实际应用。

4、为解决上述问题，本发明提供一种基于hilbert-huang变换时频分析的呼吸频率检测方法。针对上述精度低、时变相位偏移和盲点的问题，选择对环境感知敏感度较强的wifi链路，利用csi比值消除时变相位偏移的同时对csi比值进行投影，充分利用幅度和相位信息，选择“好”的子载波。而对于由于设备和网卡引入的高频噪声的问题，使用特殊的信号分解以及hilbert-huang变换时频分析去除非人体呼吸频率分量后进行重构，并利用主成分分析融合，从而达到扩大感知范围和提高呼吸频率检测精度的目的，并且具有较强的鲁棒性。

技术实现思路

1、本发明的目的是在不需要穿戴任何设备或摄像头的情况以及人体目标存在相对干扰情况下，提供一种基于hilbert-huang变换时频分析的呼吸频率检测方法，提高呼吸频率检测精度。

2、本发明所述的一种基于hilbert-huang变换时频分析的呼吸频率检测方法，具体包括以下步骤：

3、步骤一：发送设备发送射频信号至接收设备，数据通过正交频分复用(orthogonalfrequency-division multiplexing,ofdm)传输使信号能够被调制以实现多个子载波的并行传输。利用intel5300网卡获取wifi信号csi信息，若wifi设备天线个数为i，子载波个数为k，数据包个数为n。则csi数据矩阵可以写为：

4、

5、其中，hi表示第i(1≤i≤i)条wifi链路的csi数据，表示第i(1≤i≤i)条wifi链路的第k(1≤k≤k)个子载波的csi数据，表示第n(1≤n≤n)个数据包中第i(1≤i≤i)条天线的第k(1≤k≤k)个子载波的csi数据。

6、

7、其中，和分别为第n个数据包中的第i个天线，第k子载波的实部和虚部，和表示第n个数据包中的第i个天线，第k子载波的幅值和相位。

8、步骤二：利用双方差，设置阈值挑选“敏感”的wifi链路。具体算法流程如下：

9、首先，计算不同链路上各子载波的csi幅值方差。i条天线，k个子载波的幅值方差可以表示为：

10、

11、其中，vi为第i(1≤i≤i)条天线的k个子载波的幅值方差的向量，为第i(1≤i≤i)条天线的第k(1≤k≤k)个子载波的幅值方差，表示为：

12、

13、其中，μi,k为第i(1≤i≤i)条天线的第k(1≤k≤k)个子载波的n个数据包幅值的均值，表示为：

14、

15、其次，设置阈值为εi＝0.7·max{vi}(1≤i≤i)，其中，max{·}表示对{·}取最大值，选择每条wifi链路幅值方差大于阈值的子载波，子载波个数为ki，ki＝card{k|{vi}＞εi}，其中，card{·}表示计算{·}元素的个数。并对集合中子载波的幅值方差取平均值，则i条天线选择子载波的平均幅值方差记为：v＝[v1,…,vi,…,vi]。其中，vi表示第i(1≤i≤i)条wifi链路选择子载波的平均幅值方差，是一个数值。表示为：

16、

17、最后，挑选出平均幅值方差最大和次大的wifi链路，重新命名为h′1和h′2。

18、步骤二：对于步骤一得到的两条“敏感”wifi链路，首先，采用csi比值法消除时变相位偏移，将两个接收天线子载波的csi读数相除，csi比值计算公式如下：

19、

20、其中，xk表示第k个子载波的csi比值数据，和分别为选出的两个wifi链路h′1和h′2的第k个子载波的csi数据。

21、其次，csi比值仍然是复数，通过monitor模式采集的csi数据可能会丢包，丢包率约为0.1％-0.05％，会影响原始数据的准确性。进一步地，为了使采集的数据更加准确。首先，对原始数据进行线性插值。其次，将每个子载波经过hampel去除异常值，并将去除异常值后的数据再次经过savitzky-golay滤波器平滑滤波。本发明将处理过的csi比值数据重构为：

22、x′＝[x′1,…,x′k,…,x′k]

23、其中，xk′表示为xk滤波后的第k个子载波的csi比值数据。

24、

25、其中，和分别表示csi比值数据的第k个子载波x′k的实部和虚部，和分别代表csi比值数据的第k个子载波的第n个数据包的实部和虚部。

26、步骤三：投影csi比值数据，提取呼吸特征短期呼吸噪声比(breathing-to-noiseratio,bnr)，挑选满足bnr要求的子载波。具体算法流程如下：

27、对步骤二得到的滤波后的csi比值信号进行投影，结合幅度和相位信息生成代表不同呼吸模式信号的候选集。

28、首先，将滤波后的csi比值的时间序列x′k(1≤k≤k)投影到旋转轴[cosθsinθ]，从0逐渐增大参数θ，生成不同的候选序列：

29、

30、其中，yk表示子载波x′k在投影坐标轴的旋转角度为θ时生成的候选序列。设置投影角度θ范围为[0,π]，步长设为π/p，每一个子载波共生成p个候选时间序列。所有子载波投影后生成的不同候选序列表示如下：

31、

32、其中，表示第k个子载波的第p(0≤p≤p-1)个候选序列，即：第k个子载波在投影角度为θ＝(p-1)·π/p时，得到的候选时间序列。

33、其次，计算信号短期呼吸噪声比来提取呼吸特征用来衡量组合候选的性能。信号模式在一段时间内的周期性代表了其感知呼吸的能力，具体计算步骤如下：

34、首先，若采样率fs＝100hz，采用12的窗长，对应1200个样本，用6992个零值样本加零填充，得到8192个样本，并对所有样本进行傅里叶变换。其次，在人体呼吸(10-37bpm)范围内筛选出最大能量，即为呼吸信号能量，将该能量与频域的总能量进行比值，得到短期呼吸噪声比。

35、求k个子载波的所有候选时间序列的呼吸特征bnr，表示如下：

36、

37、其中，表示第k个子载波在投影角度为θ＝(p-1)·π/p时，得到的候选时间序列的bnr值，是一个数值。

38、然后，寻找第k个子载波在不同投影轴下最大的bnr值，即组成集合{bk}。取集合{bk}中最大值为b，并设置阈值为δ＝0.7b，选择集合{bk}中大于该阈值的bk所对应的子载波候选序列。选择子载波的个数为：s＝card{k|{bk}＞δ}，并寻找对应集合{k|{bk}＞δ}中第k个子载波的候选时间序列。选择结果如下：

39、y′＝[y′1,…,y′s,…,y′s]

40、其中，y′s表示选择的第s(1≤s≤s)个子载波对应的候选时间序列。

41、步骤四：信号变分模态分解(variational mode decomposition,vmd)和hilbert-huang变换时频分析，去除非人体呼吸频率分量进行重构。具体算法流程如下：

42、首先，将步骤三中挑选得到的s个子载波分别进行变分模态分解：

43、

44、其中，y′s(t)为y′内的第s个子载波，m为模态分量的数量，um(t)为第m(1≤m≤m)个模态分量信号。

45、其次，对每个模态分量构造其解析信号，解析信号是基于原信号和其希尔伯特变换的结合，根据解析信号得到瞬时频率和瞬时能量：

46、

47、其中zm表示解析信号，表示希尔伯特变换，am和θm分别表示第m个模态分量的瞬时振幅和瞬时相位，|am(t)|2和wm分别表示第m个模态分量的瞬时能量和瞬时频率。

48、然后，通过上式求得的瞬时频率和瞬时能量对每个模态分量进行hilbert-huang变换时频分析，去掉前两个与呼吸频率无关的高频分量，对剩余呼吸范围内的模态分量进行重构，重构后的第s(1≤s≤s)个子载波时间序列表示为：

49、

50、重构后的所有子载波组合记为y″：

51、y″＝[y″1,…,y″s,…,y″s]

52、步骤五：信号融合。具体算法如下：

53、利用主成分分析(principal component analysis,pca)融合步骤四中重构之后的信号，pca将原始特征向量变换成一组线性无关的主成分。该算法首先集中计算所有样本的协方差矩阵，然后通过奇异值分解求解特征值和特征向量。最后利用最大特征值和特征向量得到主呼吸信号，作为主成分分析信号：

54、

55、其中，an表示融合后的主成分分析信号的第n(1≤n≤n)个数据包。

56、步骤六：对步骤五得到的主成分呼吸信号进行峰值检测，提取呼吸速率。具体算法流程如下：

57、首先，对主成分呼吸信号获取局部峰值集：maxset＝{τj,1≤j≤j}。验证窗口设为：m；

58、其次，对于每一个峰值点τj找到位置和峰值locs:＝locatio(τj)；amp:＝amplitude(τj)；

59、再次，对该范围的每个点首先判断n是否在1＜n＜n范围内。如果在，通过对比该范围的每一个点来判断是否有比当前峰值点更大的点，若：amp＜x(n)，从marset中删除τk；并更新marset；

60、最后，由更新后的maxset＝{τl,1≤l≤l}得到真峰相邻真峰之间的总时间间隔，记为sum；由公式：得到呼吸速率，fs为采样率。

61、有益效果

62、本方法首先将得到的三条wifi链路，通过双方差挑选波动程度较大的两条wifi链路，将其进行比值消除商用wifi随时间变化的随机相位偏移问题。其次，投影csi比值，提取呼吸特征bnr，并挑选满足阈值的子载波进行vmd信号分解和hilbert-huang变换进行时频分析，进而去除非人体呼吸频率分量。在此基础上，进行重构，并利用主成分分析融合重构的信号。最后，用峰值检测提取呼吸速率。

63、本发明提供一种基于hilbert-huang变换时频分析的呼吸频率检测方法。在去除时变相位偏移及大量高频噪声、选择最优子载波的同时，不仅扩大了感知范围、还消除了“盲点”，并提高了呼吸频率检测精度，为非接触呼吸频率检测提供了更为精确的检测方法。