基于小波分析的打鼾(阻塞性睡眠呼吸暂停)检测

使用心冲击信号进行打鼾检测，可用于智能穿戴设备，共3大步骤：1. 探索性数据分析； 2.统计分析；3.打鼾(阻塞性睡眠呼吸暂停)检测

首先导入相关模块

import sys
import warnings
import itertools
warnings.filterwarnings("ignore")
### 数据分析相关模块
import math
from scipy import fftpack,signal
from mne.time_frequency import morlet # 创建Morlet小波import pandas as pd
import numpy as np###统计分析相关模块
from scipy.stats import levene, shapiro, f_oneway
from statsmodels.stats.multicomp import pairwise_tukeyhsd, MultiComparison
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm#pip install qgrid
import qgrid # 可视化相关模块
import matplotlib.pyplot as plt
import seaborn as sns

导入数据

test1 = pd.read_csv('sample1.csv', header = None, names = ['Amplitude'])
test2 = pd.read_csv('sample2.csv', header = None, names = ['Amplitude'])

对test1.csv通道数据进行分析，首先看下波形

f, ax = plt.subplots(figsize=(20, 6))
ax.set_xlim(0,1500)
sns.lineplot(x=test1.index, y=test1.Amplitude, data=test1)

从上面的波形可以观察到，在约 3 秒和约 4 秒时显示出打鼾的迹象（1s有250个点），在打鼾期间，波峰和波谷偏离正常模式。基于以上两点，可以假设在打鼾期间，振动幅值的差异会较大。每一秒的波形都细看一下

fig, axes = plt.subplots(2,3, sharey=False, sharex=False)
fig.set_figwidth(20)
fig.set_figheight(15)
axes[0][0].plot(test1.index[0:249], test1.Amplitude[0:249], label='1st Second')
axes[0][0].plot(test1.index[0:249], test1.Amplitude[0:249].diff(), label='Differenced Seconds')
axes[0][0].legend(loc='best')
axes[0][1].plot(test1.index[250:499], test1.Amplitude[250:499], label='2nd Second')
axes[0][1].plot(test1.index[250:499], test1.Amplitude[250:499].diff(), label='Differenced Seconds')
axes[0][1].legend(loc='best')
axes[0][2].plot(test1.index[500:749], test1.Amplitude[500:749], label='3rd Second')
axes[0][2].plot(test1.index[500:749], test1.Amplitude[500:749].diff(), label='Differenced Seconds')
axes[0][2].legend(loc='best')
axes[1][0].plot(test1.index[750:999], test1.Amplitude[750:999], label='4th Second')
axes[1][0].plot(test1.index[750:999], test1.Amplitude[750:999].diff(), label='Differenced Seconds')
axes[1][0].legend(loc='best')
axes[1][1].plot(test1.index[1000:1249], test1.Amplitude[1000:1249], label='5th Second')
axes[1][1].plot(test1.index[1000:1249], test1.Amplitude[1000:1249].diff(), label='Differenced Seconds')
axes[1][1].legend(loc='best')
axes[1][2].plot(test1.index[1250:1499], test1.Amplitude[1250:1499], label='6th Second')
axes[1][2].plot(test1.index[1250:1499], test1.Amplitude[1250:1499].diff(), label='Differenced Seconds')
axes[1][2].legend(loc='best')
plt.tight_layout()
plt.show()

从上图可以看出，打鼾区域的差分序列的幅值波动较大，非打鼾区域的差分序列幅值波动较小。差分序列怎么算我就不多说了，然后我们在 60 秒的窗口中查看一下差分序列

从上图可以看出，打鼾的周期性是显而易见的，差分序列的局部波形和小波的波形极其相似。

下面开始进行统计分析

test1.describe()

样本平均值为 51.34 次振动/秒，样本最小值 为 20 次振动/秒，样本最大值 为73 次振动/秒。

开始进行特征工程

def create_features(df):"""1. Avg_Amplitude 2. Min_Amplitude 3. Max_Amplitude 4. StDev_Amplitude5. Seconds6. Minute"""df2 = pd.DataFrame(data=None)array = df.Amplitude[0:59750].values.reshape(-1,250) df2['Avg_Amplitude'] = [np.mean(i) for i in array]df2['Min_Amplitude'] = [np.min(i) for i in array]df2['Max_Amplitude'] = [np.max(i) for i in array]df2['StDev_Amplitude'] = [np.std(i) for i in array]df2['Seconds'] = np.tile([int(i) for i in range(1,61)], array.shape[0])[0:array.shape[0]]df2['Minute'] = np.repeat([minute for minute in ['1st', '2nd', '3rd','4th']], 60)[0:array.shape[0]]return df2data_for_analysis_1 = create_features(df=test1)

绘制平均幅值

plt.figure(figsize=(15, 8))
plt.plot(data_for_analysis_1['Avg_Amplitude'])
plt.xlabel('Time')
plt.ylabel('Average Amplitude')

平均每秒幅值看起来没有趋势。绘制每分钟每秒的平均幅值

seconds_based_amplitude = pd.pivot_table(data_for_analysis_1, values='Avg_Amplitude',\columns='Minute', index='Seconds')
seconds_based_amplitude.plot(figsize=(25,10))

第一分钟的尖峰幅值更加明显，绘制每 15 秒间隔的平均幅值

min_winse_variable = pd.pivot_table(data_for_analysis_1, values='Avg_Amplitude',\columns='Seconds', index='Minute')
#前15秒
min_winse_variable.loc[:, 1:15].plot(figsize=(25,10))

接下来的15秒

min_winse_variable.loc[:, 16:30].plot(figsize=(25,15))

接下来15秒

min_winse_variable.loc[:, 31:45].plot(figsize=(25,15))

最后15秒

min_winse_variable.loc[:, 46:60].plot(figsize=(25,15))

综合起来

min_winse_variable.plot(figsize=(25,15))

上图观察到第 1 分钟和第 4 分钟的平均幅值变化很大。基于以上繁琐的可视化分析，可以产生几个假设：每分钟平均幅值不变；每分钟平均幅值的方差是均匀的。

data_for_analysis_1['Seconds'] = data_for_analysis_1.Seconds.astype(str)

进行统计学上的Shapiro's测试以确认正态性假设

def shapiro_test(df):"""Function to Conduct Shapiro Test.Ho : Data is Normally Distributed.Ha : Data is not normally distributed."""for val in np.unique(df['Minute']):p_val = shapiro(df.query('Minute == "{0}"'.format(val))['Avg_Amplitude'])[1]if p_val > 0.05:print('Data is Normally distributed for {0} Minute'.format(val))else:print('Data is Not Normally distributed for {0} Minute'.format(val))shapiro_test(data_for_analysis_1)

假设每分钟的平均幅值具有相等的方差

levene_pval = levene(data_for_analysis_1.query('Minute == "1st"')['Avg_Amplitude'],\data_for_analysis_1.query('Minute == "2nd"')['Avg_Amplitude'],\data_for_analysis_1.query('Minute == "3rd"')['Avg_Amplitude'],\data_for_analysis_1.query('Minute == "4th"')['Avg_Amplitude'] )[1]if levene_pval > 0.05:print('我们未能拒绝方差相等的零假设')else:print('证据不足，无法得出方差相等的结论')

我们未能拒绝方差相等的零假设。由于满足齐次方差的条件，可以进行方差分析

anova_model = ols('Avg_Amplitude ~ C(Minute)', data=data_for_analysis_1).fit()aov_table = anova_lm(anova_model, type=1)print(aov_table)

我擦，统计分析是真特么复杂。

最后，开始打鼾检测流程， 定义采样频率和时间向量

sf = 250
times = np.arange(test1.Amplitude.size)/sf

接下来，使用 Morlet 小波对信号进行卷积，至于为什么会用小波处理，也就不用我多说了吧，设置小波参数，并绘图

sf = 250
times = np.arange(test1.Amplitude.size)/sf#参数
cf = 13     # 中心频率（Hz）
nc = 12     # 振荡次数# 计算小波
wlt = morlet(sf, [cf], n_cycles=nc)[0]#绘图
t = np.arange(wlt.size) / sf
fig, ax = plt.subplots(1, 1, figsize=(7, 4))
ax.plot(t, wlt)
plt.ylim(-0.4, 0.4)
plt.xlim(t[0], t[-1])
plt.xlabel('Time')
plt.ylabel('Amplitude')

信号与小波卷积，并提取幅度和相位

def Wavelet_Formation(df, threshold, window_1, window_2):df['FirstOrderDiff'] = df['Amplitude'].diff().bfill()analytic = np.convolve(df.FirstOrderDiff, wlt, mode='same')magnitude = np.abs(analytic)phase = np.angle(analytic)# Square and normalize the magnitude from 0 to 1 (using the min and max)# Square and normalize the magnitude from 0 to 1 (using the min and max)power = np.square(magnitude)norm_power = (power - power.min()) / (power.max() - power.min())# Define the thresholdthresh = threshold# Find supra-threshold valuessupra_thresh = np.where(norm_power >= thresh)[0]# Create vector for plotting purposesval_spindles = np.nan * np.zeros(df.FirstOrderDiff.size)val_spindles[supra_thresh] = df.Amplitude[supra_thresh]# Plotfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(20, 20), sharex=True)ax1.plot(df.index, df.Amplitude, lw=1.5, label = 'Actual Value')ax1.plot(df.index, df.FirstOrderDiff, lw=1.5, label= 'Differenced Series')ax1.plot(df.index, val_spindles, color='red', alpha=.8, label='Detected Anomaly') ## Red Region marks the detected Snoring Region.ax1.set_xlim(window_1, window_2)ax1.set_ylabel('Amplitude')ax1.set_title('Comparison Graph')ax1.legend(loc='best')ax2.plot(df.index,norm_power)ax2.set_xlabel('Time-stamps')ax2.set_ylabel('Normalized wavelet power')ax2.axhline(thresh, ls='--', color='indianred', label='Threshold')ax2.fill_between(df.index, norm_power, thresh, where = norm_power >= thresh,color='indianred', alpha=.8)plt.legend(loc='best')

查看前 4000 个时间点。

Wavelet_Formation(df=test1,threshold=0.015, window_1=0, window_2=4000)

查看前 6000 个时间点。

Wavelet_Formation(df= test1,threshold=0.015, window_1=4000, window_2=10000)

基于差分项的均方误差

详细代码

基于小波分析的打鼾(阻塞性睡眠呼吸暂停)检测相关推荐

一种基于深度学习的单导联心电信号睡眠呼吸暂停检测方法
在R峰识别的基础上,加入S峰的识别,并论正了该策略对检测结果的有效性. 1.大致方法将数据集(ECG信号)划分为每五分钟的一个片段,为了减少噪声和信号伪影,首先对信号应用了一个有限脉冲响应(FIR) ...
# 睡眠3秒_小儿睡眠呼吸暂停综合征
生活中,不少家长看到孩子睡觉打呼噜会以为孩子累了睡得香,睡得沉.事实上,很有可能你的宝宝是"睡眠呼吸暂停综合征"的患者.作为家长学会基本辨别,发现异常要及时带孩子到医院就诊避免影响 ...
文献分享基于ECG和ECG呼吸信号特征的阻塞性和限制性肺病自动检测
原文链接https://www.sciencedirect.com/science/article/pii/S1746809421003888?via%3Dihub 摘要: 早期发现阻塞性和限制性肺病 ...
在数字乳腺X照片中基于小波分析和统计分析的微钙化检测新特征
摘要:计算机辅助诊断(CAD)系统在乳腺癌的早期诊断中具有重要作用.在本文中,基于对微钙化检测的新特征集,提出了一套CAD系统.新特征是受几个对一些经典特征(高阶统计特征HOS.离散小波变换DWT和小 ...
监听端口的非阻塞性不具有继承性
监听socket设置为非阻塞时,在接收到连接请求时,连接socket并没有将阻塞性继承下来.以具体的代码为例: #include <sys/socket.h> #include < ...
响应式web（四）：使用Netty作为web容器，基于注解的WebFlux阻塞式与响应式实现
目录使用 WebFlux,针对IO密集度比较高的系统,性能会有提升. 使用 Netty 作为 web 容器:注释掉spring-boot-starter-web,启动就默认用的 netty 而不是 ...
python小波分析法检测火焰_一种基于小波分析的网络流量异常检测方法
一种基于小波分析的网络流量异常检测方法杜臻 ; 马立鹏 ; 孙国梓 [期刊名称] <计算机科学> [年 ( 卷 ), 期] 2019(046)008 [摘要] 对大量网络流量数据进行高质 ...
汽车平顺性与仿真分析matlab,基于matlab的汽车平顺性的建模与仿真.docx
基于matlab的汽车平顺性的建模与仿真.docx (1)基于MATLAB的汽车平顺性的建模与仿真车辆工程专硕1601Z1604050李晨1数学建模过程11建立系统微分方程如下图所示,为车身与车轮二自 ...
基于小波分析和机器学习的时间序列分析与识别
研究对象:ECG等时间序列信号方法:小波变换,简单神经网络首先导入相关模块,需要安装尺度谱模块:pip install scaleogram 和mat4py模块:pip install mat4p ...

基于小波分析的打鼾(阻塞性睡眠呼吸暂停)检测

基于小波分析的打鼾(阻塞性睡眠呼吸暂停)检测相关推荐

最新文章

热门文章