Эта проблема связана с использованием find_peaks из scipy.signal для эффективного извлечения средней высоты пика из файлов данных. Я новичок в Python (3.7), поэтому я не уверен, что я написал свой код наиболее оптимальным образом с точки зрения скорости и качества кода.
У меня есть набор файлов измерений, содержащих одинмиллион точек данных (30 МБ) каждый. График этих данных представляет собой сигнал с пиками через равные промежутки времени и с шумом. Кроме того, базовая линия и амплитуда сигнала варьируются между частями сигнала. Я приложил изображение примера. Сигнал может быть намного менее чистым.
Моя цель - рассчитать среднюю высоту пиков для каждого файла. Для этого сначала я использую find_peaks, чтобы найти все пики. Затем я зацикливаюсь над каждым местоположением пика и обнаруживаю пик в небольшом интервале вокруг пика, чтобы убедиться, что получаю локальную высоту пика.
Затем я помещаю все эти высоты в массивы и вычисляю ихсреднее значение и стандартное отклонение после этого.
Вот базовая версия моего кода, она немного длинная, но я думаю, что это также может быть из-за того, что я делаю что-то не так.
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=peak_threshold, width=0)
#Calculate window size, make sure it is even
if round(len(time_dat)/len(peaks)) % 2 == 0:
points_per_window = int(len(time_dat)/len(peaks))
else:
points_per_window = int(len(time_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(0, len(peaks) - 2):
#Making a window around a peak_properties
sub_t = t_dat[peaks[i+1] - points_per_window // 2 : peaks[i+1] + points_per_window // 2]
sub_x = x_dat[peaks[i+1] - points_per_window // 2 : peaks[i+1] + points_per_window // 2]
#Detect the peaks (2 version, specific to the application I have)
prominence_threshold_singular = max(sub_x) - np.mean(sub_x)
sub_peaks_baseline, sub_peak_baseline_properties = find_peaks(sub_x, prominence=prominence_threshold_singular, distance=points_per_window - 1, width=0)
sub_peaks_height, sub_peak_height_properties = find_peaks(np.append(min(sub_x) - 1, sub_x), prominence=prominence_threshold_singular, distance=points_per_window - 1, width=0)
#Add the heights to the np arrays storing the heights
baseline_list[i] = sub_peak_baseline_properties["prominences"]
height_list[i] = sub_peak_height_properties["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))
Iожидается, что это будет выполняться довольно быстро, однако для его выполнения потребуется несколько десятков секунд. Это кажется мне медленным, поэтому я начал думать, что я могу проводить этот анализ неэффективно. Есть ли в моем коде некоторые аспекты, которые не выполнены должным образом? Или я недооценил, сколько времени должно занять выполнение этого вычисления?
Заранее спасибо всем, кто хотел бы взглянуть.