简介

人类心脏的电活动可以被测量为远离基线信号的振幅序列。对于一个正常的心跳周期，心电信号可分为以下几种搏动形态[1]:

ECG波形的这些区域的分割可以为评估人类心脏的整体健康状况和异常的存在提供有用的测量基础[2］．手动标注心电信号的每个区域可能是一项繁琐而耗时的任务。信号处理和深度学习方法可能有助于简化和自动化感兴趣区域注释。

本例使用的ECG信号来自公开可用的QT数据库[3.] [4］．数据包括大约15分钟的ECG记录，采样率为250hz，从总共105名患者中测量。为了获得每个记录，检查人员将两个电极放在患者胸部的不同位置，从而产生双通道信号。数据库提供由自动化专家系统生成的信号区域标签[2］．本例旨在使用深度学习解决方案，根据样本所在的区域为每个心电信号样本提供一个标签。在信号中标记感兴趣区域的过程通常被称为波形分割．

要训练深度神经网络对信号区域进行分类，可以使用长短期记忆(LSTM)网络。这个例子展示了如何使用信号预处理技术和时频分析来提高LSTM分割性能。特别地，该示例使用傅里叶同步压缩变换来表示心电信号的非平稳行为。

下载并准备数据

自动专家系统对105个双通道心电信号中的每个通道进行独立标记，并对其进行独立处理，共210个心电信号与区域标签一起存储在210个mat文件中。文件可在以下位置索取://www.tatmou.com/金宝appsupportfiles/SPT/data/QTDatabaseECGData.zip．

将数据文件下载到临时目录中，该目录的位置由MATLAB®指定tempdir命令。如果要将数据文件放在不同于tempdir，在后续提示中修改目录名称。

%下载数据dataURL =“//www.tatmou.com/金宝appsupportfiles/SPT/data/QTDatabaseECGData1.zip”；datasetFolder = fullfile(tempdir，“QTDataset”）;zipFile = fullfile(tempdir，“QTDatabaseECGData.zip”）;如果~存在(datasetFolder“dir”) websave (zipFile dataURL);解压缩(zipFile tempdir);结束

的解压缩操作创建QTDatabaseECGData文件夹中有210个mat文件。每个文件包含一个心电信号变量ecgSignal以及变量中区域标签的表signalRegionLabels．每个文件还包含信号的采样率变量Fs．在本例中，所有信号的采样率都为250hz。

创建一个信号数据存储以访问文件中的数据。的临时目录中存储了数据集QTDatabaseECGData文件夹中。如果不是这样，请更改下面代码中数据的路径。属性指定要从每个文件中读取的信号变量名称SignalVariableNames参数。

sds = signalDatastore(数据目录，“SignalVariableNames”，[“ecgSignal”，“signalRegionLabels”]）

sds = signalDatastore属性:Files:{'/tmp/QTDataset/ecg1.mat';“/ tmp / QTDataset / ecg10.mat”;“/ tmp / QTDataset / ecg100。垫”……和207更多}AlternateFileSystemRoots: [0×0 string] ReadSize: 1 SignalVariableNames: ["ecgSignal" "signalRegionLabels"]

方法时，数据存储返回一个双元素单元格数组，其中包含一个心电信号和一个区域标签表读函数。使用预览函数可以看到，第一个文件的内容是一个225,000个样本长的心电信号和一个包含3385个区域标签的表。

数据=预览(sds)

data =2×1单元格数组{225000×1 double} {3385×2 table}

查看区域标签表的前几行，观察每一行都包含区域限制索引和区域类值(P、T或QRS)。

{2}(数据)

ans =8×2表ROILimits值__________ _____ 83 117 P 130 153 QRS 201 246 T 285 319 P 332 357 QRS 412 457 T 477 507 P 524 547 QRS

将前1000个样本的标签可视化signalMask对象。

M = signalMask(data{2});plotsigroi (M,数据{1}(1:1000))

通常的机器学习分类程序如下:

该网络使用70%的数据进行训练，并使用剩余的30%进行测试。

为了获得可重复的结果，请重置随机数生成器。使用dividerand函数获取随机索引以打乱文件，而子集的函数signalDatastore将数据划分为训练和测试数据存储。

rng默认的[trainIdx，~， testdx] = diverand(数字(sds.Files)，0.7,0,0.3);trainDs =子集(sds,trainIdx);testDs =子集(sds, testdx);

在这个分割问题中，LSTM网络的输入是一个心电信号，输出是一个与输入信号长度相同的标签序列或掩码。网络任务是用它所属区域的名称标记每个信号样本。因此，有必要将数据集上的区域标签转换为每个信号样本包含一个标签的序列。使用转换后的数据存储和getmask转换区域标签的助手函数。的getmask函数添加标签类别，“n / a”，以标记不属于任何感兴趣区域的样本。

类型getmask.m

函数outputCell = getmask(inputCell) % getmask将区域标签转换为标签的掩码，大小等于输入ECG信号的%大小。inputCell是一个包含心电信号向量%和区域标签表的双元单元格数组。% % outputCell是一个包含心电信号矢量%和与信号长度相同的分类标签矢量掩码的双元单元格数组。The MathWorks, Inc. sig = inputCell{1};roiTable = inputCell{2};L =长度(sig);M = signalMask(roiTable);%获取分类掩码，并在存在重叠掩码= catmask时优先考虑QRS区域(M,L，'OverlapAction'，' priority bylist '，'PriorityList'，[2 1 3]);%设置缺失值为"n/a" mask(ismissing(mask)) = "n/a";outputCell = {sig,mask}; end

预览转换后的数据存储，观察它返回一个信号向量和一个长度相等的标签向量。绘制分类掩码向量的前1000个元素。

trainDs = transform(trainDs， @getmask);testDs = transform(testDs， @getmask);transformdata =预览(训练)

transformedData =1×2单元格数组{224993×1 double} {224993×1 categorical}

情节(transformedData {2} (1:1000))

将非常长的输入信号传递到LSTM网络可能导致估计性能下降和过多的内存使用。为了避免这些影响，可以使用转换后的数据存储和resizeDatahelper函数。helper函数创建尽可能多的5000个样本段，并丢弃剩余的样本。转换后的数据存储输出的预览显示，第一个ECG信号及其标签掩码被分成5000个样本段。注意，转换后的数据存储的预览只显示前8个元素地板(224993/5000)= 44个元素单元格数组，如果我们调用datastore读函数。

trainDs = transform(trainDs，@ resizadata);testDs = transform(testDs，@resizeData);预览(trainDs)

ans =8×2单元格数组{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}{1×5000双}{1×5000分类}

选择训练网络或下载预训练网络

本例的下一部分比较了训练LSTM网络的三种不同方法。由于数据集规模较大，每个网络的训练过程可能需要几分钟。如果您的机器具有GPU和并行计算工具箱™，那么MATLAB将自动使用GPU进行更快的训练。否则，占用CPU。

您可以跳过训练步骤，并使用下面的选择器下载预训练的网络。如果您想在示例运行时训练网络，请选择“train networks”。如果您想跳过训练步骤，请选择“下载网络”和包含所有三个预训练网络的文件rawNet，filteredNet,fsstNet -将下载到您的临时目录，其位置由MATLAB®的tempdir命令。如果您想将下载的文件放在不同于tempdir，在后续提示中修改目录名称。

actionFlag =“列车网络”；如果actionFlag = =“下载网络”下载预先训练好的网络dataURL =“https://ssd.mathworks.com/金宝appsupportfiles/SPT/data/QTDatabaseECGSegmentationNetworks.zip”；% #好< * UNRCH >modelsFolder = fullfile(tempdir，“QTDatabaseECGSegmentationNetworks”）;modelsFile = fullfile(modelsFolder，“trainedNetworks.mat”）;zipFile = fullfile(tempdir，“QTDatabaseECGSegmentationNetworks.zip”）;如果~存在(modelsFolder“dir”) websave (zipFile dataURL);解压缩(zipFile fullfile (tempdir“QTDatabaseECGSegmentationNetworks”));结束加载(modelsFile)结束

下载的网络和新训练的网络之间的结果可能略有不同，因为网络是使用随机初始权重训练的。

将原始心电信号直接输入LSTM网络

首先，使用训练数据集中的原始ECG信号训练LSTM网络。

培训前定义网络架构。指定一个sequenceInputLayer大小为1以接受一维时间序列。属性指定LSTM层“序列”输出模式为每个样本中的信号提供分类。使用200个隐藏节点以获得最佳性能。指定一个fullyConnectedLayer输出大小为4，每个波形类别对应一个。添加一个softmaxLayer和一个classificationLayer输出估计的标签。

层= [.．.sequenceInputLayer (1) lstmLayer (200“OutputMode”，“序列”) fullyConnectedLayer(4) softmaxLayer classificationLayer];

为培训过程选择确保良好网络性能的选项。请参阅trainingOptions(深度学习工具箱)每个参数的描述文档。

选项= trainingOptions(“亚当”，.．.“MaxEpochs”10.．.“MiniBatchSize”, 50岁,.．.“InitialLearnRate”, 0.01,.．.“LearnRateDropPeriod”3,.．.“LearnRateSchedule”，“分段”，.．.“GradientThreshold”, 1.．.“阴谋”，“训练进步”，.．.“洗牌”，“every-epoch”，.．.“详细”0,.．.“DispatchInBackground”,真正的);

因为整个训练数据集都适合内存，所以可以使用高如果parallel Computing Toolbox™可用，则数据存储的功能可以并行转换数据，然后将其收集到工作空间中。神经网络训练是迭代的。在每次迭代中，数据存储从文件中读取数据并在更新网络系数之前转换数据。如果数据适合计算机的内存，则将数据导入工作空间可以实现更快的培训，因为数据只读取和转换一次。注意，如果数据不适合内存，则必须将数据存储传递给训练函数，并且在每个训练阶段执行转换。

为训练集和测试集创建高数组。根据您的系统，MATLAB创建的并行池中的worker数量可能不同。

tallTrainSet =高(trainDs);

使用“本地”配置文件启动并行池(parpool)…连接到并行池(工人数:8)。

tallTestSet = tall(testDs);

现在调用收集tall数组的函数，用于计算整个数据集上的转换，并获得带有训练和测试信号和标签的单元格数组。

trainData = gather(tallTrainSet);

使用并行池“本地”评估tall表达式:-通过1 / 1:在11秒内完成评估，在12秒内完成

: trainData (1)

ans =1×2单元格数组{1×5000 double} {1×5000 categorical}

testData = gather(tallTestSet);

使用并行池“本地”评估tall表达式:-通过1 / 1:在2.9秒内完成评估，在3.1秒内完成

列车网络的

使用trainNetwork命令对LSTM网络进行训练。

如果actionFlag = =“列车网络”rawNet = trainNetwork(trainData(:，1)，trainData(:，2)，layers,options);结束

图中的训练精度和损失子图跟踪所有迭代中的训练进度。使用原始信号数据，该网络正确地将大约77%的样本分类为P波、QRS复合体、T波或未标记区域“n / a”．

测试数据分类

利用训练好的LSTM网络和神经网络对测试数据进行分类分类命令。指定一个小批量大小为50，以匹配培训选项。

(rawNet,testData(:，1))，“MiniBatchSize”, 50);

混淆矩阵提供了一种直观和信息丰富的方法来可视化分类性能。使用confusionchart命令来计算测试数据预测的总体分类精度。对于每个输入，将分类标签的单元格数组转换为行向量。指定列规范化显示，以每个类的样本百分比的形式查看结果。

confusionchart ([predTest{:}]、[testData {: 2}),“归一化”，“column-normalized”）;

使用原始心电信号作为网络的输入，只有约60%的t波样本，40%的p波样本和60%的QRS-complex样本是正确的。为了提高性能，在深度学习网络输入之前应用一些关于ECG信号特征的知识，例如由患者呼吸运动引起的基线漂移。

应用滤波方法去除基线漂移和高频噪声

三种拍子形态占据不同的频带。QRS复合体的频谱通常具有10 - 25hz左右的中心频率，其分量低于40hz。P波和T波出现的频率甚至更低:P波分量低于20hz, T波分量低于10hz [5］．

基线徘徊是由患者呼吸运动引起的低频振荡(< 0.5 Hz)。这种振荡独立于拍的形态，并且不能提供有意义的信息[6］．

设计一个通频带频率范围为[0.5,40]Hz的带通滤波器，以去除漂移和任何高频噪声。删除这些组件可以改进LSTM训练，因为网络不会学习不相关的特征。使用cellfun在高数据单元格数组上并行筛选数据集。

%带通滤波器设计hFilt = designfilt(“bandpassiir”，“StopbandFrequency1”, 0.4215,“PassbandFrequency1”, 0.5,.．.“PassbandFrequency2”现年40岁的“StopbandFrequency2”, 53.345,.．.“StopbandAttenuation1”现年60岁的“PassbandRipple”, 0.1,“StopbandAttenuation2”现年60岁的.．.“SampleRate”, 250,“DesignMethod”，“ellip”）;从转换后的数据存储创建高数组并过滤信号tallTrainSet =高(trainDs);tallTestSet = tall(testDs);filteredTrainSignals = gather(cellfun(@(x)filter(hFilt,x)，tallTrainSet(:，1)，“UniformOutput”、假));

使用并行池“本地”评估高表达式:-通过1:0%完成评估0%完成

-通过1 / 1:在13秒内完成评估

trainLabels = gather(tallTrainSet(:，2));

使用并行池“本地”评估tall表达式:-通过1 / 1:在3.6秒内完成评估，在4秒内完成

filteredTestSignals = gather(cellfun(@(x)filter(hFilt,x)，tallTestSet(:，1)，“UniformOutput”、假));

使用并行池“本地”评估tall表达式:-通过1 / 1:在2.4秒内完成评估，在2.5秒内完成

testLabels = gather(tallTestSet(:，2));

使用并行池“本地”评估tall表达式:-通过1 / 1:在1.9秒内完成评估，在2秒内完成

为一个典型案例绘制原始和过滤后的信号。

trainData = gather(tallTrainSet);

使用并行池“本地”评估tall表达式:-通过1 / 1:在4秒内完成评估，在4.2秒内完成

图subplot(2,1,1) plot(trainData{95,1}(2001:3000))“生”)网格子图(2,1,2)plot(filteredTrainSignals{95}(2001:3000))“过滤”网格)

尽管经过过滤的信号的基线可能会使习惯于在医疗设备上进行传统ECG测量的医生感到困惑，但该网络实际上将受益于漂移去除。

滤波心电信号训练网络

使用与之前相同的网络结构在过滤后的心电信号上训练LSTM网络。

如果actionFlag = =“列车网络”filteredNet = trainNetwork(filteredTrainSignals,trainLabels,layers,options);结束

对信号进行预处理后，训练准确率提高到80%以上。

对过滤后的心电信号进行分类

利用更新后的LSTM网络对预处理后的测试数据进行分类。

predFilteredTest =分类(filteredNet,filteredTestSignals，“MiniBatchSize”, 50);

将分类性能可视化为混淆矩阵。

图confusionchart ([predFilteredTest {:}], [testLabels {}):,“归一化”，“column-normalized”）;

简单预处理可提高t波分类效率约15%，QRS-complex和p波分类效率约10%。

心电信号的时频表示

对时间序列数据进行成功分类的一种常用方法是提取时频特征，并将其输入到网络中，而不是原始数据。然后，该网络同时学习时间和频率的模式[7］．

傅里叶同步压缩变换(FSST)为每个信号样本计算频谱，因此它是理想的分割问题，我们需要保持与原始信号相同的时间分辨率。使用fsst函数检查其中一个训练信号的变换。指定长度为128的Kaiser窗口以提供足够的频率分辨率。

数据=预览(训练);{1 1}图fsst(数据,250年,凯撒(128),“桠溪”）

计算训练数据集中每个信号在感兴趣的频率范围[0.5,40]Hz上的FSST。将FSST的实部和虚部作为单独的特征，并将这两个组件送入网络。此外，通过减去平均值和除以标准差来标准化训练特征。使用转换后的数据存储extractFSSTFeatures辅助函数和高函数并行处理数据。

fsstTrainDs = transform(trainDs，@(x)extractFSSTFeatures(x,250));fsstTallTrainSet = tall(fsstTrainDs);fsstTrainData = gather(fsstTallTrainSet);

使用并行池“本地”评估高表达式:-通过1:0%完成评估0%完成

-通过1 / 1:在2分35秒内完成评估

对测试数据重复此过程。

fsstTTestDs = transform(testDs，@(x)extractFSSTFeatures(x,250));fsstTallTestSet = tall(fsstTTestDs);fsstTestData = gather(fsstTallTestSet);

使用并行池“本地”评估tall表达式:-通过1 / 1:在1分4秒内完成评估

调整网络架构

修改LSTM架构，使网络接受每个样本的频谱，而不是单个值。检查FSST的大小，以查看频率的数量。

大小(fsstTrainData {1})

ans =1×240 5000

指定一个sequenceInputLayer40个输入特征。保持其余网络参数不变。

层= [.．.sequenceInputLayer (40) lstmLayer (200,“OutputMode”，“序列”) fullyConnectedLayer(4) softmaxLayer classificationLayer];

心电信号FSST训练网络

用转换后的数据集训练更新后的LSTM网络。

如果actionFlag = =“列车网络”fsstNet = trainNetwork(fsstTrainData(:，1)，fsstTrainData(:，2)，layers,options);结束

使用时频特征提高了训练精度，现在已经超过90%。

用FSST对测试数据进行分类

利用更新后的LSTM网络和提取的FSST特征，对测试数据进行分类。

predFsstTest = category (fsstNet,fsstTestData(:，1)，“MiniBatchSize”, 50);

将分类性能可视化为混淆矩阵。

confusionchart ([predFsstTest{:}]、[fsstTestData {: 2}),“归一化”，“column-normalized”）;

与原始数据结果相比，使用时频表示可以将t波分类提高约25%，p波分类提高约40%，QRS-complex分类提高30%。

使用一个signalMask对象将网络预测与单个心电信号的地面真实值标签进行比较。忽略了“n / a”标绘感兴趣的区域时的标签。

testData = gather(tall(testDs));

使用并行池“本地”评估高表达式:-通过1:0%完成评估0%完成

-通过1 / 1:在37秒内完成评估

Mtest = signalMask(testData{1,2}(3000:4000));太。speciyselectedcategories = true;太。选择edCategories = find(Mtest.Categories ~=“n / a”）;图subplot(2,1,1) plotsigroi(Mtest,testData{1,1}(3000:4000)) title(“地面实况”) Mpred = signalMask(predFsstTest{1}(3000:4000));mpr。speciyselectedcategories = true;mpr。SelectedCategories = find(Mpred。类别~ =“n / a”）;subplot(2,1,2) plotsigroi(Mpred,testData{1,1}(3000:4000)) title(“预测”）

结论

该实例说明了信号预处理和时频分析如何提高LSTM波形分割性能。带通滤波和基于傅里叶的同步压缩导致所有输出类别的平均改进从55%提高到85%左右。

参考文献

[1] McSharry, Patrick E.，等。生成合成心电图信号的动态模型。IEEE生物医学工程汇刊。第50卷，2003年第3期，第289-294页。

[2]拉古纳，巴勃罗，雷蒙Jané，和Pere Caminal。多导联心电信号中波边界的自动检测:CSE数据库的验证计算机和生物医学研究。第27卷第1期，1994年，第45-60页。

Goldberger, Ary L.， Luis A. N. Amaral, Leon Glass, Jeffery M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody，彭仲康和H. Eugene Stanley。“PhysioBank, PhysioToolkit，和PhysioNet:复杂生理信号新研究资源的组成部分。”循环。卷101，第23期，2000年，第e215-e220页。[流通电子页;http://circ.ahajournals.org/content/101/23/e215.full］．

拉古纳、巴勃罗、罗杰·g·马克、阿里·l·戈德伯格和乔治·b·穆迪。”心电图QT和其他波形间隔测量算法的评估数据库。＂心脏病学中的计算机。Vol.24, 1997, pp. 673-676。

[5] Sörnmo, Leif和Pablo Laguna。“心电图信号处理。”威利生物医学工程百科全书，2006.

[6]科勒，B-U。，Carsten Hennig, and Reinhold Orglmeister. "The principles of software QRS detection."IEEE医学与生物工程杂志．卷二，2002年第1期，页42-57。

萨拉蒙，贾斯汀，还有胡安·巴勃罗·贝罗。“用于环境声音分类的深度卷积神经网络和数据增强。”IEEE信号处理快报。Vol. 24 No. 3, 2017, pp. 279-283。