滤波算法-呼吸心跳信号检测方法（四）

最新推荐文章于 2024-03-21 09:52:56 发布

【杨（_> <_)】

最新推荐文章于 2024-03-21 09:52:56 发布

阅读量3.7k

点赞数 2

分类专栏：信号处理 # 呼吸心跳代码实现

本文链接：https://blog.csdn.net/yjh_2019/article/details/127026508

版权

信号处理同时被 3 个专栏收录

41 篇文章 161 订阅

订阅专栏

代码实现

33 篇文章 27 订阅

订阅专栏

呼吸心跳

6 篇文章 11 订阅

订阅专栏

系列文章目录

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

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

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

文章目录

前言

上文主要介绍呼吸心跳的信号模型，从人体微动的角度对心跳信号进行分离相比于直接从雷达回波信号中分离心跳信号更加简单，论文关注的是由复基带信号相位信息得到的人体体表微动信号。本文主要介绍通过滤波的方式分离出呼吸与心跳信号，并得到呼吸心跳频率的估值。且滤波器的设计基于窗函数进行设计。

一、窗函数

窗函数设计方法的基本思想为，首先根据实际需要确定滤波器的传输函数 $H_{d}\left ( e^{j\omega } \right )$ ，然后用窗函数截取该滤波器的单位脉冲响应，通过截取得到部分单位脉冲响应 $h\left ( n \right )$ 构造所需的FIR滤波器。该设计方法由于截断误差的存在，得到的FIR滤波器与待设计的滤波器是有偏差，为使设计的FIR滤波器的单位脉冲响尽量接近待设计的滤波器的单位脉冲响应，需要选择合适的窗函数。

为了确定本设计选择的窗函数的类型以及传输函数的长度，下面对一些常见的窗函数进行简要分析。由于MATLAB有自带的窗函数程序：

Boxcar 矩形窗函数

Hanning 汉宁窗函数

Hannming 海明窗函数

Blackman 布莱克曼窗函数

Kaiser 凯瑟窗函数

有关的窗函数性能分析通过MATLAB仿真实现。

1.1 矩形窗

矩形窗的时域形状是矩形的窗，其时域表达式为：

$w\left ( n \right )=R_{N}\left ( n \right )=\left\{\begin{matrix} 1 & 0\leq n\leq N-1\\ 0 & else \end{matrix}\right.$

则对应的频率响应为：

$W_{R}\left ( e^{j\omega } \right )=\frac{\sin \left ( \omega N/2 \right )}{\omega /2}e^{-j\frac{N-1}{2}\omega }$

N=101的幅频响应曲线如下图所示，从图中可以看出矩形窗的频谱副瓣约为-13.4 dB。

1.2 汉宁窗

汉宁窗的时域图为余弦函数形，所以又称升余弦窗，时域表达式为：

$w\left ( n \right )=\frac{1}{2}\left [ 1-\cos \left ( \frac{2\pi n}{N-1} \right ) \right ]R_{N}\left ( n \right )$

记 $W_{R}\left ( \omega \right )=\frac{\sin \left ( \omega N/2 \right )}{\omega /2}$ ，则对应的频率响应为：

$W_{hanning}\left ( e^{j\omega } \right )=\left [ 0.5W_{R}\left ( \omega \right ) +0.25W_{R}\left ( \omega-\frac{2\pi}{N-1} \right )+0.25W_{R}\left ( \omega +\frac{2\pi}{N-1} \right )\right ]e^{-j\frac{N-1}{2}\omega }$

N=101的幅频响应曲线如右下图所示，从图中可以看出对应的副瓣约为-31.5 dB。

左图为汉宁窗时域图，右图为汉宁窗频谱

1.3 海明窗

海明窗的时域形状也是余弦形，但与汉宁窗又有区别，时域表达式为：

$w\left ( n \right )=\left [ 0.54-0.46\cos \left ( \frac{2\pi n}{N-1} \right ) \right ]R_{N}\left ( n \right )$

对应的频率响应为：

$W_{hamming}\left ( e^{j\omega } \right )=\left [ 0.54W_{R}\left ( \omega \right ) +0.23W_{R}\left ( \omega-\frac{2\pi}{N-1} \right )+0.23W_{R}\left ( \omega +\frac{2\pi}{N-1} \right )\right ]e^{-j\frac{N-1}{2}\omega }$

N=101的幅频响应曲线如右下图所示，从图可以看出海明窗的副瓣约为-42.6 dB。

左图海明窗时域图，右图为海明窗频谱

1.4 布莱克曼窗

布莱克曼窗是利用了二阶成分的余弦函数构成的，其时域表达式为：

$w\left ( n \right )=\left [ 0.42-0.5\cos \left ( \frac{2\pi n}{N-1} \right )+ 0.08\cos \left ( \frac{4\pi n}{N-1} \right )\right ]R_{N}\left ( n \right )$

对应的频率响应为：

$W_{blackman}\left ( e^{j\omega } \right )=\left \{ 0.42W_{R}\left ( \omega \right ) +0.25\left [ W_{R}\left ( \omega-\frac{2\pi}{N-1} \right )+W_{R}\left ( \omega +\frac{2\pi}{N-1} \right ) \right ]+0.04\left [ W_{R}\left ( \omega-\frac{4\pi}{N-1} \right )+W_{R}\left ( \omega +\frac{4\pi}{N-1} \right ) \right ]\right \}e^{-j\frac{N-1}{2}\omega }$

N=101的幅频响应曲线如右下图所示，从图中可以看出布莱克曼窗的副瓣约为-58.2 dB。

左图为时域图形，右图为布莱克曼窗频谱

1.5 凯瑟窗

凯瑟窗可以通过调节参数 $\beta$ 的值改变窗的形状，时域表达式为：

$w\left ( n \right )=\frac{I_{0}\left ( \beta \sqrt{1-\left ( 1-\frac{2n}{N-1} \right )^{2}} \right )}{I_{0}\left ( \beta \right )}R_{N}\left ( n \right )$

式中 $I_{0}\left ( \cdot \right )$ 是第一类零阶变型贝赛尔函数， $\beta$ 是可以自由选择的参数。当参数 $\beta$ 为不同值时凯瑟窗的时域图和对应的频谱图如下图所示，从图中不难发现：

$\beta =0$ ，相当于矩形窗； $\beta =5.44$ ，相当于海明窗，但凯瑟窗能量更多集中在主瓣中； $\beta =7.8$ ，相当于布莱克曼窗。

此外，从图中可以看出，该参数可以同时调整主瓣宽度与旁瓣衰减。 $\beta$ 越大，则 $w\left ( n \right )$ 窗越窄，窗谱的旁瓣衰减越大，主瓣宽度也相应的增加。

不同值下的凯瑟窗时域图

对应的凯瑟窗频谱

事实上，除了上述几种窗函数，还有以下几种窗函数，其中下表是他们的性能对比表。

常见窗函数性能表

名称	滤波器过渡带宽	最小阻带衰减
矩形	1.8π/M	21dB
巴特利特	6.1π/M	25dB
汉宁	6.2π/M	44dB
海明	6.6π/M	51dB
布莱克曼	11π/M	74dB
BOHMANWIN	5.8π/M	51.5dB
NUTTALLWIN	15.4π/M	108dB
PARZENWIN	6.6π/M	56dB
FLATTOPWIN	19.6π/M	108dB
GAUSSWIN	5.8π/M	60dB
BARTHANNWIN	3.6π/M	40dB
BLACKMANHARRIS	16.1π/M	109dB
CHEBWIN	15.2π/M	113dB
TUKEYWIN	2.4π/M	22dB

从上表可以看出，高斯窗的综合性能较好。

二、滤波器类型

2.1 高通滤波器

考虑到心跳的频率大于呼吸的频率，呼吸信号强度远大于心跳信号的强度，由此可以得到以下结论，由雷达回波的相位信息得到的体表微动信号，呼吸信号分量主要集中在低频段，心跳信号分量主要在较高频段，根据呼吸与心跳频率的先验信息，可以设计高通滤波器，将呼吸信号分量滤除，剩余的信号分量则主要由心跳信号构成。这里采用窗函数方法设计高通滤波器。

设计步骤：

1）根据实际信号特点，给定待设计的滤波器的频率响应 $H_{d}\left ( e^{j\omega } \right )$ ，使其满足通带上具有单位增益和线性相位，在阻带上具有零响应。对于高通滤波器，需要确定的参数有截止频率，并且截止频率 $\omega _{c}$ 为的高通滤波器的频率响应可以由下式给定。

$H_{d}\left ( e^{j\omega } \right )=\left\{\begin{matrix} e^{-j\alpha \omega } & \omega _{c}\leq \left | \omega \right |\leq \pi\\ 0 & \left | \omega \right |< \omega _{c}\end{matrix}\right.$

2）通过对滤波器频域上的传输函数傅里叶逆变换，可以得到这个滤波器时域上的单位脉冲响应

$h_{d}\left ( n \right )=\frac{\sin \left [ \pi\left ( n-\alpha \right ) \right ]}{\pi\left ( n-\alpha \right ) }-\frac{\sin \left [ \omega _{c}\left ( n-\alpha \right ) \right ]}{\pi\left ( n-\alpha \right ) }$

为了使得到的滤波器单位脉冲响应满足因果条件和相位线性条件，可令

$\alpha =\frac{N-1}{2}$

3）由于待设计滤波器的单位脉冲响应持续时间长，但能量较为集中，因此可用窗函数截取 $h_{d}\left ( n \right )$ ，得到所设计的近似FIR滤波器h(n)

$h\left ( n \right )=h_{d}\left ( n \right )w\left ( n \right )$

根据上述分析，高通滤波器设计需要的参数包括窗函数类型，窗函数的长度以及截止频率。

2.2 带通滤波器

由于呼吸信号强度较大，检测较为容易，而心跳信号微弱，一般难以检测，因此可以考虑根据先验信息，只保留某一频段的信号成分，通过滤出后的信息来检测心跳信号的频率参数。这里同样采用窗函数方法设计带通滤波器。

设计步骤：

1）与高通滤波器设计的第一步基本一致，区别在于，带通滤波器的设计需要确定的参数主要有带通范围，带通范围为 $\left ( \omega _{p}, \omega _{u}\right )$ 的带通滤波器由下式给定。

$H_{d}\left ( e^{j\omega } \right )=\left\{\begin{matrix} e^{-j\alpha \omega } & \omega _{p}\leq \left | \omega \right |\leq \omega _{u}\\ 0 &else\end{matrix}\right.$

2）通过对滤波器频域上的传输函数傅里叶逆变换，可以得到这个滤波器时域上的单位脉冲响应

$h_{d}\left ( n \right )=\frac{\sin \left [ \omega _{u}\left ( n-\alpha \right ) \right ]}{\pi\left ( n-\alpha \right ) }-\frac{\sin \left [ \omega _{p}\left ( n-\alpha \right ) \right ]}{\pi\left ( n-\alpha \right ) }$

与高通滤波器设计的目的一样，需要令

$\alpha =\frac{N-1}{2}$

3）与高通滤波器设计时一样，用窗函数截取 $h_{d}\left ( n \right )$ ，得到所设计的近似FIR滤波器h(n)

$h\left ( n \right )=h_{d}\left ( n \right )w\left ( n \right )$

根据上述分析，带通滤波器设计需要的参数包括窗函数类型，窗函数的长度以及带通范围。

三、实测数据处理

3.1 高通滤波器

为设计高通滤波器，需要设置以下参数：截止频率，滤波器的阶数，加窗类型。由于心跳频率为0.9-1.6 Hz，所以本论文设计的高通滤波器的的截止频率为1 Hz；从下图可以看出呼吸与心跳信号强度差32dB左右，所以为了抑制呼吸信号对心跳信号的干扰，最小阻带衰减至少为32dB,理论上最小阻带衰减越大越好，本设计采用高斯窗作为滤波器设计的窗函数；考虑到呼吸与心跳频率差值大约为1Hz，根据慢时间维的采样速率为17 FPS，则滤波器的阶数M满足 $1/17\times 2\pi>5.8\pi/M$ ，由此可得滤波器阶数 $M>49.3$ ，这里选择 $M=50$ 。根据上面的参数设置可以得到下图所示的高通滤波器。

将该滤波器运用到实验得到的人体微动信息，结果如下图所示，可以看出，在滤波之前，呼吸信号很大，易检测到呼吸频率为0.1842 Hz,滤波后可以看出，呼吸信号成分得到削弱，可以很明显的检测出心跳信号的频率为 1.309 Hz。

从过滤后的频谱图可以看出，心跳信号明显强于其他干扰信号，但过滤后的信号成分中还包含较强成分的心跳的二次谐波成分，为了提高心跳频率的检测正确率，可以采用带通滤波的处理方式去除呼吸信号以及较强的心跳谐波分量。

3.2 带通滤波器

为设计带通滤波器，同样需要设置以下参数：带通范围，滤波器的阶数，加窗类型。由于心跳频率为0.9-1.6 Hz，所以本论文设计的带通滤波器的的带通范围为1-2 Hz；与高通滤波所加窗函数相同，带通滤波器所加窗函数也为高斯窗，且滤波器的阶数同样设置为50。根据上面的参数设置可以得到下图所示的带通滤波器。