基于rPPG的心率测算程序设计报告[公开版]

观前提醒：笔者非医学生，非专业学习图像处理，如有任何表述不明，有失偏驳，请在评论区提出意见。

现如今，大众越来越重视个人健康或者运动，在此基础上国内各种可穿戴设备也发展飞速，小米手环、华为Watch、OPPOwatch等设备在国内外广受好评，各种手环/手表最基础的身体健康监测使用PPG传感器（Photoplethysmography 光电容积脉搏波描迹法），体现心率同时进行健康评测（例如心脏早搏，房颤等等，PPG并非心电图，所反应心脏信息有限，一般是各种厂商联合医疗研究机构使用AI进行判断，例如华为心脏健康研究，OPPO心脏健康研究），以及通过相关算法完成各种多维数据（体能等）的测算。

在文章开始前，我们需要先知道PPG实现原理，想必大家小学都知道，心脏跳动是将血液泵送至全身，底层原理类似于老中医把脉，通过脉搏来获取心脏跳动。但是可穿戴设备无法像老中医一样自行寻找到脉搏明显起伏的地方，也不可能像血压计一样通过充气紧紧束缚手腕从而获取心脏跳动（PS：事实上早已有了充气测血压的可穿戴设备，但其主要作用是为亚健康人群提供随时可用的血压计，且作为可穿戴设备舒适性不够）。

这里我用一篇我认为解释的很浅显易懂的图来描述：

同样的光强的绿光，射到大杯和小杯红色的水里，由于红色的水对于绿色吸收率高，水量不一样反射回去的光强差别较大。
而同样光强的红光，射到大杯和小杯红色的水里，由于红色的水对于红色吸收率低，水量不一样反射回去的光强差别较小。

现在我们想象血管就是上面的水，随着心脏搏动，水量变多变少。

如果入射光是绿光，心脏搏动导致的反射光变化较大；
如果入射光是红光，心脏搏动导致的反射光变化较小。

换句话说，用绿光得到的信号会比用红光得到的信号变化幅度要大。[2]

在上述图示中我们可以很简单的理解到心率传感器使用绿光的原因，接下来，我们再来了解一下心跳规律。

心脏一次规则的跳动可以用以下心电图表示：

如果你体检时做心电图，心电图报告正是由这一个一个单位图示波形组成，不过体检心电图使用的是多导联，不同导联所展示心电图信息与图示信息可能会有差异，能够反应更多的心率信息，同时，由于肌肉电信号与仪器影响，实际心电图并不会有这样规范，当然，我不是医学生，并不能做出更多解释。

心电图其实就是相当于在你身体加了个电流计，测量电流随时间的变化，实际上更为复杂，需要处理各种电信号及滤波。

至于各种波我这里不做出过多讲解，您可以参考以下图示就能理解心脏跳动规律了。

你可以看见心电图R波非常高，实际上PPG传感器就是监测这个R波，所以叫脉搏波描迹法，再用R波单位时间出现频率就可获得人体心率（以下图示由于超出图示范围在顶峰显示为平波，实际上是一个类正弦波形，同时，心电图使用单导联以及肌肉电信号干扰，并不能表示图示中的完美示例，由于电信号速度远大于血液传输速度，因此脉搏波相比心电信号有明显滞后性，值得一提的是，由此可以推算出血液传输速率，从而进行各种健康评估，比如血管健康）

在以上图示中可以很明显的看出来PPG已经能够估算出心脏跳动频率，但是实际上PPG并不能如此完美的体现出人体脉搏波动，因为毕竟是光电传感器，会受各种干扰，比如自然光，或者测试者的运动，因此需要各种滤波算法，才能较为准确的体现心率信号。

可能此时会有读者有疑问，你这篇文章不是说的rPPG吗，为什么你这里一直讲的PPG。本质上是一个东西，PPG指Photoplethysmography，rPPG则指Remote Photoplethysmography ,多了个Remote，大家都知道可穿戴设备是接触式的，在买来时APP或说明书都会告诉你如何佩戴测量心率最准确，这里的rPPG则是远程实现心率测量。

那么，远程依靠什么远程，测量什么呢？
在众多研究中，大家基本上均采用额头和面颊来测得心率数据，一方面，除了监控外，现在数码设备的摄像头一般面朝自己脸部，另一方面额头相比其他人体裸露在外部的部位能更加明显的体现人体心脏跳动，不知道你有没有这样的经历，当压力大或者心情不好时能很明显的感受到你的太阳穴脉搏跳动，这种微小的变动会影响光强大小的改变，那么我只需要获取到额头平均光强也就可以得到P波了。

在此基础上，我们可以通过人脸识别定位出人脸，从比例上再定位额头，这里使用golang opencv（这里使用Golang实现，而不使用python以及MatLAB这种专注于图像处理的编程语言及工具，起初是想做成网页端，让每个人都可以测试，同时增加样本数据，更有助于我学习并优化算法，但是考虑到隐私及成本影响因此停止实现）库实现人脸识别定位：

// 忽略前面的代码
	classifierFace := gocv.NewCascadeClassifier()
	defer classifierFace.Close()
	if !classifierFace.Load("haar/haarcascade_frontalface_default.xml") {
		log.Fatalf("加载人脸分类器文件时出错")
	}

// 忽略中间处理
			faces := classifierFace.DetectMultiScale(img)
			if len(faces) == 0 {
				continue
			}

			maxFace = faces[0]
			for _, r := range faces {
				if r.Dx()*r.Dy() > maxFace.Dx()*maxFace.Dy() {
					maxFace = r
				}
			}

以上代码中使用了opencv haar库实现了简单且精准的haar级别库来定位人脸，并引入一个for循环来找到图像中最大的人脸，至于为什么需要引入for循环找最大人脸，是为了尽可能的排除掉非测试者走动经过摄像头拍摄区域影响。

在接下来的流程中还需要定位额头

			forehead = maxFace
			forehead.Min.Y = maxFace.Min.Y + maxFace.Dy()/10
			forehead.Max.Y = maxFace.Min.Y + maxFace.Dy()*2/10
			forehead.Min.X = maxFace.Min.X + maxFace.Dx()/2 - maxFace.Dx()/10
			forehead.Max.X = maxFace.Min.X + maxFace.Dx()/2 + maxFace.Dx()/10

			gocv.Rectangle(&img, forehead, color.RGBA{255, 0, 0, 0}, 1)
			gocv.Rectangle(&img, maxFace, color.RGBA{0, 255, 0, 0}, 1)

以上代码中定位额头我比较简单粗暴，根据选择的人脸，计算前额的位置。这里假设前额占据人脸高度的1/10，宽度为人脸宽度的1/5，位于人脸上部中间位置。但是这样并不能直接使用，因为由程序定位的人脸会时刻发生小幅度变化，在实际测量中我们应该尽可能排除掉程序影响，因此为了方便我通过代码实现了输入“5”实现固定额头位置，不再实时监测人脸位置，再按下“5”解除释放，这里就不放出代码了，仅需一个Flag+if else就可以实现。
额头如上图红框所示。

实时人脸监测活动就像这样。

实际上，在众多研究项目及论文中均使用面颊+额头计算，但是笔者算法太菜，就不过多纠结了。

前面提到，最适合的测量心率的光线颜色是绿光，那么我们将使用RGB颜色系统中的G通道表示，也有研究中使用HSV颜色模型进行测算。

那么我们先看看额头RGB颜色系统各个分量光强随时间变化的曲线。

各个通道平均光强可使用以下示例代码实现：

		greenChannel := gocv.NewMat()
		gocv.ExtractChannel(img.Region(forehead), &greenChannel, 1)

		greenChannelData := greenChannel.ToBytes()
		for i := 0; i < len(greenChannelData); i++ {
			greenChannelSum += float64(greenChannelData[i])
		}
		greenChannelCount = len(greenChannelData)

		greenChannelAvg := greenChannelSum / float64(greenChannelCount)

		blueChannel := gocv.NewMat()
		gocv.ExtractChannel(img.Region(forehead), &blueChannel, 0)
		blueChannelData := blueChannel.ToBytes()
		blueChannelSum := 0
		for i := 0; i < len(blueChannelData); i++ {
			blueChannelSum += int(blueChannelData[i])
		}
		blueChannelAvg := float64(blueChannelSum) / float64(len(blueChannelData))

		redChannel := gocv.NewMat()
		gocv.ExtractChannel(img.Region(forehead), &redChannel, 2)
		redChannelData := redChannel.ToBytes()
		redChannelSum := 0
		for i := 0; i < len(redChannelData); i++ {
			redChannelSum += int(redChannelData[i])
		}
		redChannelAvg := float64(redChannelSum) / float64(len(redChannelData))

（注，几个比较明显的波动是由于头部活动导致）

根据以上图可以看见，红光平均光强更高，这是由于人体反射红光水平更强，而且波动幅度太大，但是蓝光平均光强更低，而且人造光源（屏幕等）通常包含部分蓝色杂光，较容易影响人体心率数据采集。

通过观察波形我们可以看见会有一些小起伏，这些起伏中就可能包含脉搏波数据。

作为对比，下面放出稳定后的摄像头识别白墙区域RGB平均光强变化波形图

前面的起伏是由于程序写好的先人脸锁定后再开始记录，可以看见波动后几乎为直线，不过由于测试时为晚上，人造光源存在部分频闪，因此会出现些微规律波动（亦或是摄像头规格限制），但不明显，就像我前面提到的，人造光蓝光较多，且人造光源存在着频闪，导致蓝光光强随时间变化波动最大。

接下来，我们要进行滤波。根据相关文献，正常成年人静息心率为50-100左右，也就是心脏每分钟跳动次数。[3]

但是在部分情况下，例如运动较多人群或者刚运动完，心率会更低或更高，因此我们需要筛选掉一些数据，我将所需BPM设定为45-240，也就是心率范围0.75-4Hz的数据，代码如下所示（结合注释查看）

var BPMMin float64 = 45 // 所需最小心率
var BPMMax float64 = 240 // 所需最大心率
func calculateBPM(data []float64, times []float64) float64 {
	if len(data) < 2 {
		return 0
	}

	timeDiff := times[len(times)-1] - times[0]
	if timeDiff <= 0 {
		return 0
	}

	samplesPerSecond := float64(len(data)) / timeDiff

	// 引入带通滤波器以减少噪声
	filteredData := bandPassFilter(data, samplesPerSecond)

	// 应用汉宁窗以减少频谱泄漏
	windowedData := make([]float64, len(filteredData))
	for i, d := range filteredData {
		windowedData[i] = d * 0.5 * (1 - math.Cos(2*math.Pi*float64(i)/float64(len(filteredData)-1)))
	}

	// 将数据转换为复数
	complexData := make([]complex128, len(windowedData))
	for i, d := range windowedData {
		complexData[i] = complex(d, 0)
	}

	// 应用FFT
	fftResult := fft.FFT(complexData)

	// 计算每个频率点代表的实际频率（单位：Hz）
	freqs := make([]float64, len(data)/2+1)
	for i := range freqs {
		freqs[i] = float64(i) * samplesPerSecond / float64(len(data))
	}

	maxMag := 0.0
	maxIdx := 0
	for i, c := range fftResult[:len(data)/2+1] {
		mag := cmplx.Abs(c)
		if freqs[i] >= BPMMin/60 && freqs[i] <= BPMMax/60 {
			if mag > maxMag {
				maxMag = mag
				maxIdx = i
			}
		}
	}

	bpm := freqs[maxIdx] * 60 // 转换为每分钟的心跳次数

	return bpm
}

func bandPassFilter(data []float64, samplesPerSecond float64) []float64 {
	// 创建一个带通滤波器，将频率限制在0.75-4 Hz之间
	lowFreq := BPMMin / 60.0
	highFreq := BPMMax / 60.0

	// 创建一个带通滤波器
	filteredData := make([]float64, len(data))
	for i := range data {
		if i == 0 {
			filteredData[i] = data[i]
			continue
		}

		// 计算当前频率
		freq := float64(i) * samplesPerSecond / float64(len(data))

		// 计算带通滤波器的系数
		lowPassCoef := 1 / (1 + math.Pow(freq/lowFreq, 4))
		highPassCoef := 1 / (1 + math.Pow(highFreq/freq, 4))

		// 应用带通滤波器
		filteredData[i] = lowPassCoef * highPassCoef * data[i]
	}

	return filteredData
}

• calculateBPM函数接受时间序列的光强度数据，首先通过FFT转换得到频域数据。
• 在一定的频率范围内（对应人的心率范围），找到最大幅度的频率，并将其转换为每分钟的心率（bpm）。
• FFT（快速傅里叶变换）是信号处理中一种重要的工具，用于将信号从时域转换到频域，这里使用现成的库处理。
• bandPassFilter带通滤波器用于筛选感兴趣信号，后面会提到。这里算法过于简单，Golang本身不是用来干这行的，暂时还没找到相关数学库，因此需要自己实现，这里暂时仅简单处理。
• 汉宁窗（Hanning window）是一种窗函数，用于对信号进行加窗处理，以减少频谱泄漏（spectral leakage）和提高频谱分辨率。频谱泄漏是指当信号进行离散傅里叶变换（DFT，此处通过FFT实现）时，由于信号的周期性假设与实际信号的非周期性不匹配，导致频谱能量泄漏到其他频率 bin 中。这会使得频谱分析结果模糊，难以准确识别信号的频率成分。汉宁窗通过减少信号两端的幅值，使得信号在时域上平滑地过渡到零，从而减少频谱泄漏。具体来说，汉宁窗在每个数据点上乘以一个特定的系数，这些系数是根据汉宁窗函数计算得出的，使得数据序列的两端逐渐衰减到零。

转换到频域后，信号被表示为一系列频率和对应的幅度。每一个频率成分的幅度表明了该频率在原始时域信号中的强度。在心率测量的上下文中，我们特别关注那些与心脏跳动频率相对应的频率成分。

设定搜索范围：人的正常心率大约在每分钟45次（0.75 Hz）到240次（4 Hz）之间，这对应于频率的0.75 Hz到4 Hz。FFT分析结果中，只有这个频率范围内的成分才是我们关注的目标。
寻找最大幅度的频率：在上述范围内，程序寻找幅度最大的频率点，因为最大幅度表明了信号中最显著的周期性成分，即心跳引起的周期性变化最强烈的部分。

一旦识别出心跳对应的主要频率，将这个频率乘以60（每分钟的秒数），就可以得到心率的估计值，单位是“每分钟的心跳次数（bpm）”，在实际程序设计中需要优化，我这里非常简单并没有实现采样时间限制，实际上采样时间应该至少5S以获取更为准确的心率。

大起伏为人脸活动造成的影响，小起伏可以视为心脏跳动时的P波（实际上由于算法存在太多问题，会有干扰信号），至此，此代码就告一段落了，对于滤波相关，严格意义上需要对图像进行降噪处理并进行更复杂的滤波。但是在ECG与实验结果对比情况下已经知足了，图一为瞬时心率，图二为30S平均心率，笔者睡眠不好，有窦性心律不齐，某种意义上来讲已经符合事实了。

不过在实际应用中rPPG范围太窄已经淹没于历史的洪流中，就算是现在一线厂商的PPG传感器及算法也不能做到医疗级别的监测，rPPG由于受各种环境影响，变量更多，其能反应身体健康状态信息较少，也仅限于学术研究中，不过相信这项技术在未来更高精度的摄像头发展中，这项技术将被重新使用，在未来为智能设备用户提供更准确的身体健康反馈。

引用文献：
[1] 维基百科 光体积变化描记图法
[2] 知乎 绿or红？心率监测使用哪种光比较好？
[3] 医学百科 心率正常值
参考文献:
[1] Akselrod S, Gordon D, Ubel F A, Shannon D C, Berger A C, Cohen R J. 1981. Power spectrum analysis of heart rate fluctuation:a quantitative probe of beat-to-beat cardiovascular control. Science, 213(4504): 220-222 [DOI:10.1126/science.6166045]
[2] Niu X S, Han H, Shan S G . 2020. Remote photoplethysmography-based physiological measurement: a survey. Journal of Image and Graphics, 25(11): 2321-2336. (牛雪松, 韩琥, 山世光. 2020. 基于rPPG的生理指标测量方法综述. 中国图象图形学报, 25(11): 2321-2336.) [DOI: 10.11834/jig.200341]
[3] Chen X, Cheng J, Song R C, Liu Y, Ward R, Wang Z J. 2019. Video-based heart rate measurement:recent advances and future prospects. IEEE Transactions on Instrumentation and Measurement, 68(10): 3600-3615 [DOI:10.1109/TIM.2018.2879706]
[4] De Haan, Gerard, and Arno Van Leest. "Improved motion robustness of remote-PPG by using the blood volume pulse signature." Physiological measurement 35.9 (2014): 1913.
[5] Wang, Wenjin, et al. "Algorithmic principles of remote PPG." IEEE Transactions on Biomedical Engineering 64.7 (2016): 1479-1491.