Python分析离散心率信号（下）
如何使用动态阈值，信号过滤和离群值检测来改善峰值检测。
一些理论和背景
到目前为止，一直在研究如何分析心率信号并从中提取最广泛使用的时域和频域度量。但是，使用的信号是理想的。现在考虑这个信号：
一个挑战！这是遇到的信号质量的另一个极端。老实说，当将传感器连接到手指上时（在0到4000之间），通过测量产生了该信号。在此之后，手指中的血管需要立即适应传感器的压缩（大约4000-5000），此后信号变得稳定。在大约7500、9000和12000时，用力将传感器移到手指上方以产生额外的噪音。对传感器本身能很好地抑制噪声感到惊讶。
尽管此信号已被手动“破坏”，但实际上可能会发生部分数据包含噪音或伪影的情况。传感器可能会移动，从而产生噪声，在测量过程中可能未正确连接或分离，参与者可能打喷嚏，移动或其任何噪声诱发因素都可能干扰。
在几个阶段中看到如何处理此问题：
§ 开始，将该信号传递给算法的结果；
§ 注意：采样频率；
§ 过滤信号以去除不想要的频率（噪声）；
§ 动态阈值提高检测精度；
§ 检测错误/错过的峰；
§ 消除错误并将RR信号重构为无错误。
开始启动
首先，让看看算法如何处理该信号。
(...) line 37, in detect_peaks
maximum = max(window)
ValueError: max() arg is an empty sequence
很好，事实并非如此。怎么了？
有人可能已经注意到了这一点，detect_peaks（）函数将在一种情况下中断：当心率信号从小于移动平均线变为变为等于而连续至少两个数据点都没有移动到其上方时，发生这种情况的最可能情况是信号在一段时间内下降到零。然后，该函数将跳到else语句，并尝试检测ROI中的峰值，但是由于信号从未上升到移动平均线之上，因此未标记ROI，因此max（window）抛出错误，因为空列表没有最大值。有时候是这样的因为是个白痴因为在循环中，没有考虑到两个信号的高度相等的情况，而这种情况恰好发生在信号的早期，将传感器固定在手指上：
如果之前在代码中没有注意到这一点，不用担心，一开始也没有。现在要更新detect_peaks（）函数：
def detect_peaks(dataset):
window = []
peaklist = []
listpos = 0
for datapoint in dataset.hart:
rollingmean = dataset.hart_rollingmean[listpos]
if (datapoint <= rollingmean) and (len(window) <= 1): #Here is the update in (datapoint <= rollingmean)
listpos += 1
elif (datapoint > rollingmean):
window.append(datapoint)
listpos += 1
else:
maximum = max(window)
beatposition = listpos - len(window) + (window.index(max(window)))
peaklist.append(beatposition)
window = []
listpos += 1
measures['peaklist'] = peaklist
measures['ybeat'] = [dataset.hart[x] for x in peaklist]
现在，还注释掉第19行，稍后再讨论。
def rolmean(dataset, hrw, fs):
mov_avg = pd.rolling_mean(dataset.hart, window=(hrw*fs))
avg_hr = (np.mean(dataset.hart))
mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
#mov_avg = [x*1.2 for x in mov_avg]
dataset['hart_rollingmean'] = mov_avg
然后使用峰值检测再次绘图。
该模块不再抛出错误并检测峰值，但这几乎不是想要的结果。让从降低噪声开始，看看是否可以改善信号。
采样频率
在进行信号过滤之前，让确定文件的采样率。前两部分的文件为100Hz，这部分呢？实际上，实际的记录频率可能随不同的设备或系统而变化。也可能与设备的额定值略有不同。所有计算得出的度量的准确性取决于准确了解采样率，因此进行检查非常重要。如果组合各种来源的数据，即使100Hz和101Hz的差异也可能导致明显不同的结果。通常在每条数据行中记录“世界时间”或“自记录开始以来的时间”，以便随后计算和验证采样率。
这样就可以很容易地计算出采样率：
#Simple way to get sample rate
sampletimer = [x for x in dataset.timer] #dataset.timer is a ms counter with start of recording at '0'
measures['fs'] = ((len(sampletimer) / sampletimer[-1])*1000) #Divide total length of dataset by last timer entry. This is in ms, so multiply by 1000 to get Hz value
#If your timer is a date time string, convert to UNIX timestamp to more easily calculate with, use something like this:
unix_time = []
for x in dataset.datetime:
dt = datetime.datetime.strptime(Datum, "%Y-%m-%d %H:%M:%S.%f")
unix_time.append(time.mktime(dt.timetuple()) + (dt.microsecond / 1000000.0))
measures['fs'] = (len(unix_time) / (unix_time[-1] - unix_time[0]))
此处提供的文件的采样率实际上是117Hz！现在可以使用确定的采样率自动计算度量。
请注意，这还不是全部，除了采样率之外，还应该检查数据以了解采样范围。是否所有样本均等地间隔开，例如，数据流中是否没有“间隙”或“跳过”？如果数据包含间隙和跳跃（如果最多只有几个样本），请在处理之前考虑对缺失值进行插值。如果采样率随时间变化很大，那么请使用其数据记录设备。
现在更加准确地了解了采样频率，接下来可以进行信号滤波。
信号过滤
遇到伪像或嘈杂的信号时，应该做的第一件事就是尝试过滤信号。为什么？因为在任何现实的测量情况下，信号（无论可能是什么）都将由信号部分和误差部分组成。这是因为不存在理想的传感器，因此将始终拾取来自要测量的源以外的干扰。电力线噪声是常见干扰的一个例子。在大多数国家/ 地区中，来自墙壁触点的交流电源的频率为50Hz（某些国家（例如美国）使用60Hz）。在许多原始ECG测量中也存在这种噪声。
一种常用的降低噪音的滤波器，其特点是对指定范围内的频率响应非常均匀。充当“频率门”；抑制超出指定截止范围的频率。这个临界点不是一条绝对线。意思是，如果将截止频率设置为例如5Hz，则5.1Hz的信号仍将通过滤波器，其幅度将稍微减小（或“音量”，如果更有意义）。另一方面，10Hz的信号只会非常弱地通过（如果有的话）。这也取决于过滤器的顺序。
还是很难理解？这样考虑：正在参加一个聚会，并且同时与两个人交谈，导致无法理解其中任何一个的情况。现在在和两个人之间放置一个过滤器。过滤器将减少人1的讲话音量而不会改变人2的音量。现在可以很好地理解人2了。这是滤波器的工作，但频率除外。
无论如何，开始编码，看看信号是否可以从滤波中受益。首先下载并打开，如果还没有这么做过的数据集，通过定义过滤器scipy.signal，过滤，最后绘制信号。假设正在使用上一部分中的代码，则按如下方式定义过滤器和绘图：
from scipy.signal import butter, lfilter #Import the extra module required
#Define the filter
def butter_lowpass(cutoff, fs, order=5):
nyq = 0.5 * fs #Nyquist frequeny is half the sampling frequency
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
def butter_lowpass_filter(data, cutoff, fs, order):
b, a = butter_lowpass(cutoff, fs, order=order)
y = lfilter(b, a, data)
return y
dataset = get_data('data2.csv')
dataset = dataset[6000:12000].reset_index(drop=True) #For visibility take a subselection of the entire signal from samples 6000 - 12000 (00:01:00 - 00:02:00)
filtered = butter_lowpass_filter(dataset.hart, 2.5, 100.0, 5)#filter the signal with a cutoff at 2.5Hz and a 5th order Butterworth filter
#Plot it
plt.subplot(211)
plt.plot(dataset.hart, color='Blue', alpha=0.5, label='Original Signal')
plt.legend(loc=4)
plt.subplot(212)
plt.plot(filtered, color='Red', label='Filtered Signal')
plt.ylim(200,800) #limit filtered signal to have same y-axis as original (filter response starts at 0 so otherwise the plot will be scaled)
plt.legend(loc=4)
plt.show()
现在，这个信号似乎并没有太大改善。如果仔细观察，线条会更平滑一些，但是一开始并没有很多高振幅，高频噪声。甚至可能会争辩说，滤波会稍微降低具有较低频率噪声的部分，因为在那儿也会抑制R峰。这是一个很好的例子，说明了为什么在决定过滤信号时总是要绘制信号。过滤信号会改变，这取决于决定此改变是否更好。
在另一个项目中使用过的这种嘈杂信号就是一个Butterworth滤波器非常有用的例子：
毫无疑问，这可以改善信号以进一步处理。
通过动态阈值提高检测精度
减少次峰标记不正确的第一种也是最明显的方法是提高用作阈值的移动平均值。但是提高到什么水平？对于许多不同的信号，这将是不同的。需要采取措施来帮助确定哪个阈值水平可能是最准确的。
一些简单的措施可以提供帮助，可以：
§ 查看信号长度并计算有多少个峰与预期的峰数；
§ 确定并使用RR间隔的标准偏差（将其称为RRSD）。
检测到的峰值数量可保留有价值的信息，但只能用来去掉明显不正确的阈值。根据应用程序（大多数人都坐着不动），可以过滤不太可能的bpm值。例如，过滤阈值导致平均bpm <30bpm和> 130bpm。在不同情况下（体育锻炼），阈值可以不同。
RRSD告诉有关所有RR间隔之间的差异如何分散的一些信息。通常，如果没有太多的噪声，则最低RRSD 不为零且可信的BPM值的检测将是最合适的。RRSD不能为零，因为这意味着没有心率变异性，这强烈表明了R峰检测中的错误。
这种简单的方法之所以有效，是因为缺少节拍会导致RR间隔大约是平均RR间隔的两倍，而错误地标记一个节拍会导致RR间隔最多是平均RR间隔的一半左右，但通常较小。两种情况都会导致RRSD升高。本质上，利用心率信号包含一个相当稳定的重复模式这一事实。
为了说明这一点，在数据集的子选择中绘制四个峰值检测回合，移动平均值分别提高0％，10％，25％和35％（从上到下）。
在倒数第二个图中，所有R峰均被正确检测到，没有任何东西被错误地标记为R峰。请注意，尽管BPM本身可能对所有四个图都有效（都不是绝对混乱），但RRSD强烈指出了最正确的图。说“最正确”是因为在某些情况下，无论如何设置阈值，都会残留一些错误，稍后将对此进行更多介绍。还要注意，与第三个图相比，最后一个图中缺少单个峰如何导致RRSD大大增加。
现在如何实施呢？不能简单地运行10.000个变化，每个变化均具有略高的移动平均线。除此之外，这将使算法整体性能严重受损，也将是非常多余的，因为许多迭代会产生相同的正确结果（而其许多迭代也会产生相同的错误结果！）。可能的解决方案是定期检查，然后评估最可能的RRSD和BPM对，如下所示：
def detect_peaks(dataset, ma_perc, fs): #Change the function to accept a moving average percentage 'ma_perc' argument
rolmean = [(x+((x/100)*ma_perc)) for x in dataset.hart_rollingmean] #Raise moving average with passed ma_perc
window = []
peaklist = []
listpos = 0
for datapoint in dataset.hart:
rollingmean = rolmean[listpos]
if (datapoint <= rollingmean) and (len(window) <= 1): #Here is the update in (datapoint <= rollingmean)
listpos += 1
elif (datapoint > rollingmean):
window.append(datapoint)
listpos += 1
else:
maximum = max(window)
beatposition = listpos - len(window) + (window.index(max(window)))
peaklist.append(beatposition)
window = []
listpos += 1
measures['peaklist'] = peaklist
measures['ybeat'] = [dataset.hart[x] for x in peaklist]
measures['rolmean'] = rolmean
calc_RR(dataset, fs)
measures['rrsd'] = np.std(measures['RR_list'])
def fit_peaks(dataset, fs):
ma_perc_list = [5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 150, 200, 300] #List with moving average raise percentages, make as detailed as you like but keep an eye on speed
rrsd = []
valid_ma = []
for x in ma_perc_list: #Detect peaks with all percentages, append results to list 'rrsd'
detect_peaks(dataset, x, fs)
bpm = ((len(measures['peaklist'])/(len(dataset.hart)/fs))*60)
rrsd.append([measures['rrsd'], bpm, x])
for x,y,z in rrsd: #Test list entries and select valid measures
if ((x > 1) and ((y > 30) and (y < 130))):
valid_ma.append([x, z])
measures['best'] = min(valid_ma, key = lambda t: t[0])[1] #Save the ma_perc for plotting purposes later on (not needed)
detect_peaks(dataset, min(valid_ma, key = lambda t: t[0])[1], fs) #Detect peaks with 'ma_perc' that goes with lowest rrsd
现在对数据集运行和情节（定时，整个算法到目前为止，包括装载在上i7-4790预处理约151ms，单核性能，所以仍然还是相当快的多线程将加快这一涨了不少。）：
这已经是好了很多。并没有丢弃任何正确的R峰，但是仍然存在一些不正确的检测，还有从0到大约5000的部分，其中几乎没有心率信号。将回到这个嘈杂的部分，并在以后了解如何检测和排除噪声。
现在，让从一开始就去掉嘈杂的部分，看看如何检测和过滤离群值。
查找错误检测的峰
最后要看的是如何检测和过滤异常峰位置。一种方法是利用心率逐渐而不是突然变化的事实。例如，bpm不会在单个拍子中从60bpm跳到120bpm，反之亦然。
同样，这意味着返回到RR间隔。如果在检测中跳过了一个峰，或者标记了两个峰而不是一个峰，那么所得的RR间隔将比平均间隔大得多或小得多。可以设置一个阈值并标记超过该阈值的间隔，类似于检测峰的方式。对于收集的数据来说，这就足够了。
然而，存在一种潜在的并发症。如果一次分析很长的信号，其中心率随时间变化很大，则可能导致错误的过滤。想象一下从60 bpm开始到180bpm结束时心率逐渐增加的信号。这意味着减少RR间隔的趋势是稳定的，这表明心率加快，而不是R峰检测中的错误。仅使用基于平均RR间隔的阈值，这将导致信号的第一部分和最后一部分被过滤。如果在数据中发生这种情况，可以先降低RR_list的趋势。使用scipy.signal，这很容易：
from scipy import signal
RR_list = measures['RR_list'] #First retrieve list from dictionary
RR_list_detrended = signal.detrend(RR_list, type='linear')
但是，如果信号包含一段大幅度增加的周期，然后心率又出现类似的大幅度下降，则将需要采用其方法。解决方案超出了本教程系列的范围，但是如果遇到此问题，则可能要使用截止频率非常低的高通滤波器。另一种方法是将信号分成较小的部分（以便将上升趋势和下降趋势分开），线性下降趋势并平均测量值。如果较小部分的长度不相等，请确保在平均之前权衡小节。
自然，不要使用去趋势化的RR信号计算任何量度，而仅将其用于检测峰标记中的错误。
返回异常值过滤。对于阈值，在实践中，发现在两端具有250-300ms频带的RR差均值的阈值水平很好。使用先前的代码并将数据集限制为[5000：15000]（目前要排除嘈杂的开始），应像这样实现：
RR_list = measures['RR_list'] #Get measures
peaklist = measures['peaklist']
ybeat = measures['ybeat']
upper_threshold = (np.mean(RR_list) + 300) #Set thresholds
lower_threshold = (np.mean(RR_list) - 300)
#detect outliers
cnt = 0
removed_beats = []
removed_beats_y = []
RR2 = []
while cnt < len(RR_list):
if (RR_list[cnt] < upper_threshold) and (RR_list[cnt] > lower_threshold):
RR2.append(RR_list[cnt])
cnt += 1
else:
removed_beats.append(peaklist[cnt])
removed_beats_y.append(ybeat[cnt])
cnt += 1
measures['RR_list_cor'] = RR2 #Append corrected RR-list to dictionary
plt.subplot(211)
plt.title('Marked Uncertain Peaks')
plt.plot(dataset.hart, color='blue', alpha=0.6, label='heart rate signal')
plt.plot(measures['rolmean'], color='green')
plt.scatter(measures['peaklist'], measures['ybeat'], color='green')
plt.scatter(removed_beats, removed_beats_y, color='red', label='Detection uncertain')
plt.legend(framealpha=0.6, loc=4)
plt.subplot(212)
plt.title("RR-intervals with thresholds")
plt.plot(RR_list)
plt.axhline(y=upper_threshold, color='red')
plt.axhline(y=lower_threshold, color='red')
plt.show()
结果：
似乎遇到了所有的小问题。结果是将所有正确检测到的R峰标记为绿色的图。不正确的标记为红色。生成的列表测量['RR_list_cor']的RR列表中没有属于错误峰的那些。
以后将研究如何标记和排除噪声段以及其一些优化。
汇总
整理所有代码，并更新一些函数以接受新变量并插入峰值抑制函数。
