基于机器学习的频谱态势预测：机器学习和深度学习的心电信号ECG分析与识别

人气：416 ℃/2024-02-25 11:59:29

MATLAB环境下基于小波分析的心电信号ECG分析方法

算法程序运行环境为MATLAB R2018a，主要使用小波变换来分析心电信号，心电信号通常是非平稳的，即频率随时间而变化，由于心电信号的特征通常在时间和频率上局部化，因此在使用诸如小波变换等稀疏表示时会更简单。首先加载并绘制一个心电图波形，其中QRS复合波的R峰已由多名心脏病专家标注，ECG数据和标注取自MIT-BIH心律失常数据库，采样频率为360 Hz。

load mit200figureplot(tm,ecgsig)hold onplot(tm(ann),ecgsig(ann),'ro')xlabel('Seconds')ylabel('Amplitude')title('Subject - MIT-BIH 200')

小波变换将心电信号分解为不同频带的分量，从而实现信号的稀疏表示。由于'sym4' 小波类似于 QRS 复合波，因此选择'sym4' 小波作为分析小波。为了更清楚地说明这一点，提取一个QRS复合波，并用伸缩和平移后的 'sym4' 小波进行比较。

使用最大重叠离散小波变换增强 ECG 波形中的 R峰值，最大重叠离散小波变换相对于离散小波变换是一个未抽取的小波变换，因此可处理任意样本大小，关于最大重叠离散小波变换，我后面会详细的讲解，包括一些改进方法。首先使用默认的'sym4'小波 ECG 波形分解为5层，然后，仅使用第4和第5层的小波系数重建ECG信号，其中

第4层小波系数频率范围为[11.25, 22.5) Hz

第5层小波系数频率范围为[5.625, 11.25) Hz

从小波系数中构建信号近似的平方绝对值，并采用峰值查找算法识别R峰值，因此可以使用matlab自带的findpeaks函数定位波峰

y = abs(y).^2;[qrspeaks,locs] = findpeaks(y,tm,'MinPeakHeight',0.35,... 'MinPeakDistance',0.150);figureplot(tm,y)hold onplot(locs,qrspeaks,'ro')xlabel('Seconds')title('R Peaks Localized by Wavelet Transform with Automatic Annotations')

用小波变换计算的心率为88.60 次/分钟，而标注波形为88.72次/分钟。如果使用原始数据幅值的平方，就会发现小波变换分离R峰的能力：

figureplot(tm,ecgsig,'k--')hold onplot(tm,y,'r','linewidth',1.5)plot(tm,abs(ecgsig).^2,'b')plot(tm(ann),ecgsig(ann),'ro','markerfacecolor',[1 0 0])set(gca,'xlim',[10.2 12])legend('Raw Data','Wavelet Reconstruction','Raw Data Squared', ... 'Location','SouthEast')xlabel('Seconds')

下面看下带噪声的ECG信号

load mit203figureplot(tm,ecgsig)hold onplot(tm(ann),ecgsig(ann),'ro')xlabel('Seconds')ylabel('Amplitude')title('Subject - MIT-BIH 203 with Expert Annotations')

同样使用5层最大重叠离散小波变换分离R峰，使用findpeaks函数来确定峰位的置

看一下第4层和第5层小波系数幅度的平方值以及专家标注的R峰值时间

可以看到 4 层和 5 层小波细节系数的峰值往往同时出现

看完了心电信号，再看一下脑动力学的时变小波相干分析，反正都是生理信号。关于小波相干分析的资料也比较多，大家可自行查阅

有两个受试者的近红外光谱 (NIRS) 数据，NIRS 通过利用含氧血红蛋白和脱氧血红蛋白的不同吸收特性来测量大脑活动，记录部位是两个受试者的上额叶皮层。数据采样频率为10Hz。在实验中，受试者在一项任务上交替进行合作和竞争，任务周期是七秒。

load NIRSDatafigureplot(tm,[NIRSData(:,1) NIRSData(:,2)])legend('Subject 1','Subject 2','Location','NorthWest')xlabel('Seconds')title('NIRS Data')grid on

时域数据中，不清楚时间序列中存在哪些振荡，使用连续小波变换看一下

从CWT中可以看出2个数据集在1 Hz 左右的强调频振荡，这些是由于 2个受试者的心动周期所致。此外，在 0.15 Hz左右，2个数据集都有较弱的振荡。小波相干性可以增强对两个时间序列中共同存在的微弱振荡的检测。

在小波相干性谱图中，在 0.15 Hz 附近有很强的相关性，这代表了2个受试者中与任务相关的大脑活动的相干振荡。添加绘图时间标记，指示2个任务周期，任务之间的时间段是休息时间段

MATLAB环境下基于深层小波时间散射网络的ECG信号分类

数据描述

使用从3种 ECG 数据：心律失常数据、充血性心力衰竭数据和正常窦性心律数据，共使用来自3个 PhysioNet 数据库的162条ECG 记录：MIT-BIH心律失常数据库、MIT-BIH正常窦性心律数据库和BIDMC充血性心力衰竭数据库。共有96个心律失常患者的信号，30个充血性心力衰竭患者的信号，以及36个正常窦性心律患者的信号，目标就是训练分类器来区分心律失常（ARR）、充血性心力衰竭（CHF）和正常窦性心律（NSR）3类信号。

首先加载文件

addpath(genpath(pwd))unzip(fullfile(pwd,'physionet_ECG_data-main.zip'),pwd)unzip(fullfile(pwd,'physionet_ECG_data-main','ECGData.zip'),...fullfile(pwd,'ECGData'))

解压缩 ECGData.zip 文件后，进行数据加载

load(fullfile(pwd,'ECGData','ECGData.mat'))

ECGData是一个结构数组，包含两个字段：Data 和 Labels。数据是一个 162×65536的矩阵，其中每一行是以128Hz采样的ECG信号。每个ECG时间序列

的总持续时间为512秒。标签是一个 162×1 的标签元胞数组，每行1个数据。 3个诊断类别分别是：“ARR”（心律失常）、“CHF”（充血性心力衰竭）和“NSR”（正常窦性心律）。

创建训练和测试集

将数据随机分为训练数据集和测试数据集。将每类中70%的数据随机分配给训练集，剩下的 30%分配给测试集。

得到trainData,testData,trainLabels,testLabels

测试一下每个类别的百分比与数据集中的整体百分比一致

Ctrain = countcats(categorical(trainLabels))./numel(trainLabels).*100

Ctrain = 3×1

59.2920

18.5841

22.1239

Ctest = countcats(categorical(testLabels))./numel(testLabels).*100

Ctest = 3×1

59.1837

18.3673

22.4490

绘制几个样本的波形看看

小波时间散射

在小波时间散射网络中指定的关键参数是the scale of the time invariant，姑且称之为时间不变尺度、小波变换次数以及每个小波滤波器组中每倍频程的小波数量。在许多应用中，两个级联滤波器组足以实现良好的性能。在这个例子中，构建一个带有2个级联滤波器组的小波时间散射网络：第1个滤波器组中每倍频程有8个小波，第2个滤波器组中每倍频程有1个小波，时间不变尺度设置为 150 秒。

可视化2个滤波器组中的小波滤波器

构建小波散射网络后，获取训练数据的散射系数矩阵

scat_features_train = featureMatrix(sn,trainData');

featureMatrix的输出是 409×16×113，张量的每一“页”是一个信号的散射变换。为了得到一个与SVM分类器兼容的矩阵，将多信号散射变换重塑为一个矩阵，其中每一列对应一个散射路径，每一行是一个散射时间窗口。在这种情况下，将获得1808行，因为训练数据113个信号中的每一个都有 16 个时间窗口。

对测试数据重复该过程

创建标签来匹配窗口的数量

[sequence_labels_train,sequence_labels_test] = createSequenceLabels(Nwin,trainLabels,testLabels);

交叉验证，使用5折交叉验证估计错误率

计算损失和混淆矩阵，并显示准确率

predLabels = kfoldPredict(kfoldmodel);loss = kfoldLoss(kfoldmodel)*100;