基于小波分析的打鼾(阻塞性睡眠呼吸暂停)检测
发布时间
阅读量：
阅读量
通过心冲击信号对打鼾情况进行监测，在智能可穿戴设备中实现对该现象的实时监测与记录。该系统包含以下三个主要步骤：第一步为数据预处理阶段；第二步为特征分析阶段；第三步为阻塞性睡眠呼吸暂停事件的识别。
首先导入相关模块
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_Amplitude
5. Seconds
6. 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 df2
data_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 threshold
thresh = threshold
# Find supra-threshold values
supra_thresh = np.where(norm_power >= thresh)[0]
# Create vector for plotting purposes
val_spindles = np.nan * np.zeros(df.FirstOrderDiff.size)
val_spindles[supra_thresh] = df.Amplitude[supra_thresh]
# Plot
fig, (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)
基于差分项的均方误差
详细代码
知乎专栏
以小波分析为基础的研究对打鼾（即阻塞性睡眠呼吸暂停）进行检测具有重要意义
全部评论 (0)
还没有任何评论哟~