Multi-scale Peak and Trough Detection Optimised for Periodic and Quasi- Periodic Neuroscience Data

针对周期性和准周期性神经科学数据优化的多尺度峰谷检测
大概讲了讲bishop算法。主要是从Scholkmann算法改进的，用于检测颅内压的波峰。
生理信号中波峰和波谷的可靠检测对于医学和计算生物学的研究技术至关重要，也是许多信号处理任务的先决条件。准确峰值检测从简单的窗口阈值 window-thresholding 和小波变换技术 wavelet transform techniques 到隐马尔可夫模型HMM， k 均值聚类和基于熵的技术entropy-based techniques。
intracranial pressure 颅内压 (ICP) 波形的分析是一个挑战。适用于 ICP 峰值检测的算法必须适用于多尺度特征、随时间变化的波形形态和较差的信噪比

Scholkmann算法

Scholkmann算法，用于自动检测嘈杂的周期性和准周期性信号中的峰值，其驱动因素是他们需要在近红外光谱数据中找到峰值。 Scholkmann算法不需要任何参数，对高频和低频噪声具有相当的鲁棒性，能够准确检测准周期信号中的峰值（前提是振荡的最高频率小于或等于频率的四倍最低的信号）。该算法非常适合 ICP 波形峰值检测任务，但它受到低效的计算运行时间和相当大的内存需求的阻碍。

Original-Scholkmann 算法首先计算线性去趋势信号X，长度为 N 的局部最大值尺度图 ( local maxima scalogram LMS)，。
LMS是。
如果时间 t 属于{1 … N} ， S ∈ {1 …Smax} 是包含 0的局部最大矩阵，则包含 r + α 。r属于实数是均匀分布的随机变量，α也属于实数是常数。

LMC 可以更容易地可视化为Smax * N 矩阵，它标记每个比例（缩放级别）的最大值位置。
S为1 时，矩阵将在时间 t 编码一个局部最大值，如果信号值X大于相邻位置的信号值，以此类推S为2时，如果位置处的信号幅度 X 比 X − 2 和X + 2时的大，就编码一个局部最大值。

LMS 从第一个尺度（S= 1，最高或“最佳”分辨率）扩展到尺度S =N/2 − 1（大约是信号长度一半的“低分辨率”）。然而，它随后被裁剪为仅包括从 1 到 裁剪后的S 的尺度，裁剪后的S 是包含最大数量最大值的尺度。在算法的最后一步，计算跨尺度的列标准偏差，标准偏差为零的时间点确定最大值（峰值）的位置。

优化

通过反转原始信号并应用找峰算法来找到波谷。因此，可以以最小的额外计算成本。另外，用尺度S =N/2 − 1是多余的，可以在算法中参数化适当的比例上限，并使用特定领域的知识进行选择（也就是靠经验而不是广义搜素。。说的好像真的一样捏吗的，一个字解释就是靠玄学）
在峰值位置，LMS 的相应列将是一个零向量，可以通过线性搜索找到，从而无需计算列标准偏差和伪随机数。最后，LMS 应该只为每个信号计算一次并缓存以允许算法的后续运行在 O(n) 时间内完成。这在处理 ICP 波形时非常有用，因为可以使用该算法的递归应用来识别在线性时间内的 ICP 波形子峰P1到 p3

改进后

对标准尺度图技术的修改产生了一种具有线性计算复杂性的稳健算法，特别适合应对大型、嘈杂的生理数据集所带来的挑战。该算法可以通过优化可选参数针对特定应用进行调整，可以以最小的额外开销识别子波形特征，并且很容易适应在商品硬件上实时运行

附录了matlab代码

function [peaks,troughs,maximagram,minimagram] = new_peak_trough(data, varargin)
	N = length(data);
	if nargin == 2
		L = ceil(varargin{1}/2) - 1;
	else
		L = ceil(N/2) - 1;
	end
	%Detrend the data
	meanval = nanmean(data);
	data(isnan(data)) = meanval;
	data = detrend(data, 'linear');
	Mx = zeros(N, L);
	Mn = zeros(N, L);
	%Produce the local maxima scalogram
	for j=1:L
		k = j;
		for i=k+2:N-k+1
			if data(i-1) > data(i-k-1) && data(i-1) > data(i+k-1)
				Mx(i-1,j) = true;
			end
			if data(i-1) < data(i-k-1) && data(i-1) < data(i+k-1)
				Mn(i-1,j) = true;
			end
		end
	end
	maximagram = Mx;
	minimagram = Mn;
	%Form Y the column-wise count of where Mx is 0, a scale-dependent distribution of
	%local maxima. Find d, the scale with the most maxima (== most number
	%of zeros in row). Redimension Mx to contain only the first d scales
	Y = sum(Mx==true, 1);
	[~, d] = max(Y);
	Mx = Mx(:,1:d);
	%Form Y the column-wise count of where Mn is 0, a scale-dependent distribution of
	%local minima. Find d, the scale with the most minima (== most number
	%of zeros in row). Redimension Mn to contain only the first d scales
	Y = sum(Mn==true, 1);
	[~, d] = max(Y);
	Mn = Mn(:,1:d);
	%Form Zx and Zn the row-rise counts of Mx and Mn's non-zero elements.
	%Any row with a zero count contains entirley zeros, thus indicating
	%the presence of a peak or trough
	Zx = sum(Mx==false, 2);
	Zn = sum(Mn==false, 2);
	%Find all the zeros in Zx and Zn. The indices of the zero counts
	%correspond to the position of peaks and troughs respectively
	peaks = find(~Zx);
	troughs = find(~Zn)
end

自己改成python的就是应该这样（？

def bishop(x):
    signal = x
    N = len(signal)
    L = int(np.ceil(N / 2) - 1)

    # calculate max and min scalograms
    x = sp.detrend(signal, type="linear")
    m_max = np.full((L, N), False)
    m_min = np.full((L, N), False)

    # LMS matrices
    for k in range(1, L):  # scalogram scales
        for i in range(k + 2, N - k + 1):
            if x[i - 1] > x[i - k - 1] and x[i - 1] > x[i + k - 1]:
                m_max[k - 1, i - 1] = True
            if x[i - 1] < x[i - k - 1] and x[i - 1] < x[i + k - 1]:
                m_min[k - 1, i - 1] = True

    # find the scale with the max and min
    gamma_max = np.sum(m_max, axis=1)
    gamma_min = np.sum(m_min, axis=1)
    # find scale with max and min
    lambda_max = np.argmax(gamma_max)
    lambda_min = np.argmax(gamma_min)

    # use lambda to remove all elements of m for which k>lambda
    m_max = m_max[: (lambda_max + 1), :]
    m_min = m_min[: (lambda_min + 1), :]

    # find peaks and onsets
    m_max_sum = np.sum(m_max == False, axis=0)
    m_min_sum = np.sum(m_min == False, axis=0)
    peaks = np.asarray(np.where(m_max_sum == 0)).astype(int)
    onsets = np.asarray(np.where(m_min_sum == 0)).astype(int)
    
    return peaks