Engee 文档
Notebook

心律失常的光谱分析

该示例演示了频谱分析在具有不同类型心率的心电图信号中的应用。 考虑了从开放PhysioNet数据库获得的具有正常窦性心律、房外和心房颤动的ECG记录。 EngeeDSP库的功能用于频率分析。

功能 engee.clear() 清理工作区:

In [ ]: 
engee.clear()

两种图形渲染模式可用于可视化结果。:

*对于静态，快速渲染，使用 gr().

*对于交互式，具有缩放和查看值的能力的动态可视化，使用 plotlyjs().

通过注释掉第二个命令来选择其中一个命令。 两个后端是互斥的，最后一个被调用的后端将被激活。

In [ ]: 
# gr()
plotlyjs()
Out[0]:
Plots.PlotlyJSBackend()

我们将使用函数进行连接 include 的"read"文件。jl"用于读取文件:

In [ ]: 
include("$(@__DIR__)/read.jl")
Out[0]:
safe_read_ann (generic function with 1 method)

心电图的频谱分析是正常的

与心脏的正常窦性节律相对应的心电图信号用作所述分析的初始数据。 该记录是从[MIT-BIH正常窦性节律数据库（NSRDB）]获得的（https://physionet.org/content/nsrdb/1.0.0 /）并反映了一个健康人的心脏的工作，没有明显的节律紊乱。

使用选定的ECG记录指定文件的路径:

In [ ]: 
data = ("$(@__DIR__)/mit-bih-normal-sinus-rhythm/16272")
Out[0]:
"/user/DemoPublic/biomedical/spectral_analysis_arrhythmia/mit-bih-normal-sinus-rhythm/16272"

设置文件路径后，读取心电记录头。 从报头中提取描述信号参数的服务信息：记录的名称、采样频率、通道数和采样总数。 还显示有关记录的数据格式、增益和文本描述的信息。

In [ ]: 
hdr = read_header(data)

println("记录名称:   ", hdr.record)
println("采样率: ", hdr.fs, " 赫兹")
println("信号数量:   ", hdr.nsig)
println("计数数目:   ", hdr.nsamp)
println("数据格式:         ", hdr.fmt)
println("增益:              ", hdr.gain)
println("资料描述:              ", hdr.desc)

Наименование записи:   16272
Частота дискретизации: 128.0 Гц
Количество сигналов:   2
Количество отсчётов:   11520000
Формат данных:         212
Усиление:              0.0
Описание:              ECG1

为了演示心电图的时间表示，我们将选择一个10秒的记录片段。

In [ ]: 
t_start = 0.0
t_dur   = 10.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end   = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y1 = phys === nothing ? float.(raw[:,1])/200 : phys[:,1]

t1 = (0:length(y1)-1)./ hdr.fs.+ t_start

plot(t1, y1,
        xlabel="时间，从", 
        ylabel="的振幅", 
        title="录录玫搂鲁拢潞(hdr.录)NSRDB",
        legend=false)
Out[0]:

我们还将演示心电图记录的一分钟片段。 为此，我们将分析间隔的持续时间增加到60秒，计算信号的相应部分并构建其时间表示。

In [ ]: 
t_start = 0.0
t_dur = 60.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y1 = phys === nothing ? float.(raw[:,1])/200 : phys[:,1]

t1 = (0:length(y1)-1) ./ hdr.fs .+ t_start

plot(t1, y1,
        xlabel="时间，从", 
        ylabel="的振幅", 
        title="录录玫搂鲁拢潞(hdr.录)NSRDB",
        legend=false)
Out[0]:

要切换到信号的频率表示，我们使用EngeeDSP库的功能。 功能 EngeeDSP.Functions.fft 对时间信号执行快速傅立叶变换，将其转换为一组复数频率系数。 然后它被应用 EngeeDSP.Functions.fftshift 用于对光谱进行定心。

In [ ]: 
nsampl = length(y1)
fs = hdr.fs
df = hdr.fs / nsampl

freq_vec1 = -fs/2 : df : fs/2 - df

fft_y1 = EngeeDSP.Functions.fft(y1)     # FFT的

fft_Y1 = EngeeDSP.Functions.fftshift(fft_y1)    # 中心频谱
Out[0]:
7680-element Vector{ComplexF64}:
   0.4299999999999784 + 0.0im
    1.420074436857042 - 1.259402987504321im
  0.43216868301268896 + 2.18561269828764im
   1.4631953485700109 - 1.264690852185021im
  -0.6108928341574078 + 2.0868506337778356im
   -1.718166119515935 - 0.5277964184678581im
   -1.268020722866865 + 1.0992614404769583im
 -0.11131950936055546 + 0.48566398997173366im
  -1.3889234887099207 + 0.9665944738563965im
   0.8718361075312409 + 3.2228018323602328im
 -0.03253517186771049 + 0.3818156520900047im
  0.31213973318596056 + 1.7696508718261263im
  -0.9692764209245297 - 1.6720509877093566im
                      ⋮
  -0.9692764209245297 + 1.6720509877093566im
  0.31213973318596056 - 1.7696508718261263im
 -0.03253517186771049 - 0.3818156520900047im
   0.8718361075312409 - 3.2228018323602328im
  -1.3889234887099207 - 0.9665944738563965im
 -0.11131950936055546 - 0.48566398997173366im
   -1.268020722866865 - 1.0992614404769583im
   -1.718166119515935 + 0.5277964184678581im
  -0.6108928341574078 - 2.0868506337778356im
   1.4631953485700109 + 1.264690852185021im
  0.43216868301268896 - 2.18561269828764im
    1.420074436857042 + 1.259402987504321im

让我们计算分析信号的幅度谱，并将其绘制在0到10Hz的频率范围内。

In [ ]: 
Amp1 = (2*abs.(fft_Y1))/nsampl

 plot(
    freq_vec1, Amp1,
    xlabel="频率，赫兹",
    ylabel="的振幅",
    title="信号幅度谱△(hdr.录)NSRDB",
    grid=true,
    legend=false,
    xlim=(0, 10),
)
Out[0]:

具有正常窦性节律的ECG信号的所得幅度谱具有有序结构。 频谱清楚地显示了大约1Hz的区域中的主要频率分量，对应于心率，并且还观察到多个频率的谐波。 信号能量主要集中在高达10Hz的低频范围内。

另外，我们可以注意到大约0.1hz范围内的成分的存在，该成分与心脏活动没有直接关系，并且是由记录信号时发生的低频干扰引起的。

心电在室外的频谱分析

心室外壁是一种心律失常，其中异常收缩不是从窦房结发生，而是从心脏心室的额外激发灶发生。 与正常心动周期相比，这种过早收缩在ECG上显示为单独的，较早的，变形的QRS复合物。

胞外可以是单一的或重复的，并且取决于发生的频率和性质，它们被识别为频繁的或罕见的。 心室外壁通常不伴有患者明显的症状，但可能感觉中断或"错过"心跳。

我们使用来自[CU室性速心律失常数据库（CU VTADB）]的心室外心电图记录（https://physionet.org/content/cudb/1.0.0 /)。

In [ ]: 
data = ("$(@__DIR__)/cu-ventricular-tachyarrhythmia/cu05")
Out[0]:
"/user/DemoPublic/biomedical/spectral_analysis_arrhythmia/cu-ventricular-tachyarrhythmia/cu05"

我们将显示分析的心电图记录的主要参数。

In [ ]: 
hdr = read_header(data)

println("记录名称:   ", hdr.record)
println("采样率: ", hdr.fs, " 赫兹")
println("信号数量:   ", hdr.nsig)
println("计数数目:   ", hdr.nsamp)
println("数据格式:         ", hdr.fmt)
println("增益:              ", hdr.gain)
println("资料描述:              ", hdr.desc)

Наименование записи:   cu05
Частота дискретизации: 250.0 Гц
Количество сигналов:   1
Количество отсчётов:   127232
Формат данных:         212
Усиление:              400.0
Описание:              ECG

让我们来演示一下室外。 为此，我们将在126-128秒范围内的记录上选择一个两秒的片段，并在图形上用有色区域标记与晶外复合物相对应的区域。

In [ ]: 
t_start = 126.0
t_dur = 2.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y2 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t2 = (0:length(y2)-1) ./ hdr.fs .+ t_start

p =plot(t2, y2,
        xlabel="时间，从", 
        ylabel="阿姆利图达", 
        title="录录玫搂鲁拢潞(hdr.录)CU VTADB ",
        legend=false
        )

vline!(p, [127.1, 127.3], color=:red, alpha=0.3)

display(p)

我们将通过选择没有明显伪影的记录部分来演示心室外心电图的一个分钟片段。

In [ ]: 
t_start = 115.0
t_dur = 60.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y2 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t2 = (0:length(y2)-1) ./ hdr.fs .+ t_start

plot(t2, y2,
        xlabel="时间，从", 
        ylabel="的振幅", 
        title="录录玫搂鲁拢潞(hdr.录)CU VTADB ",
        legend=false)
Out[0]:

同样，我们使用EngeeDSP库的功能进行频率分析。

In [ ]: 
nsampl = length(y2)
fs = hdr.fs
df = hdr.fs / nsampl

freq_vec2 = -fs/2 : df : fs/2 - df

fft_y2 = EngeeDSP.Functions.fft(y2)     # FFT的

fft_Y2 = EngeeDSP.Functions.fftshift(fft_y2)    # 中心频谱
Out[0]:
15000-element Vector{ComplexF64}:
  0.24250000000006366 + 0.0im
   -0.617390662779087 - 0.9645818287347332im
  -0.5170981782064956 + 0.47210546222461947im
   0.3042996579577628 - 0.5017138577633im
  -0.6114165019779407 + 0.17445736541545964im
    0.518487255266888 + 0.8184288482754183im
  0.47517928893948813 + 0.5832875878615109im
   0.2709135407918337 - 0.0490081731987142im
   0.7785997524270325 + 0.8367519488401889im
   -0.621085745416029 - 0.12505334037132454im
 0.006377498767134782 + 0.0835351710709844im
   -0.814160502367459 + 0.584528873438301im
   -0.408193073323317 + 0.24900364154655996im
                      ⋮
   -0.408193073323317 - 0.24900364154655996im
   -0.814160502367459 - 0.584528873438301im
 0.006377498767134782 - 0.0835351710709844im
   -0.621085745416029 + 0.12505334037132454im
   0.7785997524270325 - 0.8367519488401889im
   0.2709135407918337 + 0.0490081731987142im
  0.47517928893948813 - 0.5832875878615109im
    0.518487255266888 - 0.8184288482754183im
  -0.6114165019779407 - 0.17445736541545964im
   0.3042996579577628 + 0.5017138577633im
  -0.5170981782064956 - 0.47210546222461947im
   -0.617390662779087 + 0.9645818287347332im

让我们绘制一个频率范围高达10赫兹的心室外心电图片段的振幅谱。

In [ ]: 
Amp2 = (2*abs.(fft_Y2))/nsampl

 plot(
    freq_vec2, Amp2,
    xlabel="频率，赫兹",
    ylabel="的振幅",
    title="信号幅度谱△(hdr.录)CU VTADB",
    grid=true,
    legend=false,
    xlim=(0, 10),
)
Out[0]:

与正常窦性心律相比，室外心电图的振幅谱具有较不有序的结构。 光谱能量分布的性质表明由于存在囊外复合物而违反心脏收缩的频率，这与心室囊外ECG的特征一致。

心房外心电图的光谱分析

我们使用[MIT-BIH心律失常数据库]的心电图记录（https://physionet.org/content/mitdb/1.0.0 /)，含有心律失常。 记录101通常特征在于主要是正常的窦性节律，针对该节律周期性地存在单个和重复的心房外晶体。

单次和重复的心房外壁是在正常窦性节律之外的心房中发生的过早心脏收缩。 对于单胞外，这种收缩偶尔出现，不会显着影响心脏的整体节律。 重复性心房外晶体的特征在于更频繁地发生过早冲动，这导致周期性地破坏心动周期的规律性。

In [ ]: 
data = ("$(@__DIR__)/mit-bih-arrhythmia/101")
Out[0]:
"/user/DemoPublic/biomedical/spectral_analysis_arrhythmia/mit-bih-arrhythmia/101"

我们将显示有关所选心电图记录的基本信息，包括采样参数，信号结构和记录的描述。

In [ ]: 
hdr = read_header(data)

println("记录名称:   ", hdr.record)
println("采样率: ", hdr.fs, " 赫兹")
println("信号数量:   ", hdr.nsig)
println("计数数目:   ", hdr.nsamp)
println("数据格式:         ", hdr.fmt)
println("增益:              ", hdr.gain)
println("资料描述:              ", hdr.desc)

Наименование записи:   101
Частота дискретизации: 360.0 Гц
Количество сигналов:   2
Количество отсчётов:   650000
Формат данных:         212
Усиление:              200.0
Описание:              MLII

让我们演示一个单一的心房外。 为此，我们将在记录上的间隔373-379秒中选择一个六秒片段，并用彩色垂直线标记图形上与晶外复合物相对应的区域。

In [ ]: 
t_start = 373.0
t_dur = 6.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end   = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y3 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t3  = (0:length(y3)-1) ./ hdr.fs .+ t_start

plot(t3, y3,
        xlabel="时间，从", 
        ylabel="阿姆利图达", 
        title="录录玫搂鲁拢潞(hdr.记录）MIT-BIH心律失常数据库 ",
        legend=false)
vline!([375.6, 376.2], color=:red, alpha=0.3)
Out[0]:

我们将演示包含心律失常发作的心电图记录的一分钟片段。

In [ ]: 
t_start = 370.0
t_dur = 60.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end   = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_212(data; sampfrom=sampl_start, sampto=sampl_end)
y3 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t3  = (0:length(y3)-1) ./ hdr.fs .+ t_start

plot(t3, y3,
        xlabel="时间，从", 
        ylabel="阿姆利图达", 
        title="录录玫搂鲁拢潞(hdr.记录）MIT-BIH心律失常数据库 ",
        legend=false)
Out[0]:

与前面的步骤类似，对于心电图的微小片段，我们将使用EngeeDSP库的功能进行频率分析，获得信号的中心频谱表示。

In [ ]: 
nsampl = length(y3)
fs = hdr.fs
df = hdr.fs / nsampl

freq_vec3 = -fs/2 : df : fs/2 - df

fft_y3 = EngeeDSP.Functions.fft(y3)     # FFT的

fft_Y3 = EngeeDSP.Functions.fftshift(fft_y3)    # 中心频谱
Out[0]:
21600-element Vector{ComplexF64}:
  -0.3499999999994543 + 0.0im
   -2.234148576509483 - 1.2921799198174142im
  -1.6923286675628049 - 3.0776448272429917im
   0.4223899856295361 + 4.918528645252934im
  0.06853869626723963 + 1.1634092443288893im
  0.07893171630887252 + 1.2391034578455589im
  0.23785772024982066 + 0.35861164346277974im
  0.13804578436253223 - 0.4131028992901946im
   0.9968657203876319 + 0.26269403356537424im
 -0.22749912173740938 + 0.05458454978024285im
 -0.10207517835085866 + 0.39308378482185447im
  -0.7325178393852525 + 0.02981812698969577im
 -0.11654664523090119 + 0.07441184108022014im
                      ⋮
 -0.11654664523090119 - 0.07441184108022014im
  -0.7325178393852525 - 0.02981812698969577im
 -0.10207517835085866 - 0.39308378482185447im
 -0.22749912173740938 - 0.05458454978024285im
   0.9968657203876319 - 0.26269403356537424im
  0.13804578436253223 + 0.4131028992901946im
  0.23785772024982066 - 0.35861164346277974im
  0.07893171630887252 - 1.2391034578455589im
  0.06853869626723963 - 1.1634092443288893im
   0.4223899856295361 - 4.918528645252934im
  -1.6923286675628049 + 3.0776448272429917im
   -2.234148576509483 + 1.2921799198174142im

让我们绘制频率范围高达10Hz的微小ECG片段的振幅谱。

In [ ]: 
Amp3 = (2*abs.(fft_Y3))/nsampl

 plot(
    freq_vec3, Amp3,
    xlabel="频率，赫兹",
    ylabel="的振幅",
    title="信号幅度谱△(hdr.记录）MIT-BIH心律失常数据库",
    grid=true,
    legend=false,
    xlim=(0, 10),
    ylim=(0 ,0.08)
)
Out[0]:

从MIT-BIH心律失常数据库中获得的分钟ECG片段的幅度谱与正常窦性节律相比，其特征在于明显的紊乱。 主要的低频分量存在于频谱中，但频谱峰值具有明显的展宽，并且信号能量分布在更宽的频率范围内。

心房颤动心电图的光谱分析

心房颤动是心律失常的一种形式，其中心房的电活动变得混乱和不协调。 而不是常规的窦房结冲动，多个唤醒灶出现在心房，这导致他们无效的收缩。

在心电图上，心房颤动的特征在于没有正常的P波和出现小的不规则基线波动，以及QRS复合物之间的不规则间隔。

我们从[心房颤动终止挑战数据库（AFDB）]数据库中选择与心房颤动对应的s02条目。]（https://physionet.org/content/afdb/1.0.0/)

In [ ]: 
data = ("$(@__DIR__)/af-termination-challenge/s02")
Out[0]:
"/user/DemoPublic/biomedical/spectral_analysis_arrhythmia/af-termination-challenge/s02"
In [ ]: 
hdr = read_header(data)

println("记录名称:   ", hdr.record)
println("采样率: ", hdr.fs, " 赫兹")
println("信号数量:   ", hdr.nsig)
println("计数数目:   ", hdr.nsamp)
println("数据格式:         ", hdr.fmt)
println("增益:              ", hdr.gain)
println("资料描述:              ", hdr.desc)

Наименование записи:   s02
Частота дискретизации: 128.0 Гц
Количество сигналов:   2
Количество отсчётов:   7680
Формат данных:         16
Усиление:              182.149
Описание:              ECG

让我们演示心电图记录的10秒片段，该片段显示心房颤动的不规则节律特征。

In [ ]: 
t_start = 0.0
t_dur   = 10.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end   = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_16(data; sampfrom=sampl_start, sampto=sampl_end)
y4 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t4  = (0:length(y4)-1) ./ hdr.fs .+ t_start

plot(t4, y4,
        xlabel="时间，从", 
        ylabel="阿姆利图达", 
        title="录录玫搂鲁拢潞(hdr.录)AFDB",
        legend=false)
Out[0]:

我们将展示一分钟的心房颤动心电图记录片段，反映明显的不规则心律。

In [ ]: 
t_start = 0.0
t_dur   = 60.0

sampl_start = Int(round(t_start * hdr.fs))
sampl_end   = sampl_start + Int(round(t_dur * hdr.fs))

_, raw, phys = read_dat_16(data; sampfrom=sampl_start, sampto=sampl_end)
y4 = phys === nothing ? float.(raw[:,1]) : phys[:,1]
t4  = (0:length(y4)-1) ./ hdr.fs .+ t_start

plot(t4, y4,
        xlabel="时间，从", 
        ylabel="的振幅", 
        title="录录玫搂鲁拢潞(hdr.录)AFDB",
        legend=false)
Out[0]:

与前面的步骤类似，我们将使用EngeeDSP库的功能对心房颤动心电图的微小片段进行频率分析，并获得中心信号谱。

In [ ]: 
nsampl = length(y4)
fs = hdr.fs
df = hdr.fs / nsampl

freq_vec4 = -fs/2 : df : fs/2 - df

fft_y4 = EngeeDSP.Functions.fft(y4)     # FFT的

fft_Y4 = EngeeDSP.Functions.fftshift(fft_y4)    # 中心频谱
Out[0]:
7680-element Vector{ComplexF64}:
 0.005490010925129241 + 0.0im
  -0.7423632719658406 + 1.9789740826734277im
  0.33492665860237025 - 1.5643893204067822im
  -0.7660751104496875 - 2.9121794152264586im
  -0.5678084664411713 + 0.9427265362270969im
  -0.5230071267951804 - 1.5753561875288806im
    1.182602214953216 - 0.7141930881550123im
  -0.7904291722913852 - 1.9998184339818996im
    1.542071912231664 - 0.48908491965557044im
  -1.9879774432424888 - 0.1390182820275765im
   1.0971526826344933 + 0.7966247851480954im
   -0.723195108699282 - 0.7038636138748355im
  -1.4293774193780902 - 0.9104137771854699im
                      ⋮
  -1.4293774193780902 + 0.9104137771854699im
   -0.723195108699282 + 0.7038636138748355im
   1.0971526826344933 - 0.7966247851480954im
  -1.9879774432424888 + 0.1390182820275765im
    1.542071912231664 + 0.48908491965557044im
  -0.7904291722913852 + 1.9998184339818996im
    1.182602214953216 + 0.7141930881550123im
  -0.5230071267951804 + 1.5753561875288806im
  -0.5678084664411713 - 0.9427265362270969im
  -0.7660751104496875 + 2.9121794152264586im
  0.33492665860237025 + 1.5643893204067822im
  -0.7423632719658406 - 1.9789740826734277im

让我们绘制频率范围高达10Hz的心房颤动中一分钟ECG片段的幅度谱。

In [ ]: 
Amp4 = (2*abs.(fft_Y4))/nsampl

 plot(
    freq_vec4, Amp4,
    xlabel="频率，赫兹",
    ylabel="的振幅",
    title="信号幅度谱△(hdr.录)AFDB",
    grid=true,
    legend=false,
    xlim=(0, 10),
)
Out[0]:

心房颤动中ECG的振幅谱的特征在于明显的紊乱和没有明确定义的基频分量。 信号能量分布在很宽的频率范围内，具有增加的频谱背景电平和大量不规则的频谱分量。

使用的材料

这项工作使用了来自开源PhysioNet的真实心电图记录。

  1. MIT-BIH正常窦性节律数据库（NSRDB）-具有正常窦性节律的心电图记录-<https://physionet.org/content/nsrdb/1.0.0 />
  2. 室性心律失常数据库(CUDB)-室性心律失常心电图记录-<https://physionet.org/content/cudb/1.0.0 />
  3. MIT-BIH心律失常数据库-各种心律失常的心电图记录-<https://physionet.org/content/mitdb/1.0.0 />
  4. 心房颤动数据库（AFDB）-心房颤动心电图记录-<https://physionet.org/content/afdb/1.0.0 />

结论

在示例过程中，提出了各种心律失常的光谱分析。 这项工作使用了从开放医学数据库获得的真实心电图信号。 EngeeDSP库的功能用于信号的频率分析，特别是功能 EngeeDSP.Functions.fft 来计算快速傅立叶变换和 EngeeDSP.Functions.fftshift 以形成居中的频谱表示。 这些仪器的使用使得可以直观地展示各种类型的心律失常的ECG谱的特征。