阅读论文《Removal of Motion Artifacts in Photoplethysmograph Sensors during Intensive Exercise for Accurate Heart Rate Calculation Based on Frequency Estimation and Notch Filtering》

论文原文地址：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6696451/

2019年的文章，比较新。照例先看作者，Min Wang 1 , Zhe Li 2,, Qirui Zhang 1,3 and Guoxing Wang 1, 基本上是上交的几名学生，单字的名字同名会很多，而且他们应该是学生所以没查到。

前言

基波与谐波是什么？

和该振荡最长周期相等的正弦波分量称为基波，相应于这个周期的频率称为基本频率。频率等于基本频率的整倍数的正弦波分量称为谐波。谐波就是对周期性非正弦电量进行傅立叶级数分解，除了得到与基波频率相同的分量，还得到一系列大于基波频率的分量，这部分电量称为谐波。谐波频率与基波频率的比值（n=fn/f1） 称为谐波次数。谐波实际上是一种干扰。

overview

该论文提出的信号处理框架，第一步和前一篇论文一样，使用了0.4～4Hz之间（前一篇0.4～5Hz）的带通滤波器

这里提到 The 4th order 不知道滤波器的级数是什么意思？4次吗？

在 youtube 上看到四阶滤波器由两个二阶滤波器组成，如下图所示。此外，这个视频 是讲二阶带通滤波器的设计

details

（1）cascaded adaptive noise cancellation

目的：通过三轴加速度数据来动态的去除MA

input：Band-pass filtered PPG + acceleration signals

预处理：Z-score standardization 输入信号使得这它俩scale上能匹配

if there are multiple PPG signals at the input end of the cascaded ANC, the PPG signals will be averaged into a single PPG signal before the first adaptive noise canceller. 不明白为什么这里会有多个PPG信号进来？如果PPG信号时多个的话，加速度数据理应也是多个的

这里应该指的是收集的PPG信号可能是 mutiple channel 的，所以需要平均

采用 the least mean squares-Newton (LMS-Newton) algorithm，因为它比较容易收敛。该算法的公式如下：

$$ \begin{aligned} R^{-1}(\mathrm{k})=& \frac{1}{1-\alpha}\left[R^{-1}(k-1)-\frac{R^{-1}(k-1) X(k) X^{T}(k) R^{-1}(k-1)}{\frac{1-\alpha}{\alpha}+X^{T}(k) R^{-1}(\mathrm{k}-1) X(k)}\right] \\ & W(k+1)=W(k)+2 \mu e(k) R^{-1}(\mathrm{k}) X(k) \end{aligned} $$

x(k)：加速度信号

W(k)：weight vector of the finite impulse response (FIR) filter

e(k)：the output for de-noised PPG signal

α and μ：决定收敛的速度

FIR滤波器的学习链接1链接2，看了以下这个滤波器的diagram就明白了下图为什么有一个累加，确实是一个累加的过程，其中w(k)是d(n-k)的系数。现在的问题就是牛顿最小均方算法如何初始化了，这篇硕士论文看引言有研究这个，可能有所帮助

（2）heart rate frequency tracking (HRFT)

目的：初步估计 heart rate frequency （HRF）

先对信号做FFT得到频域上的频谱图(很多点做FFT能提高分辨率)，然后采用了[15]中的scheme

i：the i-th PPG signal sequence being processed

N_init：实验对象在实验开始阶段保持静止的时长

当i = 1时，f_PHR = the largest spectral magnitude of S_PPG within the human HRF range；

当i < N_init时，Δf = Δf0；

当i > N_init时，Δf = sum(max(|previous HRF estimation - current HRF estimation|)) + bias b

前一种情况，直接end；后两种情况，f_PHR = the largest spectral magnitude within a range of ±Δf around previous fHR and also within human HRF range

f_HRlow and f_HRhigh：人类心率典型的低和高边界

f_PHR：对心率的初步估计结果

（3）heart rate frequency correction (HRFC)

目的：纠正上一步估计可能存在的错误，给出 HRF 的最终估计

assumption：如果 ｜f_HR,pre - fHR| > 阈值Th0，很有可能是出错了

1.f_HR 初始化为 f_PHR

2.｜f_HR,pre - fHR| > 阈值Th0？第三步：do nothing, end

3.是否能找到 f_N？ 第四步：第五步

4.check f_ACC is within a small range of ±Δ around f_HR？ 很大可能性是错误的，第五步：第六步

5.check S_ppg(f_N) > th1 * S_ppg(f_HR)？f_HR = f_N, end：do nothing, end

6.check f_ACC is within a small range of ±Δ around f_N？第七步：第八步

7.check S_ppg(f_N) > th2 * S_ppg(f_HR)？do nothing, end：第八步

8.check S_ppg(f_N) > th3 * S_ppg(f_HR)？f_HR = f_N, end：do nothing, end

9.After final estimation of HRF, the HR value (60 f_HR bpm) for current PPG sequence is also calculated as an output of the proposed method.

f_N：the largest current PPG spectral peak between f_HR and f_PHR；如果能找到的话，f_N 很大可能性是心率的正确估计；如果没有找到

f_ACC：the frequency of the largest spectral magnitude of S_ACC

S_ACC：the magnitude spectrum of acceleration signal calculated from the average of the three-axis acceleration sequences

f_HR：对心率的最终估计结果

f_N 的存在是说这部分的心率很大程度是由 MA 造成的，不是准确的心率。这里很多细节没搞清楚，不如为什么th2~th3之间处理，大于或小于就不变？

（4）two tunable notch filters

目的：得到纯净的PPG信号

notch filters or comb filters can be used to eliminate unwanted components and keep only the PPG-related ones.

$$ \begin{array}{c}{\omega_{\mathrm{n}}=\frac{2 \pi f_{n}}{F_{s}}} \\ {H(z)=\frac{r^{2}-\left(1+r^{2}\right) \cos \left(\omega_{n}\right) z^{-1}+z^{-2}}{1-\left(1+r^{2}\right) \cos \left(\omega_{\mathrm{n}}\right) z^{-1}+r^{2} z^{-2}}}\end{array} $$

其中，第二个公式表示 frequency domain amplitude response of the notch filter

fn: the notch frequency Fs: signal sampling rate r: controls its bandwidth

滤波器的作用：抑制了 f_n 附近的信号，同时保持其他频率的分量接近其原始幅度。

每次迭代，把上一步得到的HRF和它的二次谐波频率作为 notch frequencies 来构建两个 notch filters，为什么这样设置？

because the main features of time domain PPG signal in one period are its main pulse and dicrotic pulse

使用两个 notch filters 对级联ANC级的输出进行滤波，在其中保留了MA的同时删除了PPG组件。通过最终从ANC去噪的PPG信号中减去陷波滤波输出，可以实现恢复的PPG信号。

results

目标：estimated HR value + the restored PPG signal

实验设计：

参数设置：

公开数据集

自建数据集

误差：前一篇论文中的error1，分子没有除

分析方法：相关系数（MA-removed PPG and reference PPG signal）/ Pearson correlation plot / Bland-Altman plot

related work

  1. ICA[6] 不满足假设无效[7]
  2. wavelet transform 小波变换[8] 需要凭经验设置阈值
  3. empirical mode decomposition (EMD)[9] 会受到间歇性（intermittency）噪声的影响，也叫 mode-mixing problem 模式混合问题
  4. combines Fourier series reconstruction and frequency-domain independent component analysis (FD-ICA) [10]
  5. combines multiple signal processing techniques and adaptive noise cancellation (ANC), which uses fast Fourier transform (FFT), singular value decomposition (SVD) and ICA [11]
  6. combine sparse signal reconstruction (SSR) and a special HR tracking scheme [13]
  7. ensemble empirical mode decomposition (EEMD) algorithm [14]
  8. combines the phase vocoder technique, a HR tracking and a smoothing stage [15]