Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

使用深度学习进行波形分割

此示例说明如何使用递归深度学习网络和时频分析来分割人体心电图 (ECG) 信号。

简介

人类的心电图活动可以作为远离基线信号的幅值序列来测量。对于单个正常心跳周期,ECG 信号可分为以下几种心跳形态 [1]:

  • P 波 - 表示心房去极化的 QRS 复波前的小偏移

  • QRS 复波 - 心跳的最大幅值部分

  • T 波 - 表示心室复极化的 QRS 复波后的小偏移

ECG 波形的这些区域的分割可作为基础测量数据用于评估人类心脏整体健康和异常状况 [2]。手动注释 ECG 信号的每个区域可能是一项乏味且耗时的任务。信号处理和深度学习方法有助于简化注释并自动对关注区域进行注释。

此示例使用来自公开可用的 QT 数据库的 ECG 信号 [3] [4]。这些数据由总共 105 名患者大约 15 分钟的 ECG 录音片段组成,采样率为 250 Hz。为了获得每个录音片段,检查人员将两个电极放置在患者胸部的不同位置,以产生双通道信号。该数据库提供由自动专家系统生成的信号区域标签 [2]。此示例旨在使用深度学习解决方案根据样本所在的区域为每个 ECG 信号样本提供标签。这种对信号中的关注区域加标签的过程通常称为波形分割

为了训练深度神经网络来对信号区域分类,可以使用长短期记忆 (LSTM) 网络。此示例说明如何使用信号预处理方法和时频分析来提高 LSTM 分割性能。具体而言,该示例使用傅里叶同步压缩变换来表示 ECG 信号的非平稳行为。

下载并准备数据

105 个双通道 ECG 信号的每个通道由自动专家系统独立标注并独立处理,总共 210 个 ECG 信号,它们与区域标签一起存储在 210 个 MAT 文件中。这些文件可在以下位置获得:https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData.zip

将数据文件下载到您的临时目录中,临时目录的位置由 MATLAB® 的 tempdir 命令指定。如果要将数据文件放在不同于 tempdir 的文件夹中,请在后续说明中更改目录名称。

% Download the data
dataURL = 'https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData1.zip';
datasetFolder = fullfile(tempdir,'QTDataset');
zipFile = fullfile(tempdir,'QTDatabaseECGData.zip');
if ~exist(datasetFolder,'dir')
     websave(zipFile,dataURL);
     unzip(zipFile,tempdir);
end

unzip 操作会在您的临时目录中创建 QTDatabaseECGData 文件夹,其中包含 210 个 MAT 文件。每个文件在变量 ecgSignal 中包含一个 ECG 信号,在变量 signalRegionLabels 中包含一个区域标签表。每个文件还在变量 Fs 中包含信号的采样率。在此示例中,所有信号的采样率均为 250 Hz。

创建一个信号数据存储来访问文件中的数据。此示例假设数据集已存储在临时目录中的 QTDatabaseECGData 文件夹下。如果不是这样,请更改下面代码中数据的路径。使用 SignalVariableNames 参数指定要从每个文件中读取的信号变量名称。

sds = signalDatastore(datasetFolder,'SignalVariableNames',["ecgSignal","signalRegionLabels"])
sds = 
  signalDatastore with properties:

                       Files:{
                             '/tmp/QTDataset/ecg1.mat';
                             '/tmp/QTDataset/ecg10.mat';
                             '/tmp/QTDataset/ecg100.mat'
                              ... and 207 more
                             }
    AlternateFileSystemRoots: [0×0 string]
                    ReadSize: 1
         SignalVariableNames: ["ecgSignal"    "signalRegionLabels"]

每次调用 read 函数时,数据存储都会返回一个包含 ECG 信号和区域标签表的二元素元胞数组。使用数据存储的 preview 函数,可以看到第一个文件的内容是长度为 225000 个采样的 ECG 信号和一个包含 3385 个区域标签的表。

data = preview(sds)
data=2×1 cell array
    {225000×1 double}
    {  3385×2 table }

查看区域标签表的前几行,注意每行都包含区域范围索引和区域类值(P、T 或 QRS)。

head(data{2})
ans=8×2 table
    ROILimits     Value
    __________    _____

     83    117     P   
    130    153     QRS 
    201    246     T   
    285    319     P   
    332    357     QRS 
    412    457     T   
    477    507     P   
    524    547     QRS 

使用 signalMask 对象可视化前 1000 个采样的标签。

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

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

  1. 将数据库分成训练数据集和测试数据集。

  2. 使用训练数据集训练网络。

  3. 使用经过训练的网络对测试数据集进行预测。

用 70% 的数据对网络进行训练,用其余的 30% 对网络进行测试。

为了获得可重现的结果,请重置随机数生成器。使用 dividerand 函数获得随机索引来对文件进行乱序处理,使用 signalDatastoresubset 函数将数据分成训练数据存储和测试数据存储。

rng default
[trainIdx,~,testIdx] = dividerand(numel(sds.Files),0.7,0,0.3);
trainDs = subset(sds,trainIdx);
testDs = subset(sds,testIdx);

在此分割问题中,LSTM 网络的输入是 ECG 信号,输出是与输入信号长度相同的标签序列或标签掩膜。网络任务将用信号采样点所属区域的名称来标注每个信号采样点。因此,有必要将数据集中的区域标签变换为包含针对每个信号采样点的各个标签的序列。使用变换后的数据存储和 getmask 辅助函数来变换区域标签。getmask 函数添加一个标签类别 "n/a",用于标注不属于任何关注区域的样本。

type getmask.m
function outputCell = getmask(inputCell)
%GETMASK Convert region labels to a mask of labels of size equal to the
%size of the input ECG signal.
%
%   inputCell is a two-element cell array containing an ECG signal vector
%   and a table of region labels. 
%
%   outputCell is a two-element cell array containing the ECG signal vector
%   and a categorical label vector mask of the same length as the signal. 

% Copyright 2020 The MathWorks, Inc.

sig = inputCell{1};
roiTable = inputCell{2};
L = length(sig);
M = signalMask(roiTable);

% Get categorical mask and give priority to QRS regions when there is overlap
mask = catmask(M,L,'OverlapAction','prioritizeByList','PriorityList',[2 1 3]);

% Set missing values to "n/a"
mask(ismissing(mask)) = "n/a";

outputCell = {sig,mask};
end

预览变换后的数据存储,观察它是否返回长度相等的一个信号向量和一个标签向量。绘制分类封装向量的前 1000 个元素。

trainDs = transform(trainDs, @getmask);
testDs = transform(testDs, @getmask);

transformedData = preview(trainDs)
transformedData=1×2 cell array
    {224993×1 double}    {224993×1 categorical}

plot(transformedData{2}(1:1000))

将非常长的输入信号传入 LSTM 网络会导致估计性能下降和内存使用量过多。为了避免这些影响,请使用变换后的数据存储和 resizeData 辅助函数来拆分 ECG 信号及其对应的标签掩膜。辅助函数创建尽可能多的包含 5000 个采样的信号段,并丢弃其余采样。变换后的数据存储的输出预览显示,第一个 ECG 信号及其标签掩膜分成了若干包含 5000 个采样的信号段。请注意,变换后的数据存储的预览仅显示 8 个元素,它们是在我们调用数据存储 read 函数时会生成的包含 floor(224993/5000) = 44 个元素的元胞数组的前 8 个元素。

trainDs = transform(trainDs,@resizeData);
testDs = transform(testDs,@resizeData);
preview(trainDs)
ans=8×2 cell array
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}

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

此示例的下列各节比较三种不同的 LSTM 网络训练方法。由于数据集很大,每个网络的训练过程可能需要几分钟。如果您的机器同时有 GPU 和 Parallel Computing Toolbox™,则 MATLAB 会自动使用 GPU 以加快训练速度。否则将使用 CPU。

您可以跳过训练步骤,使用以下选择器下载预训练网络。如果您要像示例中那样训练网络,请选择“Train Networks”。如果您要跳过训练步骤,请选择“Download Networks”,一个文件(其中包含所有三个预训练网络 - rawNetfilteredNetfsstNet-)将下载到您的临时目录中,其位置由 MATLAB® 的 tempdir 命令指定。如果要将下载的文件放在不同于 tempdir 的文件夹中,请在后续说明中更改目录名称。

actionFlag = "Train networks";
if actionFlag == "Download networks"
    % Download the pre-trained networks
    dataURL = 'https://ssd.mathworks.com/supportfiles/SPT/data/QTDatabaseECGSegmentationNetworks.zip'; %#ok<*UNRCH>
    modelsFolder = fullfile(tempdir,'QTDatabaseECGSegmentationNetworks');
    modelsFile = fullfile(modelsFolder,'trainedNetworks.mat');
    zipFile = fullfile(tempdir,'QTDatabaseECGSegmentationNetworks.zip');
    if ~exist(modelsFolder,'dir')
        websave(zipFile,dataURL);
        unzip(zipFile,fullfile(tempdir,'QTDatabaseECGSegmentationNetworks'));
    end
    load(modelsFile)
end

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

将原始 ECG 信号直接输入 LSTM 网络中

首先,使用来自训练数据集的原始 ECG 信号训练 LSTM 网络。

在训练前定义网络架构。指定大小为 1 的 sequenceInputLayer,以接受一维时序。用 'sequence' 输出模式指定一个 LSTM 层,以便为信号中的每个采样提供分类。使用 200 个隐含节点以获得最佳性能。指定输出大小为 4 的 fullyConnectedLayer,对每个波形类指定一个层。添加一个 softmaxLayer 和一个 classificationLayer 以输出估计的标签。

layers = [ ...
    sequenceInputLayer(1)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

为训练过程选择选项,以确保获得良好的网络性能。有关每个参数的说明,请参阅 trainingOptions 文档。

options = trainingOptions('adam', ...
    'MaxEpochs',10, ...
    'MiniBatchSize',50, ...
    'InitialLearnRate',0.01, ...
    'LearnRateDropPeriod',3, ...
    'LearnRateSchedule','piecewise', ...
    'GradientThreshold',1, ...
    'Plots','training-progress',...
    'shuffle','every-epoch',...
    'Verbose',0,...
    'DispatchInBackground',true);

由于整个训练数据集可放入内存,因此,如果 Parallel Computing Toolbox™ 可用,可以使用数据存储的 tall 函数以并行方式变换数据,然后将其收集到工作区中。神经网络训练是迭代进行的。在每次迭代中,数据存储从文件中读取数据,变换数据,然后更新网络系数。如果数据可放入计算机的内存中,将数据导入工作区可以加快训练速度,因为数据只需读取和变换一次。请注意,如果数据无法放入内存,您必须将数据存储传递给训练函数,并且在每轮训练中执行变换。

为训练集和测试集创建 tall 数组。根据您的系统,MATLAB 创建的并行池中的工作进程数量可能会有所不同。

tallTrainSet = tall(trainDs);
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 8).
tallTestSet = tall(testDs);

现在调用 tall 数组的 gather 函数来计算整个数据集上的变换,并获得具有训练和测试信号及标签的元胞数组。

 trainData = gather(tallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 11 sec
Evaluation completed in 12 sec
 trainData(1,:)
ans=1×2 cell array
    {1×5000 double}    {1×5000 categorical}

 testData = gather(tallTestSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 2.9 sec
Evaluation completed in 3.1 sec

训练网络

使用 trainNetwork 命令训练 LSTM 网络。

if actionFlag == "Train networks"
     rawNet = trainNetwork(trainData(:,1),trainData(:,2),layers,options);
end

图窗中的训练准确度和损失子图跟踪所有迭代的训练进度。使用原始信号数据,网络将大约 77% 的采样正确分类为属于 P 波、QRS 复波、T 波或不带标签的区域 "n/a"

对测试数据分类

使用经过训练的 LSTM 网络和 classify 命令对测试数据进行分类。指定小批量大小为 50 以匹配训练选项。

predTest = classify(rawNet,testData(:,1),'MiniBatchSize',50);

混淆矩阵提供了一种直观的分类性能可视化方式。使用 confusionchart 命令计算用于测试数据预测的总体分类准确度。对于每个输入,将分类标签的元胞数组转换为行向量。指定列归一化显示,以每个类的采样百分比形式查看结果。

confusionchart([predTest{:}],[testData{:,2}],'Normalization','column-normalized');

使用原始 ECG 信号作为网络输入,则只有大约 60% 的 T 波采样、40% 的 T 波采样和 60% 的 QRS 复波采样是正确的。为了提高性能,请在输入到深度学习网络之前应用一些 ECG 信号特性的知识,例如由患者呼吸运动引起的基线漂移。

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

这三种心跳形态占据不同频带。QRS 复波的典型频谱以大约 10-25 Hz 为中心频率,并且其分量低于 40 Hz。发生 P 波和 T 波的频率甚至更低:P 波分量低于 20 Hz,T 波分量低于 10 Hz [5]。

基线漂移是由患者呼吸运动引起的低频 (< 0.5 Hz) 振荡。这种振荡与心跳形态无关,不会提供有意义的信息 [6]。

设计一个通带频率范围为 [0.5, 40] Hz 的带通滤波器,以消除漂移和任何高频噪声。消除这些分量可改进 LSTM 训练,因为网络不会学习不相关特征。对 tall 数据元胞数组使用 cellfun 以并行方式对数据集进行滤波。

% Bandpass filter design
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');

% Create tall arrays from the transformed datastores and filter the signals
tallTrainSet = tall(trainDs);
tallTestSet = tall(testDs);

filteredTrainSignals = gather(cellfun(@(x)filter(hFilt,x),tallTrainSet(:,1),'UniformOutput',false));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 13 sec
Evaluation completed in 14 sec
trainLabels = gather(tallTrainSet(:,2));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 3.6 sec
Evaluation completed in 4 sec
filteredTestSignals = gather(cellfun(@(x)filter(hFilt,x),tallTestSet(:,1),'UniformOutput',false));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 2.4 sec
Evaluation completed in 2.5 sec
testLabels = gather(tallTestSet(:,2));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1.9 sec
Evaluation completed in 2 sec

对一种典型情况下的原始信号和滤波后的信号绘图。

trainData = gather(tallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 4 sec
Evaluation completed in 4.2 sec
figure
subplot(2,1,1)
plot(trainData{95,1}(2001:3000))
title('Raw')
grid
subplot(2,1,2)
plot(filteredTrainSignals{95}(2001:3000))
title('Filtered')
grid

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

使用滤波后的 ECG 信号训练网络

使用与以前相同的网络架构基于滤波后的 ECG 信号训练 LSTM 网络。

if actionFlag == "Train networks"
    filteredNet = trainNetwork(filteredTrainSignals,trainLabels,layers,options);
end

信号预处理将训练准确度提高到 80% 以上。

对滤波后的 ECG 信号进行分类

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

predFilteredTest = classify(filteredNet,filteredTestSignals,'MiniBatchSize',50);

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

figure
confusionchart([predFilteredTest{:}],[testLabels{:}],'Normalization','column-normalized');

简单的预处理将 T 波分类提高了约 15%,将 QRS 复波和 P 波分类提高了约 10%。

ECG 信号的时频表示

时序数据成功分类的常见方法是提取时频特征并将其馈送到网络而不是原始数据。然后,网络同时跨时间和频率学习模式 [7]。

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

data =  preview(trainDs);
figure
fsst(data{1,1},250,kaiser(128),'yaxis')

基于关注的频率范围 [0.5, 40] Hz 计算训练数据集中每个信号的 FSST。将 FSST 的实部和虚部视为单独的特征,并将两个分量都馈送到网络中。而且,通过减去均值并除以标准差来标准化训练特征。使用变换后的数据存储、extractFSSTFeatures 辅助函数和 tall 函数来并行处理数据。

fsstTrainDs = transform(trainDs,@(x)extractFSSTFeatures(x,250));
fsstTallTrainSet = tall(fsstTrainDs);
fsstTrainData = gather(fsstTallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 2 min 35 sec
Evaluation completed in 2 min 35 sec

对测试数据重复此过程。

fsstTTestDs = transform(testDs,@(x)extractFSSTFeatures(x,250));
fsstTallTestSet = tall(fsstTTestDs);
fsstTestData = gather(fsstTallTestSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1 min 4 sec
Evaluation completed in 1 min 4 sec

调整网络架构

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

size(fsstTrainData{1,1})
ans = 1×2

          40        5000

指定一个包含 40 个输入特征的 sequenceInputLayer。保持其余网络参数不变。

layers = [ ...
    sequenceInputLayer(40)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

使用 ECG 信号的 FSST 训练网络

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

if actionFlag == "Train networks"
    fsstNet = trainNetwork(fsstTrainData(:,1),fsstTrainData(:,2),layers,options);
end

使用时频特征提高了训练准确度,现在已超过 90%。

用 FSST 对测试数据进行分类

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

predFsstTest = classify(fsstNet,fsstTestData(:,1),'MiniBatchSize',50);

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

confusionchart([predFsstTest{:}],[fsstTestData{:,2}],'Normalization','column-normalized');

与原始数据结果相比,使用时间频率表示法将 T 波分类提高了约 25%,将 P 波分类提高了约 40%,将 QRS 复波分类提高了 30%。

使用 signalMask 对象将网络预测与单个 ECG 信号的真实值标签进行比较。绘制关注区域时忽略 "n/a" 标签。

testData = gather(tall(testDs));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 37 sec
Evaluation completed in 37 sec
Mtest = signalMask(testData{1,2}(3000:4000));
Mtest.SpecifySelectedCategories = true;
Mtest.SelectedCategories = find(Mtest.Categories ~= "n/a");

figure
subplot(2,1,1)
plotsigroi(Mtest,testData{1,1}(3000:4000))
title('Ground Truth')

Mpred = signalMask(predFsstTest{1}(3000:4000));
Mpred.SpecifySelectedCategories = true;
Mpred.SelectedCategories = find(Mpred.Categories ~= "n/a");

subplot(2,1,2)
plotsigroi(Mpred,testData{1,1}(3000:4000))
title('Predicted')

结论

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

参考资料

[1] McSharry, Patrick E., et al."A dynamical model for generating synthetic electrocardiogram signals."IEEE® Transactions on Biomedical Engineering.Vol. 50, No. 3, 2003, pp. 289–294.

[2] Laguna, Pablo, Raimon Jané, and Pere Caminal."Automatic detection of wave boundaries in multilead ECG signals:Validation with the CSE database."Computers and Biomedical Research.Vol. 27, No. 1, 1994, pp. 45–60.

[3] Goldberger, Ary L., Luis A. N. Amaral, Leon Glass, Jeffery M. Hausdorff, Plamen Ch.Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-Kang Peng, and H. Eugene Stanley."PhysioBank, PhysioToolkit, and PhysioNet:Components of a New Research Resource for Complex Physiologic Signals."Circulation.Vol. 101, No. 23, 2000, pp. e215–e220. [Circulation Electronic Pages; http://circ.ahajournals.org/content/101/23/e215.full].

[4] Laguna, Pablo, Roger G. Mark, Ary L. Goldberger, and George B. Moody."A Database for Evaluation of Algorithms for Measurement of QT and Other Waveform Intervals in the ECG."Computers in Cardiology.Vol.24, 1997, pp. 673–676.

[5] Sörnmo, Leif, and Pablo Laguna."Electrocardiogram (ECG) signal processing."Wiley Encyclopedia of Biomedical Engineering, 2006.

[6] Kohler, B-U., Carsten Hennig, and Reinhold Orglmeister."The principles of software QRS detection."IEEE Engineering in Medicine and Biology Magazine.Vol. 21, No. 1, 2002, pp. 42–57.

[7] Salamon, Justin, and Juan Pablo Bello."Deep convolutional neural networks and data augmentation for environmental sound classification."IEEE Signal Processing Letters.Vol. 24, No. 3, 2017, pp. 279–283.

另请参阅

| | | | (Signal Processing Toolbox) | (Signal Processing Toolbox)

相关主题