-1
提示：仅作个人学习用
前言
提示：仅作个人学习用：
麻醉深度BIS模块的matlab脚本代码复现，可能还有些bug和细节需要调整，各位如果拿去用的话参照自己的使用需求调整规划吧。
提示：以下是本篇文章正文内容，下面案例可供参考
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 用公共数据集完成BIS模块的matlab脚本调试 %
% 完整的模块代码 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;
fclose('all');
clc;
load('...\0006.mat');
% 删掉nan的数值，防止出现滤波导致的全部nan值的问题
% EEG1 = double(EEG1);
nan_mask = isnan(EEG1);
EEG1(nan_mask) = [];
Fs = 128; % 采样频率 128 Hz
total_points = length(EEG1); % 总点数 120000
% 计算窗口参数
window_duration = 2; % 窗长 2秒
window_length = Fs * window_duration; % 窗长点数: 256点
stride = 0.5; % 步长 0.5秒
nStride = Fs*stride; % 窗口移动一个步长所移动的点数
step_length = Fs * nStride; % 步长点数: 64点
N=floor((length(EEG1)-Fs)/nStride)-10; % 当前数据长度，窗口总共可移动多少次
overlap_length = window_length - step_length; % 重叠点数: 192点
%% 计算BSR & BC
%%%%%%%%%%%%%%%%%%%%%%%%%% 计算抑制比——计算窗长1s，步长0.5s
BSRmap = zeros(N, 1);
for n = 1:N % 循环取每个步长
x = EEG1((n + 6.5)*nStride : (n + 6.5 + (1/stride))*nStride); % 取每个步长取1s的数据
v = (1:length(x))' .^(0:1);
baselineX = v * (v\x); % a=v\x是方程va=x的解。即v的逆矩阵*矩阵x
BSRmap(n) = all(abs(x-baselineX) <= 5); % ≤5，是抑制段记录为1
end
BSR = 100 * movmean(BSRmap, [(63/stride) - 1, 0]);
%%%%%%%%%%%%%%%%%%%%% 整理参考另一篇论文中的方式进行了复现，都可以用，看怎么取舍吧
% 抑制阈值参数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');
%% 计算BIS结果
%%%%%%%%%%%%%%%%%%%%%%%%%% 过高通滤波
% [B, A] = butter(2, 0.65/(Fs/2), 'high'); % 创建0.65 Hz的二阶高通巴特沃斯滤波器
A = [1,-1.95488437262939,0.955879873967924];
B = [0.977691061649329, -1.95538212329866,0.977691061649329];
eegHiPassFiltered = filter(B, A, EEG1); % 过高通滤波器
psd = nan(N, 4 * nStride/2); % 分配空间来存储功率谱密度，每移动一次窗口都有一组结果（7509个结果），每组结果是0到63.5 Hz， 0.5 Hz的频率分辨率（128长度）
%%%%%%%%%%%%%%%%%%%%%%%%%% 计算主要参数
xp = [0, 3, 6];
yp = [0, 0.25, 1];
x = 0:0.5:63.5;
xLimit = min(max(x, 0), 6); % 将输入信号x限制在 0~6 范围内
suppressionFilter = interp1(xp, yp, xLimit).^2; % 定义一个0 - 6hz范围内的抑制滤波器
components = nan(N, 3); % 建个空矩阵
for n = 1:N % 对可用于计算psd的epochs数量进行评估，每个epoch窗长2s，步长0.5s，所以每个窗口间有75%的数据是重叠的
% 每0.5s输出一组计算结果
if (n>=(window_duration/stride)) && (sum(BSRmap(n-((window_duration/stride)-1):n))==0) % 如果这段2s的数据，不包含任何抑制段
x = eegHiPassFiltered((n+(window_duration/stride))*nStride + (1:(window_duration/stride)*nStride));
v = (1:length(x))' .^(0:1);
baselineX = v * (v\x);
FFTRet = fft(blackman(length(x)) .* (x - baselineX)); % 执行布莱克曼窗口FFT，消除基线信号漂移（减去baseline）
psd(n,:) = 2 * abs(FFTRet(1:length(x)/2)').^2 / (length(x)*sum(blackman(length(x)).^2)); % 将频率幅值转换为频率功率
RawEEG = EEG1((n+(window_duration/stride))*nStride + (1:(window_duration/stride)*nStride));
if sawtoothDetector(RawEEG, nStride) % 如果检测到存在锯齿形k成分
psd(n, :) = suppressionFilter .* psd(n, :); % 相应地抑制低频成分的占比（防止k成分对低频结果造成影响，导致最终输出的BIS值偏低？？）
end
end
thirtySec = max(1, (n -(30/stride)+1)):n; % 最近30s的数据，窗口是0.5s长度，所以是找最近60个点的索引
vhighPowerConc = sqrt(mean(psd(thirtySec, (((39.5/0.5):(46.5/0.5)) +1)) .* psd(thirtySec, (((40/0.5):(47/0.5)) +1)), 2));
% 39.5-46.5Hz频率成分和40-47Hz频率成分，对应相乘后求平均，再做根号
wholePowerConc = sqrt(mean(psd(thirtySec, (((0.5/0.5):(46.5/0.5)) +1)) .* psd(thirtySec, (((1/0.5):(47/0.5)) +1)), 2)); % 截取0.5-46.5Hz的频段成分
% 0.5-46.5Hz频率成分和1-47Hz频率成分，对应相乘后求平均，再做根号
% 上述都是一次计算0.5s的值
MidPSD = psd(thirtySec, (((11/0.5):(20/0.5)) +1)); % 截取出来11-20Hz频率的成分
midBandPower = prctmean(nanmean(10*log10(MidPSD),1),50, 100); % 计算11-20Hz频率的均值
Psd30To47 = psd(thirtySec, (((30/0.5):(47/0.5)) +1)); % 取出30-47Hz部分的psd
components(n ,1) = mean(10*log10(v(~isnan(v)))) - midBandPower;
% 10*log10转化dB，之后计算均值（排除掉nan值）
components(n, 2) = trimmean (10 * log10(vhighPowerConc./ wholePowerConc), 50); % 全麻指数（40-47Hz/1-47Hz）
Psd05To4 = psd(thirtySec, (((0.5/0.5):(4/0.5)) +1)); % 取出0.5-4Hz部分的psd
components(n ,3) = mean(10*log10(v(~isnan(v)))) - midBandPower;
% 10*log10转化dB，之后计算均值（排除掉nan值）
end
%%%%%%%%%%%%%%%%%%%%%%%%%% 加权计算得出最终的麻醉深度值结果
%%% 将成分1映射到logistic s曲线上，获得镇静评分
sedationScore = 104.4 - 49.4./(1+exp((components(:, 1)-(-13.9))/5.29));
%%% 将成分2映射到线性曲线上，并映射到logistic s曲线，获得全麻指数
componentsBoud = min(max(components(:,2), -60.89), -30);
% 在这里把成分2限制在[-60.89, -30]的区间内，因为后需要插值，x不能超出去
generalScore1 = interp1([-60.89, -30], [-40, 43-1], componentsBoud);
% interp1()函数————插值用（默认使用线性插值），按照x=[-60.89, -30]，y=[-40,43-1]，这个函数曲线，对刚刚限制后的成分2进行插值
generalScore2 = 61.3 - 72.6./(1+exp((components(:,2)-(-24.0))/3.55)); % % 在logistic s曲线上将成分2映射
generalScore2 = generalScore2.*(components(:,2)>= -30);
generalScore = generalScore1 + generalScore2;
bsrScore = piecewise(BSR, [0, 100], [50, 0]); % 使用分段线性函数将BSR转换为BSR分数
generalWeight = piecewise( components(:,3), [0,5], [0.5, 1]) .* (generalScore < sedationScore); % 将分量3转换为权重
bsrWeight = piecewise(BSR, [10, 50], [0, 1]); % 将BSR转换为权重
x = (sedationScore .* (1-generalWeight)) + (generalScore .* generalWeight); % 将镇静和全麻评分一起加权
MyOpeniBis = piecewise(x, [-40, 10, 97, 110], [0, 10, 97, 100]).*(1-bsrWeight) + bsrScore.*bsrWeight; % 用BSR压缩和加权
%%%%%% 最终的结果是0.5s一个值
MyBIS.BIS = MyOpeniBis(1:2:end);
figure,
plot(BIS,'b');hold on
plot(MyBIS.BIS,'r');
legend('BIS','MyBIS');
%% 计算SEF & MF & PPF & TP & EMG & 各频段占比 & Alpha/Delta（边界简单验证过）
%%%%%%%%%%%%%%%%%%%%%%%%%%% 生成调试用信号 %%%%%%%%%%%%%%%%%%%%%%%%%%%
% 参数设置
amplitude = 100; % 幅值
Fs = 128; % 采样率 128Hz
f = 30; % 频率 5Hz
duration = 100; % 持续时间 1秒
% 生成时间向量
t = 0:1/Fs:duration-1/Fs;
% 生成信号
eegHiPassFiltered = amplitude * sin(2*pi*f*t)';
N=floor((length(eegHiPassFiltered)-Fs)/nStride)-10; % 当前数据长度，窗口总共可移动多少次
%%%%%%%%%%%%%%%%%%%%%%%%%%% 生成调试用信号 %%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 1:N % 对可用于计算psd的epochs数量进行评估，每个epoch窗长2s，步长0.5s，所以每个窗口间有75%的数据是重叠的
if (n>=(window_duration/stride)) && (sum(BSRmap(n-((window_duration/stride)-1):n))==0) % 如果这段2s的数据，不包含任何抑制段
x = eegHiPassFiltered((n+(window_duration/stride))*nStride + (1:(window_duration/stride)*nStride));
v = (1:length(x))' .^(0:1);
baselineX = v * (v\x);
FFTRet = fft(x - baselineX); % 执行布莱克曼窗口FFT，消除基线信号漂移（减去baseline）
% 计算幅度谱
magnitude = abs(FFTRet);
% 计算功率谱 (|FFT|^2)
power_spectrum = magnitude .^ 2;
% 归一化：除以(采样点数 * 采样频率)
PSD_full = power_spectrum / (N * Fs);
% 取单边谱（由于实信号的对称性）
num_unique_pts = length(x)/2 + 1;
PSD = PSD_full(1:num_unique_pts);
PSD(2:end-1) = 2 * PSD(2:end-1); % 乘以2补偿负频率部分
%% 计算TP
freq_range = 0.5/stride+1:30/stride+1; % matlab索引从1开始，所以都+1 (2:61)对应频率(0.5:30)Hz
TP(n) = sum(PSD(freq_range));
%% 计算SEF 95%
cumulative_power = cumsum(PSD(freq_range)); % 计算累积功率
cumulative_percentage = cumulative_power / TP(n); % 计算累积百分比
idx = find(cumulative_percentage >= 0.95, 1, 'first'); % 找到第一个超过95%的位置
% indices = find(condition, k, direction);
% condition：逻辑条件表达式; k：要找到的元素数量; direction：搜索方向（'first'或'last'）
SEF(n,1) = (idx)*0.5; % matlab索引从1开始，C索引从0开始，应该要+1，C化的时候注意这个细节，调试确认是否需要+1
%% 计算MF
idx = find(cumulative_percentage >= 0.5, 1, 'first'); % 找到第一个超过50%的位置
MF(n,1) = (idx)*0.5; % matlab索引从1开始，C索引从0开始，应该要+1，C化的时候注意这个细节，调试确认是否需要+1
%% 计算PPF
[max_val, max_idx] = max(PSD(freq_range)); % 找到最大值的位置
PPF(n,1) = (max_idx)*0.5; % matlab索引从1开始，C索引从0开始，应该要+1，C化的时候注意这个细节，调试确认是否需要+1
%% 计算EMG(70-100Hz)
% freq_range = 70/stride:100/stride;
% EMG(n,1) = sum(PSD(freq_range));
%% 计算各频段占比 & Alpha/delta
% Alpha waves: 8 to 13 Hz
% Beta waves: 13 to 30 Hz
% Theta waves: 4 to 8 Hz
% Delta waves: 0.5 to 4 Hz
Alpha(n,1) = 100*sum(PSD(8/stride+1:13/stride))/sum(PSD(freq_range)); % 这些索引都是左闭右开的区间，防止频点在边缘处重合
Beta(n,1) = 100*sum(PSD(13/stride+1:30/stride+1))/sum(PSD(freq_range)); % 只有beta是左闭右闭
Theta(n,1) = 100*sum(PSD(4/stride+1:8/stride))/sum(PSD(freq_range));
Delta(n,1) = 100*sum(PSD(0.5/stride+1:4/stride))/sum(PSD(freq_range));
BandSum(n,1) = Alpha(n)+Beta(n)+Theta(n)+Delta(n); % 各频段占比总和应该是100%，这里做个验证
AtoD(n,1) = 100*Alpha(n,1)/Delta(n,1);
end
end
MyBIS.TP = TP(1:2:end);
MyBIS.SEF = SEF(1:2:end);
MyBIS.MF = MF(1:2:end);
MyBIS.PPF = PPF(1:2:end);
MyBIS.EMG = EMG(1:2:end);
MyBIS.Alpha = Alpha(1:2:end);
MyBIS.Beta = Beta(1:2:end);
MyBIS.Delta = Delta(1:2:end);
MyBIS.Theta = Theta(1:2:end);
MyBIS.AtoD = AtoD(1:2:end);
%% 计算SQI
%% 函数定义
% 递归方差估计函数，用于检测爆发-抑制模式
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;
res(i) = ~(sigsquare < theta);
end
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
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;
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
总结
上述文章中的计算方式均来自公开论文和公开数据集。
信号质量SQI的计算没写，暂时空余，后续会专门研究后发一篇博文，感兴趣的话可以关注后续更新情况。
扫一扫
1244
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
