论文阅读和分析：A Novel Time-Varying Spectral Filtering Algorithm for Reconstruction of Motion Artifact

论文阅读和分析：A Novel Time-Varying Spectral Filtering Algorithm for Reconstruction of Motion Artifact Corrupted Heart Rate Signals During Intense Physical Activities Using a Wearable Photoplethysmogram Sensor

运动伪迹的来源和难点：

1、MAs主要是由于环境光泄漏到PPG传感器表面和皮肤表面之间的间隙而引起；

2、运动引起的血流量变化是另一个MA来源；

3、去除运动伪迹的难点：MAs没有预先知道的的窄频带，而且它们的频谱经常与期望信号的频谱重叠；

数据集：

数据集用的是IEEE Signal Process和Chon lab的开源数据：

数据集信息：

数据集包含的运动类型：

1、activity involved walking or running on a treadmill for intervals of 0.5 min-1 min1 min-1 min-1 min-0.5 min with speeds of 1–2 km/h, 6–8 km/h, 12–15 km/h, 6–8 km/h, 12–15 km/h, 1–2 km/h, respectively. The subjects were asked to purposely move the hand with the wristband to generate motion artifacts.

2、activity included various forearm and upper arm exercise which are common arm motions (e.g., shaking hands, stretching, pushing objects, running, jumping, and push-ups).

3、 activity consisted of intensive forearm and upper arm movements (e.g., boxing).

4、activity involved 1 min rest, 1 min walking (3 mph), 1 min rest, 2 min jogging (5 mph), 1 min rest, 2 min running (7 mph), 1 min rest, 1 min arbitrary movement.

参考信号：The ECG-based reference HR was recorded in order to assess the performance of the algorithms being tested.

SpaMA算法总流程：

Stage 1. Time-Varying Spectral analysis
​ 1.1. Down sample the PPG and Accelerometer signal to 20 Hz
​ 1.2. Compute the power spectral density of both PPG and Accelerometers (0–10 Hz).

Stage 2. Spectral Filtering
​ 2.1. Assume HR to be in the frequency range of (0.5 Hz–3 Hz), this accounts for both low and high heart rates.
​ 2.2. The first highest three peaks and their corresponding frequencies in the PPG filtered spectrum are assumed to have HR information.
​ 2.3. Only the largest frequency peak of the accelerometers’ spectra is used for MA detection in stage 3

Stage 3. Motion Artifact Detection
​ 3.1. Compare the frequencies of the three peaks in the PPG spectrum with the frequency of the largest peak in the accelerometers’ spectra. If the first or second largest peaks in the PPG spectrum are similar to that of the accelerometers’ peaks, then motion artifact is present in the PPG.
​ 3.2. If motion artifact is detected from 3.1, then the corresponding frequency peak (usually the first or second largest peak) in the PPG spectrum should be discarded.

Stage 4. Heart Rate Tracking and Extraction from PPG Spectrum
​ Case (1): From 3.1—if the spectrum is corrupted by movement and only the first largest peak is corrupted, then the HR frequency should be the frequency of the second peak in the spectrum.
​ Case (2): From 3.1—if the spectrum is corrupted by movement and both the first and second largest peaks have similar frequencies to those of the accelerometers’ peaks, then the HR frequency should be the frequency of the third peak in the spectrum.
​ Case (3): Due to a gap between the pulse oximeter and a subject’s skin, the HR frequency cannot be extracted from the spectrum and in this case the previous HR frequency is used or for offline implementation a cubic spline interpolation can be applied to fill in the missing HR information.

Stage 5. PPG Signal Reconstruction
​ 5.1. The PPG signal is reconstructed by using the amplitude, frequency and phase information corresponding to the HR components (extracted in stage 4) that are calculated from the spectrum at each window

Heart Rate Variability Analysis
By using a sample-by-sample windowing strategy, HR can be extracted, from which dynamics of heart rate variability analysis can be obtained on the motion artifact-removed reconstructed HR time series.

注：实现第一步(PSD计算)的最简单方法是使用periodogram。然而，它的缺点是它是一个不一致的频谱估计器，具有高方差，并有泄漏效应。因此，一个主导谱峰可以导致在不应该有功率的频带中包含功率的估计频谱。然而，这两个问题都可以通过对原始信号进行下采样，然后通过设置大量频率点来使用足够小的频率步长来解决。因此，将信号从原始采样频率重新采样到它的1/4，然后应用频率分辨率为0.001的periodogram算法。

时域的分析：

对时频谱进行分析：

研究三个轴的时频谱和PPG的时频谱，绿色圈出来的是运动伪迹。

对频谱的峰值对应的频率进行分析：

画出PPG信号最高三个峰对应的心率的散点图：

画出参考心率曲线，观察曲线落在最高三个峰对应的心率的散点上的情况：根据加速度的频谱去除。

去除运动伪迹后，重构出时域图，再去计算HRV：
$$Re{c}_{Signal}\left(k\right)=A_{HR}\left(k\right)\times \sin \left(2\pi t\left(k\right)f_{HR}\left(k\right)+{\phi }_{HR}\left(k\right)\right)$$
式中，$k=1,\dots ,N$和$N$是信号样本数和窗口的总数；$A_{HR}(k)$和${\phi }_{HR}(k)$根据信号能量、在PSD上HR频率和对应HR频率的FFT矩阵的复数元素的相位角；

重构的结果：

对算法结果的分析：

画出参考心率和算法结果的对比图：

SpaMA算法性能对比：

其中的error计算方法如下：
$$Error\left(1\right)=\frac{1}{W}\sum _{k=1}^w\left|H{R}_{SpaMA}\left(k\right)-H{R}_{ref}\left(k\right)\right|$$

$$Error\left(2\right)=\frac{1}{W}\sum _{k=1}^w\frac{\left|H{R}_{SpuMA}\left(k\right)-H{R}_{ref}\left(k\right)\right|}{H{R}_{ref}\left(k\right)}\times 100\%$$

重构信号的PSD和ECG的PSD对比，以及计算相关度0.99

计算SpaMA重构信号的HRV和参考ECG信号的HRV，进行对比：