python - 在循环中使用 find_peaks 从数据中查找峰值高度

Question

这个问题是关于使用scipy.signal.find_peaks从数据文件中有效提取平均峰高。我是 Python (3.7) 的初学者，因此我不确定在速度和代码质量方面我是否以最佳方式编写了代码。

我有一组测量文件，每个文件包含一百万个数据点 (30MB)。该数据的图表是一个具有规则间隔的峰值和噪声的信号。此外，信号的基线和幅度随信号的不同部分而变化。我附上了一个例子的图片。信号可能不太干净。

我的目标是计算每个文件的峰的平均高度。为了做到这一点，我首先使用 find_peaks 来定位所有峰值。然后，我循环遍历每个峰位置并在峰周围的小间隔内检测峰，以确保获得峰的局部高度。

然后我将所有这些高度放入 numpy 数组中，然后计算它们的平均值和标准差。

这是我的代码的准系统版本，它有点长，但我认为这也可能是因为我做错了什么。

import numpy as np
from scipy.signal import find_peaks

# Allocate empty lists for values
mean_heights = []
std_heights = []
mean_baselines = []
std_baselines = []
temperatures = []

# Loop over several files, read them in and process data
for file in file_list:
    temperatures.append(file)
    # Load in data from a file of 30 MB
    t_dat, x_dat = np.loadtxt(file, delimiter='\t', unpack=True)
    # Find all peaks in this file
    peaks, peak_properties = find_peaks(x_dat, prominence=prom, width=0)

    # Calculate window size, make sure it is even
    if round(len(t_dat)/len(peaks)) % 2 == 0:
        n_points = len(t_dat) // len(peaks)
    else:
        n_points = len(t_dat) // len(peaks) + 1

    t_slice = t_dat[-1] / len(t_dat)

    # Allocate np arrays for storing heights
    baseline_list = np.zeros(len(peaks) - 2)
    height_list = np.zeros(len(peaks) - 2)

    # Loop over all found peaks, and re-detect the peak in a window around the peak to be able
    # to detect its local height without triggering to a baseline far away
    for i in range(len(peaks) - 2):
        # Making a window around a peak_properties
        sub_t = t_dat[peaks[i+1] - n_points // 2: peaks[i+1] + n_points // 2]
        sub_x = x_dat[peaks[i+1] - n_points // 2: peaks[i+1] + n_points // 2]
        # Detect the peaks (2 version, specific to the application I have)
        h_min = max(sub_x) - np.mean(sub_x)
        _, baseline_props = find_peaks(
            sub_x, prominence=h_min, distance=n_points - 1, width=0)
        _, height_props = find_peaks(np.append(
            min(sub_x) - 1, sub_x), prominence=h_min, distance=n_points - 1, width=0)
        # Add the heights to the np arrays storing the heights
        baseline_list[i] = baseline_props["prominences"]
        height_list[i] = height_props["prominences"]

    # Fill lists with values, taking the stdev and mean of the np arrays with the heights
    mean_heights.append(np.mean(height_list))
    std_heights.append(np.std(height_list))
    mean_baselines.append(np.mean(baseline_list))
    std_baselines.append(np.std(baseline_list))

执行大约需要 30 秒。这是正常的还是太慢了？如果可以的话可以优化吗？

Answer 1

与此同时，我通过消除使用 Python 分析器发现的各种低效率问题，提高了速度。我将在此处列出按速度重要性排序的优化:

• 使用 pandas pd.read_csv() 进行 I/O 而不是 np.loadtxt() 会减少大约 90% 的运行时间。正如还提到的here ，这样可以节省很多时间。这意味着改变这个:

  t_dat, x_dat = np.loadtxt(file, delimiter='\t', unpack=True)
  

  对此:

  data = pd.read_csv(file, delimiter = "\t", names=["t_dat", "x_dat"])
  t_dat = data.values[:,0]
  x_dat = data.values[:,1]
  

• 删除多余的 len() 调用。我注意到 len() 被调用了很多次，然后注意到这是不必要的。改变这个:

  if round(len(t_dat) / len(peaks)) % 2 == 0:
      n_points = int(len(t_dat) / len(peaks))
  else:
      n_points = int(len(t_dat) / len(peaks) + 1)
  

  对此:

  n_points = round(len(t_dat) / len(peaks))
  if n_points % 2 != 0:
      n_points += 1
  

事实证明这也是一个重大改进。

• 最后，内置 Python 函数 min()、max() 和 sum() 使用了不成比例的计算时间部分(约 20%)。由于我已经在使用 numpy 数组，因此切换到这些函数的 numpy 等效项后，这部分的性能提升了 84%。这意味着例如将 max(sub_x) 更改为 sub_x.max()。

这些都是不相关的优化，我仍然认为对像我这样的 Python 初学者可能有用，而且它们确实有很大帮助。