麻醉深度监护Openibis学习笔记

zxdnjp

已于 2025-10-28 09:48:22 修改

阅读量456

点赞数 14

CC 4.0 BY-SA版权

分类专栏：脑电信号处理学习笔记文章标签：学习笔记信号处理

于 2025-10-25 16:00:04 首次发布

本文链接：https://blog.youkuaiyun.com/zxdnjp/article/details/153870524

脑电信号处理学习笔记专栏收录该内容

5 篇文章

订阅专栏

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章目录

前言
一、文献引用和介绍
二、代码及注释（仅作个人学习记录）
总结

前言

提示：仅作个人学习记录用：

利用脑电信号监测麻醉深度模块的代码学习记录。代码中有详细注释对运算步骤进行说明

提示：以下是本篇文章正文内容，下面案例可供参考

一、文献引用和介绍

参考文献：Open Reimplementation of the BIS Algorithms for Depth of Anesthesia
作者：Christopher W Connor, MD PhD [Associate Professor of Anaesthesia, Attending Anesthesiologist, Adjunct Associate Professor]
上述论文是本篇代码学习笔记记录的基础，仅对文中的代码做解读和注释，不进行任何修改，以博主自己浅显易懂的方式做记录，请各位使用的时候也请引用原文献，并声明出处！！

二、代码及注释（仅作个人学习记录）

1.读取数据并定义信号相关基础参数

代码如下（示例）：

clear all;
close all;
clc;
fclose('all');

load('...\BIS数据库\data\case1.mat');
Fs = 128;       % 采样频率
stride = 0.5;   % 步长
nStride = Fs*stride;    % 窗口移动一个步长所移动的点数
N=floor((length(EEG)-Fs)/nStride)-10;   % 当前数据长度，窗口总共可移动多少次

2.计算爆发抑制比

BSRmap = zeros(N, 1);
for n = 1:N     % 循环取每个步长
    % x = EEG((n + 6.5)*nStride + (1:2*nStride));     % 取每个步长取1s的数据
    x = EEG((n + 6.5)*nStride : (n + 6.5 + 2)*nStride);     % 取每个步长取1s的数据
    v = (1:length(x))' .^(0:1);   
    baselineX = v * (v\x);      % a=v\x是方程va=x的解，关于矩阵左除的运算，可以deepseek看具体示例能够更好的明白
    BSRmap(n) = all(abs(x-baselineX) <= 5);     % ≤5，是抑制段记录为1
end
BSR = 100 * movmean(BSRmap, [(63/stride) - 1, 0]);

3.过高通滤波

%% 过高通滤波
[B, A] = butter(2, 0.65/(Fs/2), 'high');     % 创建0.65 Hz的二阶高通巴特沃斯滤波器
eegHiPassFiltered = filter(B, A, EEG);       % 过高通滤波器

4.计算主要参数

xp = [0, 3, 6];
yp = [0, 0.25, 1];
x = 0:0.5:63.5;
xLimit = min(max(x, 0), 6);     % 将输入信号x限制在 0~6 范围内
% if x<0, x=0; if x>6, x=6;
suppressionFilter = interp1(xp, yp, xLimit).^2;      % 定义一个0 - 6hz范围内的抑制滤波器
% interp1()函数是对已知的(xp, yp)数据点（xp yp画条曲线），计算在xLimit位置处的函数值（输出预估的y值，相当于进行了一个线性插值运算），xlimit超过xp的值会返回NaN；
% 上面已经把x限定在了0-6范围内，所以最高输出的y值就是1，前面没到6的时候都是进行插值的结果

%%% 写C语言的时候可以直接把suppressionFilter这个系数写死，节省计算资源

components = nan(N, 3);     % 建个空矩阵
for n = 1:N    % 对可用于计算psd的epochs数量进行评估，每个epoch窗长2s，步长0.5s，所以每个窗口间有75%的数据是重叠的
    % 每0.5s输出一组计算结果

    if (n>=4) && (sum(BSRmap(n-3:n))==0)       % 如果这段2s的数据，不包含任何抑制段
        x = eegHiPassFiltered((n+4)*nStride + (1:4*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 = EEG((n+4)*nStride + (1:4*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频率的均值（舍弃掉了部分小值，不知道为什么？？）
    % fft的结果是uV^2，10*log10转换成dB的单位
    % 之后按列计算均值（舍弃掉nan值），其实就是每个频点的近30s结果计算均值————相当于近30s信号的平均psd结果
    % nanmean是忽略掉nan值，计算其他剩下的平均值，除以的总数也是除去nan之外的
    % 之后找到中位数和最大值，并计算psd中，在中位数和最大值之间所有数字的均值（包含中位数和最大值）

    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 ,1) = meanBandPower(psd(thirtySec, :), 30, 47, 0.5) - midBandPower; % 镇静指数（30-47Hz - 11-20Hz）

    components(n, 2) = trimmean (10 * log10(vhighPowerConc./ wholePowerConc), 50);   % 全麻指数（40-47Hz/1-47Hz）
    % 剔除掉最大和最小值的各25%的占比（总数是50%）后，再求平均（trimmean函数用法写在下面了，但其实可以C代码里写死，因为数据点数是固定的）
    % trimmean(x,b)，其中b是%（b%为总数，上下各占一半 (b/2)%），且会剔除个数会向下取2的整数倍
    % 剔除点数的计算：需要剔除的总数据点数：数据长度 × b%。
    % 但这个计算结果会向下舍入到最接近的2的倍数（例如，计算结果为1，则舍入为0；计算结果为3，则舍入为2），
    % 以确保能从头部和尾部对称地剔除相同数量的数据点
    
    
    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值）
    %%% 改写自：components(n, 3) = meanBandPower(psd(thirtySec, :), 0.5, 4, 0.5) - midBandPower;
end
    
end

其中sawtoothDetector(RawEEG, nStride) 这个函数，我单独写了一个来研究他的计算原理：

function y = sawtoothDetector(eeg, nStride)         % 确定该脑电图片段是否包含强锯齿形k组分
    saw     = [zeros(1, nStride -5), 1:5]'; 
    saw     = (saw - mean(saw)) / std(saw, 1);      % 构建归一化锯齿波形波形
    r       = 1:(length(eeg) - length(saw)); 
    
    %% 计算移动方差
    v       = movvar(eeg, [0 length(saw)-1], 1);    % 滑动计算窗口内数据方差，每次移动一个点
    % [0 length(saw)-1]，窗口包含前0个点+当前点+后63个点(共64个点)
    % 最后的参数 1 表示计算的是总体方差，方差计算中的除数是窗口长度n(样本方差是除以n-1)
    % 相当于每移动一个点计算一次 var(eeg(1:64),1)
    % 相当于下面这段代码：
    % for i=1:length(eeg)
    %     if i<=length(eeg)-63
    %         v(i) = var(eeg(i:i+63),1);
    %     else
    %         v(i) = var(eeg(i:end),1);   % movvar函数当剩余数组长度不够时，会选择剩余的数据计算
    %     end
    % end
    
    %% 
    m       = ([conv(eeg, flipud(saw), 'valid') conv(eeg, saw, 'valid')] / length(saw)) .^2;    % 通过卷积匹配锯齿波和EEG波形
    % flipud函数将数组上下翻转
    % conv(x,b,'valid'); 仅计算的没有补零边缘的卷积部分
    y       = max([(v(r)>10).*m(r, 1)./v(r); (v(r)>10).*m(r, 2) ./v(r) ]) >0.63;        % 存在强匹配性的话，就返回true
end

上述过程中还有个 prctmean 函数，我也单独写出来方便理解和代码简洁：

function y = prctmean(x, lo, hi)
    v = prctile(x, [lo hi]);  
    % prctile——计算矩阵x中的[lo hi]所指代的分位数，比如[25 50]，就是计算1/4分位数和中位（50%）数

    y = mean(x(x>=v(1) & x<=v(2))); 
    % 先计算出来中位数和最大值，找出来在中位数和最大值之间的数据，计算其均值

end  % 计算百分比

5.加权计算最终结果

说实话，最后这段我也梳理不出来什么计算的逻辑，为什么要这样计算，这样计算代表了什么，不过反正是看懂代码并加以运用，就看懂就得了：

%% 加权计算得出最终的麻醉深度值结果
%%% 将成分1映射到logistic s曲线上，获得镇静评分
sedationScore = 104.4 - 49.4./(1+exp((components(:, 1)-(-13.9))/5.29));
% 改写自sedationScore   = scurve(components(:, 1), 104.4, 49.4, -13.9, 5.29);       % 在logistic s曲线上将成分1映射到镇静评分

%%% 将成分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进行插值
% 改写自 generalScore    = piecewise(components(:,2), [-60.89, -30], [-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;
% 改写自 generalScore    = generalScore + scurve(components(:, 2), 61.3, 72.6, -24.0, 3.55) .*(components(:,2)>= -30);  % 将组件2映射到一般分数，s曲线区域

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压缩和加权

上述计算过程中反复用到了piecewise这个函数，本来博主想全程少用函数来写，但是全是一样的操作还是封装函数方便点吧。具体函数内容如下：

function y = piecewise(x, xp, yp)  % 分段线性fn的响应
    xBound = min(max(x, xp(1)), xp(end));
    % 限制x矩阵的范围，为下面进行线性插值准备

    y = interp1(xp, yp, xBound); 
    % interp1()函数————插值用（默认使用线性插值）；
    % 按照x=[-60.89, -30]，y=[-40,43-1]，这个函数曲线，对刚刚限制后的xBound进行插值
end