使用深度神经网络对心电图ECG 信号进行分类

具体的软硬件实现点击MCU-AI技术网页_MCU-AI人工智能

此示例说明如何使用深度学习和信号处理对来自 PhysioNet 2017 Challenge 的心电图 (ECG) 数据进行分类。具体而言，该示例使用长短期记忆网络和时频分析。

ECG 记录一段时间内人体心脏的电活动。医生使用 ECG 检测患者的心跳是正常还是不规则。心房颤动 (AFib) 是一种不规则的心跳，当心脏的上腔（即心房）与下腔（即心室）失去协调时就会发生心房颤动。

此示例使用的 ECG 数据来自 PhysioNet 2017 Challenge [1]、[2]、[3]，可从 AF Classification from a Short Single Lead ECG Recording: The PhysioNet/Computing in Cardiology Challenge 2017 v1.0.0 获得。数据由一组以 300 Hz 采样的 ECG 信号组成，由一组专家分成四个不同类：正常 (N)、AFib (A)、其他心律 (O) 和含噪记录 (~)。此示例说明如何使用深度学习来自动化分类过程。该过程使用一个二类分类器，该二类分类器可以将正常 ECG 信号和显示 AFib 符号的信号区分开来。

长短期记忆 (LSTM) 网络是一种循环神经网络 (RNN)，非常适合研究序列和时间序列数据。LSTM 网络可以学习序列的时间步之间的长期相关性。LSTM 层 (lstmLayer) 可以前向分析时间序列，而双向 LSTM 层 (bilstmLayer) 可以前向和后向分析时间序列。此示例使用双向 LSTM 层。

     此示例说明在解决人工智能 (AI) 问题时使用以数据为中心的方法的优势。最初尝试使用原始数据训练 LSTM 网络得到的结果未达到标准。使用提取的特征训练相同的模型架构使分类性能得到显著提高。实例为MATLAB语言。

运行 ReadPhysionetData 脚本，从 PhysioNet 网站下载数据，并生成包含适当格式的 ECG 信号的 MAT 文件 (PhysionetData.mat)。下载数据可能需要几分钟。使用一个条件语句，限定仅在当前文件夹中不存在 PhysionetData.mat 时才运行该脚本。

if ~isfile('PhysionetData.mat')

ReadPhysionetData         

end

load PhysionetData

加载操作会向工作区添加两个变量：Signals 和 Labels。Signals 是保存 ECG 信号的元胞数组。Labels 是分类数组，它保存信号的对应真实值标签。

可视化每个类中一个信号的一段。AFib 心跳间隔不规则，而正常心跳会周期性发生。AFib 心跳信号还经常缺失 P 波，P 波在正常心跳信号的 QRS 复合波之前出现。正常信号的绘图会显示 P 波和 QRS 复合波。

要设计分类器，请使用上一节中生成的原始信号。将信号分成一个训练集（用于训练分类器）和一个测试集（用于基于新数据测试分类器的准确度）。使用 summary 函数显示 AFib 信号与正常信号的比率为 718:4937，约为 1:7。

由于约 7/8 的信号是正常信号，因此分类器会发现通过简单地将所有信号分类为正常信号就可达到高准确度。为了避免这种偏置，需要通过复制数据集中的 AFib 信号来增加 AFib 数据，以便正常信号和 AFib 信号的数量相同。这种复制通常称为过采样，是深度学习中使用的一种数据增强形式。

afibX = Signals(Labels=='A');
afibY = Labels(Labels=='A');

normalX = Signals(Labels=='N');
normalY = Labels(Labels=='N');

  接下来，使用 dividerand 将每个类的目标随机分为训练集和测试集。
[trainIndA,~,testIndA] = dividerand(718,0.9,0.0,0.1);
[trainIndN,~,testIndN] = dividerand(4937,0.9,0.0,0.1);

XTrainA = afibX(trainIndA);
YTrainA = afibY(trainIndA);

XTrainN = normalX(trainIndN);
YTrainN = normalY(trainIndN);

XTestA = afibX(testIndA);
YTestA = afibY(testIndA);

XTestN = normalX(testIndN);
YTestN = normalY(testIndN);

   现在有 646 个 AFib 信号和 4443 个正常信号用于训练。要在每个类中获得相同数量的信号，请使用前 4438 个正常信号，然后
使用 repmat 对前 634 个 AFib 信号重复七次。对于测试集，现在有 72 个 AFib 信号和 494 个正常信号。使用前 490 个正常
信号，然后使用 repmat 对前 70 个 AFib 信号重复七次。默认情况下，神经网络会在训练前随机对数据进行乱序处理，以确保相
邻信号不都有相同的标签。

XTrain = [repmat(XTrainA(1:634),7,1); XTrainN(1:4438)];
YTrain = [repmat(YTrainA(1:634),7,1); YTrainN(1:4438)];

XTest = [repmat(XTestA(1:70),7,1); XTestN(1:490)];
YTest = [repmat(YTestA(1:70),7,1); YTestN(1:490);];

  现在，正常信号和 AFib 信号在训练集和测试集中的分布均衡。
  
     LSTM 网络可以学习序列数据的时间步之间的长期相关性。此示例使用双向 LSTM 层 bilstmLayer，因为它前向和后向检测序列。
  由于输入信号各有一个维度，将输入大小指定是大小为 1 的序列。指定输出大小为 100 的一个双向 LSTM 层，并输出序列的最后
  一个元素。此命令指示双向 LSTM 层将输入时间序列映射到 100 个特征，然后为全连接层准备输出。最后，通过包含大小为 2 
  的全连接层，后跟 softmax 层和分类层，来指定两个类。
     layers = [ ...
    sequenceInputLayer(1)
    bilstmLayer(100,'OutputMode','last')
    fullyConnectedLayer(2)
    softmaxLayer
    classificationLayer
    ]
    
       接下来指定分类器的训练选项。将 'MaxEpochs' 设置为 10，以允许基于训练数据对网络进行 10 轮训练。'MiniBatchSize' 
    为 150 指示网络一次分析 150 个训练信号。'InitialLearnRate' 为 0.01 有助于加快训练过程。指定 'SequenceLength' 
    为 1000 以将信号分解成更小的片段，这样机器就不会因为一次处理太多数据而耗尽内存。将 'GradientThreshold' 设置为 1 
    以防止梯度过大，从而稳定训练过程。将 'Plots' 指定为 'training-progress'，以生成显示训练随迭代次数的增加而变化的
    进度图。将 'Verbose' 设置为 false 以隐藏对应于图中所示数据的表输出。如果您要查看此表，请将 'Verbose' 设置为 
    true。此示例使用自适应矩估计 (ADAM) 求解器。与默认的具有动量的随机梯度下降 (SGDM) 求解器相比，ADAM 在使用 LSTM 
    之类的 RNN 时性能更好。
    
    options = trainingOptions('adam', ...
    'MaxEpochs',10, ...
    'MiniBatchSize', 150, ...
    'InitialLearnRate', 0.01, ...
    'SequenceLength', 1000, ...
    'GradientThreshold', 1, ...
    'ExecutionEnvironment',"auto",...
    'plots','training-progress', ...
    'Verbose',false);
    
    现在用相同的网络对测试数据进行分类。
    testPred = classify(net,XTest,'SequenceLength',1000);
    计算测试准确度。
    LSTMAccuracy = sum(testPred == YTest)/numel(YTest)*100
    LSTMAccuracy = 81.2245
    
    从数据中提取特征有助于提高分类器的训练和测试准确度。本示例采用的方法是计算瞬时频率 (instfreq)、谱熵 (pentropy)。
    可视化每个信号类型的瞬时频率。

可视化每个信号类型的谱熵。

   现在每个信号都有两个维度，就有必要通过将输入序列大小指定为 2 来修改网络架构。指定输出大小为 100 的一个双向

 LSTM 层，并输出序列的最后一个元素。通过使用一个大小为 2 的全连接层，后跟 softmax 层和分类层，来指定两个类。

layers = [ ...
    sequenceInputLayer(2)
    bilstmLayer(100,'OutputMode','last')
    fullyConnectedLayer(2)
    softmaxLayer
    classificationLayer
    ]
    
    指定训练选项。将最大轮数设置为 30，以允许基于训练数据对网络进行 30 轮训练。

  

options = trainingOptions('adam', ...
    'MaxEpochs',30, ...
    'MiniBatchSize', 150, ...
    'InitialLearnRate', 0.01, ...
    'GradientThreshold', 1, ...
    'ExecutionEnvironment',"auto",...
    'plots','training-progress', ...
    'Verbose',false);

      使用更新后的网络对测试数据进行分类。

  testPred2 = classify(net2,XTestSD);

  LSTMAccuracy = sum(testPred2 == YTest)/numel(YTest)*100
  
  LSTMAccuracy = 92.1020