数字开关鉴相算法在血氧仪的设想
m24h2022/06/04原创 电子技术 IP:上海

已经开源在github 如果使用STC-ISP来烧录 请仔细阅读XXXXXXXXX因为STC-ISP在EEPROM的地址分割上似乎有问题 需要手动纠正一下 

.................................

最近忽然想弄个小巧的血氧仪 目前只是想法 期待快更的 需要考虑我可能就不宣而停了(我未必有空或者有必要去实践 而且一直封控在家 快递也未正常 手上的笔记本还时好时坏 心情糟透了)

关键的元器件选型以小巧廉价为主:

  1. MAX30102作为传感器(约4 RMB)

  2. STC8g1k08A作为MPU(约1 RMB)

  3. 0.5寸OLED 48p*88p(约4 RMB)

血氧的算法可以使用美信公司提供的经验公式

而心率的算法目前有两种 一种是时域的 检测交流成分的过零或者峰谷周期 另一种是频域的 对信号做FFT/STFT然后找出峰值所在

STC8g1k08A虽然小巧便宜功能多如牛毛(1块钱 带ADC PWM PCA 比较器 无需晶振无需复位无需LDO 还支持USB下载 以至于我一直想找个机会用它)但是内存小 没有硬件乘法(更正一下:08A是带16位c51核心外的硬件乘除法器的 08没有A的话 8脚的是带的 其他脚数是不带的 而且操作复杂) 如果做FFT/STFT感觉有点力不从心 因此我这里设想了一个算法

  1. 对于输入的信号首先用卡尔曼滤波或者移动平均法去除高频噪声

  2. 移动平均法获取信号DC分量 AC幅值等低频信号 并检出其中的高频成分进行心率计算

  3. 设置160对变量 对应于40-200/分钟的每个频点 每个频点有两个变量 存储两个90度正交的相应频率的开关量对信号的乘积的带时间窗口加权的积分(真抱歉说得这么拗口 但是明白正交鉴相原理的应该会明白我在说什么 不就是I方加Q方再低通嘛 等仿真检验以后我再用算法语言描述吧)

  4. 对每对变量做平方和 找到最大值对应的频点就是心率

相比于传统的FFT算法 有下面好处

  1. 无需计算不需要的比如0-40/分钟处的频谱

  2. 除了最后一步求平方和 中间无需乘法计算 时间窗口加权需要的乘法也可以使用合适的权重通过位移等手段进行 而平方和的相互比较 在很大概率也可以避免乘法(当频谱明显时候 用绝对值之和也许就足够了)

  3. 计算是增量式积分式的 而不是帧式的 不光减轻了计算量和算法繁琐度(和复杂度的定义区分开) 还大大减轻了数据缓冲区对内存的要求 可以采用更长的窗口也没有带来更大压力

目前我还只是设想或者叫做“给可能需要的人出个点子” 而且网上也没找到类似的设计(也可能早有人做过我不知道) 

有时间的话 我先网上找一下MAX30102的原始采样数据 做一下算法的防真再说 虽然我现在已经在想成品大概什么样

IMG_20220604_145547.jpg

(是MAX30102 图上写错了懒得改了 我设想是按住后面的开关 直接接通电源并测量 松手关机)

如果利用起STC8g1k08A的剩余引脚 可以自己用红绿LED和光传感器取代MAX30102 可能还能压缩成本 但是标定和拟合公式对个人太麻烦了 总成本可能也就是10块和11块的区别

。。。。。。

补充说明一下 从数学上说 使用“开关”鉴相 结果中会是会混有奇次谐波的成分的 但是就这个应用而言 并不是要真正得到严谨的傅氏分解的频谱 只是想得到心跳频率 最在意的是 在得到的频谱上 心跳频率点是不是更能明显分辨出来

  1. 奇次谐波分量少 比如3次的混有1/3 5次的有1/5 平方后只有1/9 1/25 也许可以忽略

  2. 得到“积分和” 后可以进行修正 比如60频点的功率 减去1/9在180频点的功率 虽然不算严谨 但是有一定的实用性

  3. 注意我不是要从心跳中提取正弦波 我甚至考虑如果有时间的话 观察原始数据 构造更接近心跳波形的小波基函数代替“开关”方波 进行连续小波分析 可能频谱上心跳频率点会更突出

总之 有空先做个算法仿真 然后再看怎么改进算法 贴合实际

[修改于 2年4个月前 - 2022/07/10 07:14:40]

来自:电子信息 / 电子技术
16
6
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
m24h 作者
2年6个月前 修改于 2年6个月前 IP:上海
904762

从网上弄到了数据 今天用Excel做了算法仿真 感谢数据来源 是ghub上的glucose-PPG-data-set项目 据说来自某PPG开发板(不得不说 这个板子的原始数据也够漂亮  可能实际上没有那么好 或者已经滤波处理过)

IMG_20220605_155207.jpg

先看看原始数据 大概采样率是每秒8个点 (我移动平均点数选太大了 不过不影响分析) 

a1.png

如果采用包含多个周期的 时间窗口较长的FFT 结果如下图 

fft-long.png

如果采用只有2个波形的 短时间窗口的FFT 结果如下

fft-short.png

结果很明显也符合经验 长窗口带来好的分辨能力 但是非常大部分的频点都是不需要的 但又不能从蝶形分解运算量中忽略 该用的缓冲区内存也不能省(数量估计可以用结果频点数 乘窗口内需要容纳的最慢时候的心跳波包数得到 一般需要1000多个缓存点 还不包括蝶形分解再加一倍的缓存 才能满足包含几个波包的窗口)

而短窗口可以减少对不需要的频点的运算 但是结果的分辨率不高

而且 如果不是excel自带FFT分析 我不可能在上面用公式模拟 太烦了

但是正交鉴相就是用正交的两个波形去对采样值做简单的乘法和累加 可以容易的自写公式实现 开发起来也简单

下面是使用传统的正弦余弦波进行鉴相的结果 中间的曲面图揭示了心跳频点随着采样点增加 逐渐凸显的过程 

opd-sincos.png

计算简单 而且不需要对不在乎的频率进行计算 即使计算再长的时间窗口 也不需要多余的缓冲区内存 仅需要保留每个频点的积分值 32个采样点就有明显结果 256个采样点后 结果漂亮 分辨能力绝对够高 

如果采用开关方波进行鉴相 更不用乘法 只是加或者减 结果如下

opd-sign.png

与精密的正弦余弦鉴相对比 有所区别但是区别不大 对于找到心跳频率这一种民用产品 不需要更多的处理 也足够用了

再继续精简和去除剩下的乘法 用IQ两值的绝对值的和 代替需要计算平方和的矢量长 结果如下

opd-sign-no-multi.png

也是够用了 (对于每秒8个点而言 也许不必省下这不多的乘法 我可能会采用每秒20个点 应该也不用省下)

引用
评论
1
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
WernerPleischner
2年6个月前 IP:广东
904768

直接用数字PLL锁住 行不行

引用
评论(1)
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
三水合番
2年6个月前 IP:四川
904773

STC8g1k08A虽然小巧便宜功能多如牛毛…… 如果做FFT/STFT感觉有点力不从心

这个问题是切实存在的吗?STC8G 手册上说时钟频率能跑到30 MHz 多,处理比如“每秒8个点”的信号,会出现算力不足的情况吗?

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
黑体辐射
2年6个月前 修改于 2年6个月前 IP:新疆
904777
引用三水合番发表于3楼的内容
STC8g1k08A虽然小巧便宜功能多如牛毛…… 如果做FFT/STFT感觉有点力不从心这个问题是切...

🤔

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年6个月前 修改于 2年6个月前 IP:上海
904784
引用三水合番发表于3楼的内容
STC8g1k08A虽然小巧便宜功能多如牛毛…… 如果做FFT/STFT感觉有点力不从心这个问题是切...

算力并不算关键 虽然FFT不能用增量的方式每秒算8个点这么说 而是刷新一次算一套 就算慢一点也就是刷新慢一点 关键是内存小啊 才1kB

心跳量程128肯定不够的 用256吧 FFT结果只能用一半的(另一半是共轭幻象)就是得做512点

如果用8位精度 每层运算误差为1的话 算完这么多层 积累和扩大的误差甚至很容易超过量程了 要16位就需要2个字节 再加上这是复数 就是4个字节 这么算下来就2kB了 超过一倍了

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
warmonkey
2年5个月前 IP:广东
904850

这种信号用过中位触发去检测就足够了

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 IP:上海
904853
引用warmonkey发表于6楼的内容
这种信号用过中位触发去检测就足够了

看信号和心率波形了 如果是样本这样的波形没有问题 如果波形太复杂 过零法很难奏效

比如下面的波形 

IMG_20220609_132510.jpg

就属于频域法好用 时域法会很头疼的

再如下面的波形 

IMG_20220609_134000.jpg

如果不能完全精确地跟踪大起伏 也会影响结果的精度

还有像下面忽然的波形缺失(另一种情况是忽然的杂波)也很难处理

IMG_20220609_134321.jpg

开发上 时域法需要考虑波形各种不完美 打各种补丁 这让算法流程变得复杂不流畅 而频域法对抗以上几种情况都比较容易 算法写起来舒服


引用
评论
2
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 修改于 2年5个月前 IP:上海
904873

我已经下单了元件 并把原理图初版做出来了

sch.png

其实我的理想是用PW5100升压干电池到3.3或者5伏 因为烦了身边好多需要充电的小东西 长久遗忘又可能损伤电池 但是考虑有一点  血氧仪一个比较有用的用途 是去高原旅行用 干电池不如可以用充电宝充电方便 而且更小更轻对旅行者更友好

反正是小电池 所以采用电阻限流充电法 锂电池需要买带过压过充保护板的 并且安装前实测充满后是否自动停止

如果采用普通USB线只是充电 如果采用类似OTG线的ID接地的USB线 就可以USB下载 或者稳妥点接外面的TTL线

最右下角分压电路是为了更精确检测电池电压而复用 也许会精简掉 仅仅在检测MPU供电3.3V是否足够 如果不足就显示警告 实际上 锂电即使掉到3V 都能带警告地使用

此外..... 我在淘宝看的有那种可能是破产流出的成品血氧仪 低到6块钱一个 透射式的 ... 所以感觉动力不足

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 IP:上海
905039

修改了电池充电和供电线路的问题

sch.png

并且 把PCB做出来了 

pcb-photo.jpg

元件也到手了 有空就继续做下去

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 IP:上海
905078

好吧。。。我遇到一个坑

已经把裸板做好 但没安装按键 电池打算用一个4mm*11mm*15mm的 安装在显示屏下方 整体会非常小巧

IMG_20220620_230841.jpg

但尝试点亮OLED发现问题 这坑也是怪我粗心 用的8pin STC8G与非同系列的其他STC8G相比 SDA和SCL正好互换了!!! 我原本用的手册没有指出这个MCU的区别 下了一个比原来大了8倍的新手册 才发现这个区别 

IMG_20220620_232301.jpg

已经通过悬浮IC脚 重新交叉引线的方式解决 点屏和访问max30102都成功了

虽然如果最后成功 开源出的原理图和PCB都会是修正过的 不过原理图先更新一下

sch.png

引用
评论
1
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 修改于 2年5个月前 IP:上海
905166

基本已经实现 这大概是我目前做的最小巧的东西了 心跳和憋气测试有效果 (目前还没到考虑标定的那步)

1656061886800.jpg

从机器内通过串口导出的内部真实算出频谱结果如下 应该算具备足够的分辨能力了

pd-real.png

中间用采到的心跳数据 使用格兹尔算法也模拟了一遍 也有类似的效果

Goertzel.png

这个算法同样是非帧的 可以累进处理的 对内存要求也一样 但是多了cos查表和乘法 对运算误差有些敏感 而且还没有想到怎么进行收敛 (我现在开关鉴相算法中使用一阶低通来代替积分实现收敛)

这里把算法步骤放一下 运算相当简单 利用了整数的累加溢出后等效于取模的周期性

u16 idata g_pdi;   
u16 idata g_pdq;   
u16 idata g_pdf;   
i16 xdata g_hr_i[HR_MAX-HR_MIN+1]; // I信号数组 下标是心跳频率(从最小频率开始)
i16 xdata g_hr_q[HR_MAX-HR_MIN+1]; // I信号数组 下标是心跳频率(从最小频率开始)

// (65536*HR_MIN/60/SEN_RATE-6)=322, -6 is a fixed value for smallest error in whole heart beat rate scope 
#define PD_TADD  322
// (65536/(60*SEN_RATE))=11
#define PD_FADD  11

清零 g_hr_i[] g_hr_q[];
g_pdi=0;  // 0 degree
g_pdq=16384; // 90 degree

对每次采样结果(每秒100次):	{	
	pdi=g_pdi;
        pdq=g_pdq;
	依次对频率点(30-180/分): {
		if (pdi & 0x8000)
			g_hr_i[i]-=ac;
		else 
			g_hr_i[i]+=ac;

		if (pdq & 0x8000)
			g_hr_q[i]-=ac;
		else
			g_hr_q[i]+=ac;
              .....低通收敛一下....;
		pdi+=g_pdf; 
		pdq+=g_pdf;		
	}

	g_pdi+=PD_TADD;
	g_pdq+=PD_TADD;
	g_pdf+=PD_FADD;
}

我期间看了一下美信再mbed上的算法 其复杂性令人害怕 完全时域上的各种处理 找峰值 屏蔽干扰 鉴别等等等等.... 而我求心跳峰值的算法 借鉴了二极管检波求峰的原理和现有峰值自然衰退的方法

// max/min parameters decay
red_maxd -= (red_maxd - red32 +64)>>7;

// calculate new max/min
if (red32>red_maxd) {
	g_red_max += (red32 - g_red_max + 16) >> 5;
	red_maxd=red32;
} 		

 虽然极简 但是效果不错

max_maxd.png

目前程序离8k还差几百字节(已经包括了接近检测 光功率自动调整 串口输出和命令下载等所有的杂项功能) 还没有用任何浮点运算 (否则就只能把字库和查表放到eeprom里面省出足够空间) 可以算可用的初版了 下面等电池到货 设计一下外壳就差不多了 可能还想再优化一下鉴相算法 尤其是低通收敛这块 看看能不能让结果更平滑点

....................................

试了一下内建余弦表 使用余弦鉴相代替开关鉴相 得到心跳频谱图如下

hr_cos.png

有所改观但是区别不大 但是算力恰恰够 我的显示刷新是在空闲时候进行(并未打开中断而是轮询) 使用余弦鉴相后 刷新率下降到了2秒一次 而使用开关鉴相 刷新率基本符合0.6秒一次 代码仅仅改动如下 可以略微优化但是余量不大 现在数据采集率100Hz 如果降到50Hz可能会轻松点  但是衡量一下 还是保持采集率好处更大

i16 cos(u8 x)
{
	if (x<=64)  return COS_LUT[x];
	if (x<=128) return -COS_LUT[128-x];
	if (x<=192) return -COS_LUT[x-128];
	return COS_LUT[(u8)-(i8)x];
}

#ifdef USE_COS
		g_hr_i[i]+=(i16)(((i32)ac*(i32)cos((u8)(pdi>>8)))>>16);
#else		
		if (pdi & 0x8000)
			g_hr_i[i]-=ac;
		else 
			g_hr_i[i]+=ac;
#endif

..................................

还是优化了 减少了32位乘法 表中直接使用8位数减少位移 

i16 cos(i8 x)
{
	if (x<0) x=-x;
	if (x>64) return -COS_LUT[128-(u8)x];
  return COS_LUT[x];	
}
...............
.............
#ifdef USE_COS
		g_hr_q[i]+=(ac*(cos(*(i8*)&pdq)))>>7;
#else		

刷新率到了1秒1次 再三衡量 还是用这个罢

所以 题目虽然是开关鉴相 但是 最后还是做了余弦鉴相 .... 不过 也不妨碍开关鉴相好用可用

引用
评论
1
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年5个月前 修改于 2年4个月前 IP:上海
905472

做了个壳子 这下不怕OLED排线被扯断了(已经扯断过了 而且修复失败 因为焊几根短铜线不难 难的是焊接时受热又脱开 用阻焊UV胶固定也不好用 只能再下单买了一块)

不过如果再设计一次 我应该长宽方向换一下 不做长条做宽条状 太小了感觉手指按得不舒服

IMG_20220707_182412.jpg

目前感觉还可以优化一下 改用sdcc来编译(避免版权问题 )周末如果弄好了再开源

。。。。。。。。。

本来的想法是按下开关开始测量 食指正好压在后面的传感器 但是现在发现一个问题

食指会过于用力压在传感器和外壳上 影响了正常的脉搏血流 原始数据很清晰地反映了这一点 首先起伏变小 原理类似血压计 其次是刚开始测量 手指刚受压 需要很长时间血流和脉搏才进入稳定的状态

早知道就用拨动开关接通电源 传感器放在上方 手指轻轻搭上去就好 但是来不及了

这是个教训 也希望给别人一个经验

我这个不打算再改了 只能使用时候注意技巧了

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
irly01
2年0个月前 IP:安徽
910341

楼主您好,那个通过开关或正余弦正交鉴相算法  , 经过256次采样后积分的计算结果的峰值应该不是频率值吧,而是相位值吧。

计算频率值是不是要把那个图求一次导数才行?还是峰值位置减一进行近似计算?我也不知道我说的对不对。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
2年0个月前 IP:上海
910354
引用irly01发表于13楼的内容
楼主您好,那个通过开关或正余弦正交鉴相算法  , 经过256次采样后积分的计算结果的峰值应该不是频率...

你再算算 I或者Q都和相位有关 但两者平方和就消去相位了

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
zjsx8192
1个月14天前 IP:广东
938234

请问楼主,算法中把USE_COS这个定义去掉,是不是就是开关鉴相的算法?

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
m24h作者
1个月14天前 IP:上海
938235
引用zjsx8192发表于15楼的内容
请问楼主,算法中把USE_COS这个定义去掉,是不是就是开关鉴相的算法?

哈哈 不太记得了 应该是吧

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

所属专业
所属分类
上级专业
同级专业
m24h
进士 学者 机友
文章
52
回复
895
学术分
1
2020/01/22注册,18分28秒前活动

个人开源项目: XXXXXXXXXXXXXX

主体类型:个人
所属领域:无
认证方式:手机号
IP归属地:上海
文件下载
加载中...
{{errorInfo}}
{{downloadWarning}}
你在 {{downloadTime}} 下载过当前文件。
文件名称:{{resource.defaultFile.name}}
下载次数:{{resource.hits}}
上传用户:{{uploader.username}}
所需积分:{{costScores}},{{holdScores}}下载当前附件免费{{description}}
积分不足,去充值
文件已丢失

当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}
视频暂不能访问,请登录试试
仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。
音频暂不能访问,请登录试试
支持的图片格式:jpg, jpeg, png
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
收藏
取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}