脑电图（EEG）数据的处理
该示例显示了使用专门的EngeeDSP库实现的脑电图（EEG）数字信号处理过程。 使用以标准EDF格式记录的真实数据示例，一致地说明了分析的关键阶段。
功能 engee.clear() 清理工作区:
engee.clear()
两种图形渲染模式可用于可视化结果。:
*对于静态，快速渲染，使用 gr().
*对于具有缩放和查看值功能的交互式动态可视化，请使用 plotlyjs().
通过注释掉第二个命令来选择其中一个命令。 两个后端是互斥的，最后一个被调用的后端将被激活。
# gr()
plotlyjs()
读取EDF文件
我们将使用函数进行连接 include 文件"edfread。jl"用于读取EDF文件:
include("$(@__DIR__)/edfread.jl")
功能 edfread 它旨在读取EDF（欧洲数据格式）格式的数据，这是存储生物医学信号的标准格式。 结构 hdr 此函数返回的信息包含有关记录的完整元信息:
一般记录参数:
-
ver-EDF格式版本 -
patientID-病人ID -
recordID-记录ID -
startdate和starttime-开始录制的日期和时间 -
bytes-以字节为单位的标头大小 -
records-文件中的数据块数 -
duration-以秒为单位的一个区块的持续时间 -
ns-记录中的通道数
每个通道的参数:
-
labels-频道名称 -
transducers-传感器类型 -
physicalDims-物理测量单位 -
physicalMins和physicalMaxs-最小和最大物理值 -
digitalMins和digitalMaxs-最小和最大数字值 -
prefilters-在录制过程中应用的过滤器 -
samples-每个通道的一个块中的样本数
脑电图（EEG）数据将作为生物医学信号处理的一个例子。 让我们从上传文件开始 S004R02.edf 来研究其结构并构建信号的波形。
hdr, record = edfread("$(@__DIR__)/S004R02.edf")
下一步是分析EDF文件的元数据。 我们将提取用于进一步处理的关键参数包括采样率，通道数和记录持续时间。
println("格式版本: ", hdr.ver)
println("病人ID: ", strip(hdr.patientID))
println("记录的描述: ", strip(hdr.recordID))
println("开始日期/时间: ", hdr.startdate, " ", hdr.starttime)
println("频道数目: ", hdr.ns)
println("参赛作品数目: ", hdr.records)
println("总工期: ", hdr.records * hdr.duration, " 证券交易委员会")
为了清楚起见，我们还介绍了前三个通道的关键特性。
for i in 1:min(3, size(record, 1))
channel_data = record[i, :]
fs = hdr.samples[i] / hdr.duration
t = (1:length(channel_data)) ./ fs
println(" $(i). $(hdr.labels[i]):")
println(" 频率是可自由支配的。: ", round(fs, digits=1), " 赫兹")
println(" 信号长度: ", length(channel_data), " 计数")
println(" 持续时间: ", round(t[end], digits=2), " 证券交易委员会")
println(" 最小/最大振幅: ", round(minimum(channel_data), digits=2), " .. ", round(maximum(channel_data), digits=2))
end
接下来，所选择的脑电通道的信号被可视化。 波形的构建允许您直观地评估大脑活动的性质和伪影的存在。
# 渠道选择
Nk1 = 1
Nk2 = 17
eegNk1 = record[Nk1, :]
eegNk2 = record[Nk2, :]
# 信号参数
fs = hdr.samples[1] / hdr.duration # 采样率
Tn = 60.0 # 观察时间
N1 = Int64(Tn * fs) # 样本长度
dims = hdr.physicalDims[Nk1]
# 信号提取
y1 = eegNk1[1:N1]
y2 = eegNk2[1:N1]
t = range(0, step=1/fs, length=N1)
p1 = plot(t, y1,
title="脑电图波形（通道№Nk1）",
xlabel="时间，从",
ylabel="振幅，△暗淡",
grid=true)
p2 = plot(t, y2,
title="脑电图波形($Nk2通道)",
xlabel="时间，从",
ylabel="振幅，△暗淡",
grid=true)
脑电信号处理
脑电信号的频谱分析
接下来，我们使用EngeeDSP库切换到频率数据分析。 对于选定的通道，应用快速傅立叶变换，将信号从时间表示转换为频率表示。 所获得的光谱被居中，以便随后对频率分量进行可视化和分析。
nsamp = length(y1)
df = fs / nsamp
freq_vec = -fs/2 : df : fs/2 - df
fft_y1 = EngeeDSP.Functions.fft(y1)
fft_y2 = EngeeDSP.Functions.fft(y2)
# 中心频谱
fft_Y1 = EngeeDSP.Functions.fftshift(fft_y1)
fft_Y2 = EngeeDSP.Functions.fftshift(fft_y2)
在接收到信号的频率表示后，其幅度谱被可视化。
Amp1 = abs.(fft_Y1)./ nsamp
p3 = plot(
freq_vec, Amp1,
xlabel="频率，赫兹",
ylabel="振幅，△暗淡",
title="振幅谱(通道△Nk1)",
grid=true,
legend=false
)
对于另一个信道，执行频谱到分贝的附加转换。
Amp2 = abs.(fft_Y2)./nsamp
# 幅度单位为dB（相对于1mv）
Amp2_dB = 20 .* log10.(Amp2 .+ eps()) # 添加eps以避免-Inf
p4 = plot(
freq_vec, Amp2_dB,
xlabel="频率，赫兹",
ylabel="幅度，dB(rel. 1元)",
title="振幅谱以dB为单位(通道∶Nk2)",
grid=true,
legend=false
)
我们使用函数计算第一个通道的周期图 periodogram 来估计频谱功率密度(SPM)。 曲线图示出了信号能量的频率分布。
power,freq = EngeePhased.Functions.periodogram(y1, fs=fs,out=:data)
p5 = plot(vec(freq), power, title="脑电图周期图($Nk1通道)", xlabel="频率，赫兹", ylabel="SPM MV^2/Hz", grid=true)
在光谱上观察到明显的网络干扰（60Hz）。 在分析脑电节律之前，要求通过过滤来消除它。
对脑电信号进行滤波
频谱分析结果显示，频率为60Hz时有一个明显的峰值，对应于来自交流电网的尖端。 为了使用EngeeDSP库消除它，使用了10阶Butterworth陷波滤波器。 滤波器配置为抑制57Hz至63Hz的频段，谐振频率为60Hz，从而消除网络干扰。
f_center = 60.0 # 谐振频率，Hz
bw = 6.0 # 带宽，Hz
f1 = f_center - bw/2 # 下限频率，Hz
f2 = f_center + bw/2 # 上限频率，Hz
wn1 = f1 / (fs/2)
wn2 = f2 / (fs/2)
order = 10 # 过滤顺序
b, a = EngeeDSP.Functions.butter(order, [wn1, wn2], "stop")
两个通道的信号都由陷波滤波器通过函数处理 filter EngeeDSP图书馆。 这使得能够抑制60Hz的网络干扰。
filt_y1 = EngeeDSP.Functions.filter(b, a, y1)
filt_y2 = EngeeDSP.Functions.filter(b, a, y2)
网络干扰抑制后，构建脑电波形。
p6 = plot(
t, filt_y1,
xlabel = "时间，c",
ylabel = "振幅，MV",
title = "滤波后的脑电图波形（通道△Nk1）",
grid = true,
legend = false
)
使用函数 EngeeDSP.Functions.fft 进行快速傅立叶变换，然后以频谱为中心。 构建的振幅谱允许您在60Hz频率下直观地评估不存在峰值。
fft_filt_y1 = EngeeDSP.Functions.fft(filt_y1)
fft_filt_y1_shift = EngeeDSP.Functions.fftshift(fft_filt_y1)
Amp_y1_filt = abs.(fft_filt_y1_shift)./nsamp
p7 = plot(
freq_vec, Amp_y1_filt,
xlabel = "频率，Hz)",
ylabel = "振幅，△暗淡",
title = "滤波后的脑电信号的频谱(通道△Nk1)",
grid = true,
legend = false
)
为了直观地评估陷波滤波器的操作，在干扰抑制之前和之后比较信号的幅度谱。 一张图表显示了原始频谱和滤波后的频谱。
p8 = plot(freq_vec, Amp1,
label = "没有过滤"
)
plot!(freq_vec, Amp_y1_filt,
label = "过滤后"
)
xlabel!("频率，赫兹")
ylabel!("振幅，△暗淡")
title!("振幅谱(△Nk1通道)")
脑节律的隔离
带通FIR滤波器（有限脉冲响应）用于从纯化的脑电信号中提取α节律（8-12Hz）。 使用函数 EngeeDSP.Functions.fir1 正在设计具有给定带宽的200阶滤波器。 然后函数 EngeeDSP.Functions.filter 将此滤波器应用于信号 filt_y1 通过隔离alpha节奏范围内的分量并形成新信号 alpha_y1.
f_low = 8.0 # 下限频率，Hz
f_high = 12.0 # 上限频率，Hz
wn = [f_low, f_high] ./ (fs/2)
fir_order = 200 # 过滤顺序
b_alpha = EngeeDSP.Functions.fir1(fir_order, wn)
alpha_y1 = EngeeDSP.Functions.filter(b_alpha, [1.0], filt_y1)
构建所选alpha节奏的示波图 - alpha_y1.
# 阿尔法节奏波形
p9 = plot(
t, alpha_y1,
xlabel = "时间，从",
ylabel = "振幅，MV",
title = "Alpha节奏（8-12Hz），（通道№Nk1）",
legend = false,
grid = true
)
display(p9)
为了进一步分析，设置了大脑活动的主要节律的特征频率范围。:
*三角洲节奏： 0.5 – 4 赫兹
*Theta节奏： 4 – 8 赫兹
*阿尔法节奏： 8 – 12 赫兹
*Beta节奏： 13 – 30 赫兹
这些范围将用于调整带通滤波器，以便将相应的韵律分量与EEG信号隔离。
delta_band = (0.5, 4.0)
theta_band = (4.0, 8.0)
alpha_band = (8.0, 12.0)
beta_band = (13.0, 30.0)
为了更容易区分不同的脑电图节奏，已经创建了一个函数 band_fir_filter，其中使用EngeeDSP库实现了典型的FIR带通滤波操作
function band_fir_filter(sig, fs, band, order)
f1, f2 = band
wn = [f1, f2] ./ (fs/2)
b = EngeeDSP.Functions.fir1(order, wn)
y = EngeeDSP.Functions.filter(b, [1.0], sig)
return y
end
使用先前创建的函数 band_fir_filter 从脑电信号中依次区分三个主要节奏分量：delta节奏（0.5-4Hz），theta节奏（4-8Hz）和beta节奏（13-30Hz）。
# 滤波第一信号
delta_y1 = band_fir_filter(filt_y1, fs, delta_band, fir_order)
theta_y1 = band_fir_filter(filt_y1, fs, theta_band, fir_order)
beta_y1 = band_fir_filter(filt_y1, fs, beta_band, fir_order)
同样，对于第二个信号，使用相同的带通滤波功能分离大脑的主要节律。 从信号 filt_y2 获得相应的分量：delta，theta，alpha和beta节奏。
# 滤波所述第二信号
delta_y2 = band_fir_filter(filt_y2, fs, delta_band, fir_order)
theta_y2 = band_fir_filter(filt_y2, fs, theta_band, fir_order)
alpha_y2 = band_fir_filter(filt_y2, fs, theta_band, fir_order)
beta_y1 = band_fir_filter(filt_y2, fs, beta_band, fir_order)
第一信号的脑节律的可视化。
p10 = plot(
plot(t, delta_y1,
xlabel = "时间，从",
ylabel = "Amlituda，mkV",
title = "三角洲节奏(0.5-4赫兹)",
legend = false,
grid = true),
plot(t, theta_y1,
xlabel = "时间，从",
ylabel = "Amlituda，mkV",
title = "Θ节奏(4-8赫兹)",
legend = false,
grid = true),
plot(t, alpha_y1,
xlabel = "时间，从",
ylabel = "Amlituda，mkV",
title = "阿尔法节奏(8-12赫兹)",
legend = false,
grid = true),
plot(t, beta_y1,
xlabel = "时间，从",
ylabel = "Amlituda，mkV",
title = "贝塔节奏(13-30赫兹)",
legend = false,
grid = true),
layout = (4, 1),
size = (900, 900),
margin = 40*Plots.px
)
display(p10)
在分析脑电图时，信号被划分为标准频率范围，每个频率范围都与特定的大脑状况相关联。:
-
Delta节奏（0.5-4Hz）-最慢的节奏。 在深度睡眠期间占主导地位。
-
Θ节律（4-8赫兹）-与嗜睡，放松和记忆过程的状态有关。
-
阿尔法节奏（8–12赫兹）-平静清醒的节奏在眼睛闭合时最明显。
-
贝塔节律（13-30赫兹）-与活跃的觉醒，注意力和精神活动有关。
相关性分析
让我们继续到相关数据处理部分。 相关性分析用于量化所选节奏的内部结构和频率。
自相关函数
使用函数 EngeeDSP.Functions.xcorr 信号的自相关来计算 alpha_y1. 函数返回两个数组:
-
R_aa-每个移位的相关系数的值。 -
lag_aa-计数中的时移的相应值。
maxlag = length(alpha_y1) - 1
R_aa, lag_aa = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y1, maxlag)
tau_aa = vec(lag_aa) ./ fs
Alpha节奏的自相关函数的可视化
R_aa_norm = R_aa ./ maximum(R_aa)
plot(
tau_aa, R_aa_norm,
xlabel = "转变，从",
ylabel = "ACF，mkV^2",
title = "阿尔法节奏ACF(通道№Nk1)",
grid = true,
legend = false
)
互相关函数
为了分析不同大脑区域中α活动的关系，计算了第一和第二通道的α节律信号之间的互相关函数（VCF）。
R_ab, lag_ab = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y2, maxlag)
tau_ab = vec(lag_ab) ./ fs
Α节律相互关联的可视化
R_ab_norm = R_ab ./ maximum(R_ab)
plot(
tau_ab, R_ab_norm,
xlabel = "转变，从",
ylabel = "VKF，mkV^2",
title = "Α节律的互相关函数",
grid = true,
legend = false
)
结论
在本例中，演示了使用EngeeDSP库以标准EDF格式记录的EEG数据的处理。 所有的工作都是基于它的功能，并以逻辑顺序执行。:
-
**数据下载：**数据是从EDF文件中读取的。
-
**光谱分析：**使用功能 **
EngeeDSP.Functions.fft**及 **EngeeDSP.Functions.fftshift**信号被转换到频域，这揭示了网络干扰的存在。 -
**滤波：**为了抑制60Hz干扰，设计并应用了Butterworth陷波滤波器，使用以下功能 **
EngeeDSP.Functions.butter**及EngeeDSP.Functions.filter. -
**突出节律：**大脑的主要节律（delta，theta，alpha，beta）使用该函数创建的带通FIR滤波器进行隔离 **
EngeeDSP.Functions.fir1**并通过EngeeDSP.Functions.filter. -
**相关性分析：**为了研究信号的时间结构和关系，使用以下方法计算了自相关（ACF）和互相关（VCF）函数
EngeeDSP.Functions.xcorr.