时频分析-呼吸心跳信号检测方法(六)

发布于:2022-12-02 ⋅ 阅读:(621) ⋅ 点赞:(0)

系列文章目录

《呼吸心跳信号检测方法(一)》

《数据采集-呼吸心跳信号检测方法(二)》

《信号建模-呼吸心跳信号检测方法(三)》

《滤波算法-呼吸心跳信号检测方法(四)》

《MUSIC算法-呼吸心跳信号检测方法(五)》


文章目录

 前言

一、时频分析工具

1.1 短时傅里叶变换

1.2 小波变换

二、实测数据仿真

2.1 短时傅里叶分析

2.2 连续小波变换

总结


 前言

       由于呼吸与心跳在较长时间内可能并不是平稳的,通过傅里叶分析得到的频谱图失去了时间维度上的信息,并不能完全反映呼吸与心跳的变化特征,为了检测呼吸与心跳频率随时间的变化规律,需要考虑利用时频分析工具对呼吸心跳信号的时频特性进行分析。


一、时频分析工具

1.1 短时傅里叶变换

       短时傅里叶变换的基本原理是通过窗函数的滑动截取待处理信号不同时间段的信号成分,并对截取的信号进行傅里叶分析得到不同时间段的频率信息,以此来观测呼吸与心跳频率变化的规律。短时傅里叶变换公式如下:

W\left ( f,\tau \right )=\int w\left ( t-\tau \right )f\left ( t \right )e^{-j\omega t}dt

其中w\left ( t \right )是窗函数。可以看出,短时傅里叶变换需要确定的参数主要包括:窗函数类型,窗函数长度以及窗函数滑动步长。

1.2 小波变换

       虽然短时傅里叶变化可以得到微动信号的时频特征信息,但由于时频信息的不确定性原理,短时傅里叶变换得到的结果并不能同时满足较高的时间分辨率和频率分辨率的要求。

       小波变换在一定程度可以解决这个问题。虽然小波变换仍然摆脱不了不确定性原理的桎梏,但其可以提供多个分辨率尺度下的尺度-平移图像,并且由小波变换得到的尺度-平移图像,在低尺度(高频率)处有较高的时间分辨率,在高尺度(低频)处有较高的细节(频率)分辨率。而一般的信号(包括本实验得到的微动信号)处理都需要在低尺度处有较高的时间分辨率,在高尺度处有较高的细节分辨率。下面简单介绍小波变换的基本概念,并根据本论文待处理的信号特点,合理利用小波变换分析实验数据。 

  •  连续小波变换  

        与傅里叶变换类似,小波变换也涉及到基函数的概念,但与傅里叶变换的基函数不同的是:傅里叶变换的基函数覆盖范围为整个时间轴,而小波变换的基函数的持续时间有限(这也是小波名称的由来),由此,小波的变换的基函数需要两个自变量来描述:一个是尺度因子(用a​表示),用来描述基函数的震动快慢的量(对应于短时傅里叶变换的频率量);另一个是平移因子(用b​表示),用来描述基函数在时间轴的位置。由此,连续小波变换的定义如下:

 W_{f}\left ( a,b \right )=\left \langle f\left ( t \right ),\Psi _{a,b}\left ( t \right ) \right \rangle=\int_{-\infty }^{\infty}f\left ( t \right )\Psi _{a,b}^{*}\left ( t \right ) dt​ 

其中a>0,f\in L^{2}\left ( R \right ),式中系列函数\Psi _{a,b}\left ( t \right )=\frac{1}{\sqrt{a}}\Psi\left ( \frac{t-b}{a} \right ) 称为小波函数,它是由母函数\Psi \left ( t \right )经过不同的时间尺度伸缩和不同的时间平移得到的。对应的小波逆变换为:

f\left ( t \right )=C_{\Psi }^{-1}\int_{0}^{\infty }\frac{1}{a^{2}}\int_{-\infty}^{\infty}W_{f}\left ( a,b \right )\Psi\left ( \frac{t-b}{a} \right )db\, da

式中 C_{\Psi }=\int_{-\infty }^{\infty}\frac{\left | \Psi \left ( \omega \right ) \right |^{2}}{\left | \omega \right |}d \omega为小波容许性常数,满足0<C_{\Psi }<\infty,\Psi \left ( \omega \right )为\Psi \left ( t \right )的傅里叶变换。

       从定义上可以看出,连续小波变换的基函数是由母小波函数的平移和收缩得到的,因此母小波的性能将直接影响到小波变换的质量,为此母小波函数必须满足一定的条件。下面对一些经典小波进行简单介绍。

    1)Haar小波

       Haar函数是小波分析中最早用到的小波函数,也是最简单的一个小波函数,表达式如下:

\Psi _{H}\left ( t \right )=\left\{\begin{matrix} 1 & 0\leq t\leq \frac{1}{2}\\ -1 & \frac{1}{2}\leq t\leq 1\\ 0 & else \end{matrix}\right.

        对应的波形如下图所示,从图中可以看出,该小波在时域上虽然是紧支撑的,但是不连续的,在本项目中作为基本小波性能不好。

 haar小波

2)Mexico草帽小波

          Mexican Hat函数为Gauss函数的二阶导数。表达式如下:

\Psi \left ( t \right )=\frac{2}{\sqrt{3}}\pi^{-1/4}\left ( 1-t^{2} \right )e^{-t^{2}/2}

 Mexico草帽小波

         上图为对应的波形,可以看出,它在时域有较好的局部化。

3)Morlet小波

          Morlet小波是高斯包络下的单频率正弦函数,表达式如下:

\Psi \left ( t \right )=\pi^{-1/4}e^{-t^{2}/2}\cos \left ( 5t \right )

 Morlet小波

4)Gaussian 小波

        高斯小波是由一基本高斯函数分别求导而得到的,表达式为:

\Psi \left ( t \right )=c\frac{d^{k}}{dt^{k}}e^{-t^{2}/2}

        当求导阶次k取偶数时,\Psi \left ( t \right )正对称,当求导阶次k取奇数时,\Psi \left ( t \right )反对称。下图为k=4时的的时域波形。

 Gaussian小波

        除了上述经典类小波,还有通过尺度函数的加权组合产生的正交小波,以及放宽正交性条件下构造的双正交小波。下面是各小波在正交性、紧支撑性以及消失矩阶数方面对比的结果。

常见小波函数性能对比表

小波函数

Haar

Daubechies

Biorthogoal

Coiflets

Symlets

Morlet

Mexican Hat

Meyer

正交性

有

有

无

有

有

无

无

有

紧支撑性

有

有

有

有

有

无

无

无

连续小波变换

可以

可以

可以

可以

可以

可以

可以

可以

离散小波变换

可以

可以

可以

可以

可以

不可以

不可以

可以但无FWT

支撑长度

1

2N-1

重构:2Nr+1分解:2Nd+1

6N-1

2N-1

有限长

有限长

有限长

小波函数消失矩阶数

1

N

Nr-1

2N

N

-

-

-

二、实测数据仿真

2.1 短时傅里叶分析

        由于最终处理的信号是采样后的离散信号,所以需要确定窗函数的窗长以及窗函数类型。对于滤波前的微动信号,由于呼吸信号成分为主要信号成分,考虑到呼吸频率在0.2-0.5 Hz,采样频率为17 Hz。所以处理滤波前的微动信号时,窗函数的长度定为128。窗函数类型选择为矩形窗。对于滤波后的微动信号,由于心跳信号成分为主要信号成分,考虑心跳频率为0.9-1.6 Hz,所以处理滤波后的微动信号时,窗函数的长度定为32。窗函数类型选择为汉宁窗。

        这里只针对单人单次数据进行时频分析,采用的是短时傅里叶变换算法,为了使得到的时频信息更加直观,本文通过对时频图像进行处理,得到最大频率变化过程以及频谱质心图。频谱质心图可以通过下式计算得到。

centroid\left ( t \right )=\frac{\int_{-\infty }^{\infty}x\left ( \omega \right )\left | W\left ( f,t \right ) \right |^{2}d\omega }{\int_{-\infty }^{\infty}\left | W\left ( f,t \right ) \right |^{2}d\omega}

其中x\left ( \omega \right )为加权函数,本文取x\left ( \omega \right )=\omega,W\left ( f,t \right )为时频数据。

        对未经滤波的微动信号进行短时傅里叶变换,得到下图所示的时频图像。

滤波前由短时傅里叶变换得到的时频图

       从上述时频图可以看出,滤波前信号的频谱成分随时间是变化的,一方面频谱的能量随时间是变化的,另一面频谱的最大频率值也是变化的。下图是对应的最大频率变化曲线以及频谱质心图。

最大频率变化曲线                                                 频谱质心曲线

       从最大频率变化曲线看,呼吸频率在0.16-0.21 Hz频率范围内无规律变化。从频谱质心图看,呼吸频率在0.16-0.32 Hz范围内无规律变化。

       对滤波后的微动信号进行短时傅里叶变换,得到下图所示的时频图像。

滤波后由短时傅里叶变换得到的时频图

        从时频图可以看出,滤波前信号的频谱成分随时间是变化的,一方面频谱的能量随时间是变化的,另一面频谱的最大频率值也是变化的。下图是对应的最大频率变化曲线以及频谱质心图。

最大频率变化曲线                                                         频谱质心曲线

        从最大频率变化曲线看,呼吸频率在1-2 Hz频率范围内波动。从频谱质心图看,呼吸频率在1.1-1.8 Hz范围内波动。

2.2 连续小波变换

       考虑到待处理的信号是连续变化的信号,所以选择的小波函数应该是连续的,另外考虑到基函数的性能标准,本文采用Symlets系列函数(这里考虑sym5)作为母小波函数。由于小波变换得到的尺度-平移三维图像,根据信号频率与尺度成反比,再加上一般待处理的信号都是离散的,为了从连续小波变换得到的尺度-平移图像中得到信号的时频图像,需要对尺度s做如下变换得到频率信息。

f=\frac{f_{s}}{s}

式中s为尺度因子,f_{s}为采样频率,f为s对应的频率点。

         对未经滤波的微动信号进行连续小波变换,再经过上式映射处理,可以得到下图所示的时频图。

                                         时频图                                                                          最大频率变化曲线

        从上图可以看出呼吸频率大概在0.2 Hz附近波动。

        对滤波后的微动信号进行连续小波变换,再经过上式映射处理,可以得到下图所示的时频图。

                                               时频图                                                                   最大频率变化曲线

 从上图可以看出,心跳频率大概在1.5 Hz左右波动。

代码见:《智慧医疗+生物雷达+呼吸心跳检测(包含采集的数据)》


总结

       本文主要借用短时傅里叶变换以及小波变换对呼吸心跳信号的时频特性进行了分析。至此,呼吸心跳信号检测方法系列文章已从数据的采集、信号的建模、滤波算法、MUSIC算法、时频分析等角度进行了阐述,此内容为本人毕设研究内容,如有争议,请在评论区留言。转载请附上链接:【杨(_> <_)】的博客_CSDN博客-信号处理,SAR,代码实现领域博主


微信公众号

今日签到

点亮在社区的每一天
去签到