提示：本贴只做个人学习记录用
文章目录
前言一、爆发抑制比是什么？二、计算公式及代码1.论文研究方式 & 计算公式：2.遗忘因子 β3.最优阈值 θ & 分类误差 e4.matlab脚本
总结
前言
脑电信号的爆发抑制比(Burst Suppression Ratio, BSR)是麻醉深度监测中的重要指标，用于量化脑电图中爆发活动与抑制活动的比例。
本文内容主要来源于对论文《Real-time segmentation of burst suppression patterns in critical care EEG monitoring》的阅读和理解，相当于本论文的阅读笔记。
提示：以下是本篇文章正文内容，下面案例可供参考
一、爆发抑制比是什么？
爆发抑制特征：
爆发(Burst)：高振幅、相对高频的脑电活动
抑制(Suppression)：低振幅、近乎平坦的脑电活动段
BSR = 抑制时间 / 总时间 × 100%
爆发抑制模式在新生儿发育过程中以及各种新生儿脑损伤以及所有年龄段的弥漫性缺氧性脑损伤患者中都会自发出现。爆发抑制模式也可能由麻醉药物的使用诱导，例如在手术麻醉期间或用于治疗难治性癫痫持续状态，或者通过控制性低温疗法，这在心脏手术中的完全循环停止手术中很常见。
二、计算公式及代码
1.论文研究方式 & 计算公式：
实际采集临床脑电信号，并医生手动标记出抑制段和爆发段，最后用代码处理后的结果和医生手动标记的结果Y1 Y2做对比。
文章中使用的是递归方差分析（阈值+局部信号方差），具体的公式如下图：
上述公式中，μ是当前信号段的均值，σ是当前信号段的方差，x是信号数据，其中t、t-1表示不同的时间点的数据点索引，当前的计算结果与上一刻的数据有关。公式（1）和（2）就表示了自适应均值与自适应方差的计算过程。
其中δ表示某函数关系，即方差小于阈值θ的时候，就幅值为1，表示抑制段；否则（≥的时候），就赋值为0，表示爆发段。
上述公式中最关键的两个参数，遗忘因子 β 和分类阈值 θ。下面单独介绍两个参数的计算思路和流程。
2.遗忘因子 β
遗忘因子有两个特点：
（1）介于0~1之间；
（2）与遗忘时间 T 成反比：
其中Fs是信号采样率。
根据论文思路，进行了如下步骤操作：
（1）临床依据表明合适的遗忘时间在 0.01s~3s 之间；
（2）在（1）中的范围内进行逐点计算（精细化网格搜索），计算出各个 T 值对应的 β值，并参与 1 中的递归方差分析方程；
（3）求出各个 β值对应的最优阈值 θ；
（4）看最终结果，分类误差最小的 T（104.7ms），就是最优遗忘因子β（0.9534）。
3.最优阈值 θ & 分类误差 e
根据论文思路，进行了如下步骤操作：
（1）求出临床实际数据的每段的方差，按照方差降序排列；
（2）每个方差都视作一个阈值，进行 δ 函数关系的分类；
（3）每个 θ值都对应计算如下参数：
a. 识别为 Burst 的点数 Tb；
b. 识别为 suppression 的点数 Ts；
c. 实际的 Burst 的点数 Bb；
d. 实际的 suppression 的点数 Bs。
（4）当前遗忘因子 β值 和阈值 θ值对应的分类误差为：
找到误差 e值最小的阈值 θ （一般临床常用固定的5uV做阈值）和遗忘因子 β，即是最优结果。
论文中还有更进一步的优化结果，但是本贴中不进行细说，性能优化不大，下面上matlab脚本代码实例。
4.matlab脚本
%% EEG爆发-抑制模式检测脚本 (Matlab版本)
% 本脚本演示了如何在ICU脑电图(EEG)通道上检测爆发-抑制模式(burst-suppression patterns)
% 使用了一种基于局部方差阈值的简单方法
% 该方法在以下论文中被描述：
% 'Real-time segmentation of burst suppression patterns in critical care EEG monitoring',
% by Westover et al.
clear;
% close all;
fclose('all');
clc;
%% 参数设置
% 使用模拟数据进行演示（这样不需要实际的EDF文件）
use_simulation = 0;
if use_simulation
% 生成模拟的EEG信号（包含爆发-抑制模式）
fprintf('使用模拟数据进行演示...
');
samplingrate = 250; % 采样率（Hz）
duration = 60; % 模拟信号时长（秒）
% 生成模拟信号
sigbuf = generate_simulation_eeg(samplingrate, duration);
% 要绘制的时间间隔（单位：秒）
plot_interval = [5, 25]; % 绘制5-25秒的信号
else
% 读取真实脑电信号
load('...Annotations_1.mat');
load('...p1_data.mat')
samplingrate = Fs; % 采样率（Hz）
sigbuf = data';
plot_interval = [5, 25]; % 绘制5-25秒的信号
end
% 抑制阈值参数theta
theta = 50; % 按照5uV设置即可，考虑ADC增益，视情况而定
min_suppr_len = samplingrate/2; % 最小抑制段长度-0.5s
%% 使用递归方差估计对信号进行分段处理
[res, sig_segmented, sig_recvar] = recursive_variance_estimation(sigbuf(:,1), 0.9633, theta, 50);
sigbuf2 = (sigbuf - mean(sigbuf)) / (5 * rms(sigbuf)); % 信号归一化，画图方便
%%% Y1 Y2都是0是爆发段，1是抑制段，要处理一下
Y1AFT = ~Y1;
Y2AFT = ~Y2;
figure,
plot(sig_segmented,'r'); hold on
plot(sigbuf2(:,1),'b');hold on
plot(Y1AFT,'y');
plot(Y2AFT,'c');
%% 函数定义部分
% 递归方差估计函数，用于检测爆发-抑制模式
function [res, res_smoothed, res_var] = recursive_variance_estimation(sig, beta, theta, min_suppr_len)
% 输入：
% sig - 输入信号
% beta - 遗忘因子（默认：0.9633）
% theta - 抑制阈值（默认：50）
% min_suppr_len - 最小抑制段长度（默认：50）
% 输出：
% res - 原始分段结果（1=爆发，0=抑制）
% res_smoothed - 平滑后的分段结果（按照最小抑制段长度判断后，保留下来的抑制结果）
% res_var - 递归估计的方差值
% 参数默认值设置
if nargin < 2 || isempty(beta)
beta = 0.9633;
end
if nargin < 3 || isempty(theta)
theta = 50;
end
if nargin < 4 || isempty(min_suppr_len)
min_suppr_len = 50;
end
% 确保sig是列向量
sig = sig(:);
n = length(sig);
% 初始化结果数组
res = zeros(n, 1); % 原始分段结果（1=爆发，0=抑制）
res_var = zeros(n, 1); % 递归估计的方差值
res_smoothed = zeros(n, 1); % 平滑后的分段结果
% 初始均值和方差设为0
mu = 0; % 均值估计
sigsquare = 0; % 方差估计
% 对每个信号点进行递归处理
for i = 1:n
% 更新均值估计（指数加权移动平均）
mu = beta * mu + (1 - beta) * sig(i);
% 更新方差估计
sigsquare = beta * sigsquare + (1 - beta) * (sig(i) - mu)^2;
% 保存方差值
res_var(i) = sigsquare;
% 如果方差小于阈值theta，则标记为抑制状态(0)，否则为爆发状态(1)
res(i) = ~(sigsquare < theta);
end
%%%%%%% 上述这段就是论文中的公式：
% u(t) = β*u(t-1) + (1-β)*x(t)
% σ2(t) = β*σ2(t-1) + (1-β)*(x(t) - u(t))2
% z(t) = δ(σ2(t)<θ)，＜θ是1，抑制段；≥θ是0，爆发段
%%%%%%%%%%%%%%%%%%%%% 20251103 zxd 这的逻辑其实应该改一下（自己的想法），
%%%%%%%%%%%%%%%%%%%%% 实际使用的时候是以0.5s为单位判断一次，但下面是抑制段如果0.5s结束后，会一直判断到最后一个抑制点才转为爆发段
% 与原始方法相比的额外处理：
% 平滑处理：如果一个'抑制'期短于min_suppr_len样本点，则将其删除（视为伪抑制）
% is_suppressed_keep = false; % 标记当前抑制段是否需要保留
%
% for i = 1:n
% if res(i) == 0 % 当前点被标记为抑制
% if is_suppressed_keep
% % 已经确定这个点所在的抑制段足够长，直接保留它，不进行后续判断
% res_smoothed(i) = 0;
% else
% % 检查后面是否有足够多的抑制样本点
% % 如果从当前点开始的min_suppr_len个点都是抑制状态，则认为是有效的抑制
% if i + min_suppr_len - 1 <= n && max(res(i:i+min_suppr_len-1)) == 0 % (i + min_suppr_len - 1)是防止超出数据索引界限
% res_smoothed(i) = 0;
% is_suppressed_keep = true; % 标记为有效抑制段
% else
% % 抑制段太短，可能是伪影，视为非抑制状态
% res_smoothed(i) = 1;
% end
% end
% else % 当前点被标记为非抑制（爆发）
% % 重置标记
% is_suppressed_keep = false;
% res_smoothed(i) = 1;
% end
% end
%%%%%%%%%%%%%%%%%%%%% 抑制从当前点开始，长度维持超过min_suppr_len，则这一段都标记为抑制段，
%%%%%%%%%%%%%%%%%%%%% 否则就是爆发段
is_suppressed_keep = false; % 标记当前抑制段是否需要保留
nn = 0; % 计数
for i = 1:n
if res(i) == 0 % 当前点被标记为抑制
if is_suppressed_keep
% 已经确定这个点所在的抑制段足够长，直接保留它，不进行后续判断
res_smoothed(i) = 0;
nn = nn + 1;
if nn>=min_suppr_len % 当前段判断超过50个点之后，重置重新计数判断
is_suppressed_keep = false;
nn = 0;
end
else
% 检查后面是否有足够多的抑制样本点
% 如果从当前点开始的min_suppr_len个点都是抑制状态，则认为是有效的抑制
if i + min_suppr_len - 1 <= n && max(res(i:i+min_suppr_len-1)) == 0 % (i + min_suppr_len - 1)是防止超出数据索引界限
res_smoothed(i) = 0;
is_suppressed_keep = true; % 标记为有效抑制段
nn = nn + 1;
else
% 抑制段太短，可能是伪影，视为非抑制状态
res_smoothed(i) = 1;
nn = 0;
end
end
else % 当前点被标记为非抑制（爆发）
% 重置标记
is_suppressed_keep = false;
res_smoothed(i) = 1;
end
end
% 确保输出为双精度
res = double(res);
res_smoothed = double(res_smoothed);
end
%% 生成模拟的EEG信号（包含爆发-抑制模式）
function sig = generate_simulation_eeg(samplingrate, duration)
% 输入：
% samplingrate - 采样率（Hz）
% duration - 信号时长（秒）
% 输出：
% sig - 生成的模拟EEG信号
% 计算总样本数
n_samples = round(samplingrate * duration);
% 初始化信号为列向量
sig = zeros(n_samples, 1);
% 设置模拟参数
burst_duration = 0.5; % 爆发期持续时间（秒）
suppression_duration = 1.0; % 抑制期持续时间（秒）
% 生成爆发-抑制模式
i = 1;
while i <= n_samples
% 随机选择模式持续时间
current_burst_duration = burst_duration + (rand() - 0.5) * 0.3; % 添加一些随机性
current_suppression_duration = suppression_duration + (rand() - 0.5) * 0.4;
% 计算爆发期和抑制期的样本数
n_burst = round(current_burst_duration * samplingrate);
n_suppression = round(current_suppression_duration * samplingrate);
% 生成爆发期信号（包含多种频率成分）
burst_end = min(i + n_burst - 1, n_samples);
if burst_end >= i
% 混合多种频率的正弦波来模拟EEG爆发活动
freq1 = 5 + rand() * 10; % 5-15 Hz (alpha和beta活动)
freq2 = 0.5 + rand() * 4.5; % 0.5-5 Hz (delta和theta活动)
% 生成噪声和信号
noise_amp = 2; % 噪声振幅
signal_amp = 10 + rand() * 10; % 信号振幅
% 计算实际的信号索引
actual_idx = i:burst_end;
segment_length = length(actual_idx);
% 为当前段创建局部时间向量 - 确保是列向量
t_segment = (0:segment_length-1)';
t_segment = t_segment / samplingrate;
% 生成信号 - 确保是列向量
burst_signal = signal_amp * (sin(2*pi*freq1*t_segment) + 0.5*sin(2*pi*freq2*t_segment)) + noise_amp * randn(segment_length, 1);
% 验证维度是否匹配
if length(burst_signal) ~= length(actual_idx)
error('信号维度不匹配: burst_signal长度=%d, actual_idx长度=%d', length(burst_signal), length(actual_idx));
end
% 赋值到信号中
sig(actual_idx) = burst_signal;
end
% 生成抑制期信号（主要是噪声）
suppression_end = min(burst_end + 1 + n_suppression - 1, n_samples);
if suppression_end >= burst_end + 1
actual_idx = (burst_end + 1):suppression_end;
segment_length = length(actual_idx);
% 直接生成列向量噪声
suppression_signal = 0.5 * randn(segment_length, 1); % 非常低的振幅噪声
sig(actual_idx) = suppression_signal;
end
% 更新索引
i = suppression_end + 1;
end
% 添加一些人为的伪影（模拟运动伪影等）
n_artifacts = randi(3); % 随机添加1-3个人为伪影
for k = 1:n_artifacts
artifact_pos = randi(n_samples);
artifact_len = round(0.1 * samplingrate); % 100ms的伪影
artifact_start = max(1, artifact_pos - round(artifact_len/2));
artifact_end = min(n_samples, artifact_start + artifact_len - 1);
if artifact_end >= artifact_start
% 生成线性伪影 - 确保是列向量
artifact_amp = 50 + rand() * 50; % 50-100的振幅
artifact = artifact_amp * linspace(0, 1, artifact_end - artifact_start + 1)';
sig(artifact_start:artifact_end) = sig(artifact_start:artifact_end) + artifact;
end
end
% 标准化最终信号
max_abs_sig = max(abs(sig));
if max_abs_sig > 0
sig = sig / max_abs_sig * 100; % 缩放到合适的振幅范围
end
end
总结
以上是所有文章的阅读笔记梳理，欢迎大家使用。
© 版权声明
文章版权归作者所有，未经允许请勿转载。
相关文章
没有相关内容!
暂无评论...
