Использование find_peaks в цикле for для поиска пиковых высот из файлов данных для вычисления среднего - PullRequest
1 голос
/ 08 ноября 2019

Эта проблема связана с использованием find_peaks из scipy.signal для эффективного извлечения средней высоты пика из файлов данных. Я новичок в Python (3.7), поэтому я не уверен, что я написал свой код наиболее оптимальным образом с точки зрения скорости и качества кода.

У меня есть набор файлов измерений, содержащих одинмиллион точек данных (30 МБ) каждый. График этих данных представляет собой сигнал с пиками через равные промежутки времени и с шумом. Кроме того, базовая линия и амплитуда сигнала варьируются между частями сигнала. Я приложил изображение примера. Сигнал может быть намного менее чистым.

this is what the peaks might look like

Моя цель - рассчитать среднюю высоту пиков для каждого файла. Для этого сначала я использую 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ожидается, что это будет выполняться довольно быстро, однако для его выполнения потребуется несколько десятков секунд. Это кажется мне медленным, поэтому я начал думать, что я могу проводить этот анализ неэффективно. Есть ли в моем коде некоторые аспекты, которые не выполнены должным образом? Или я недооценил, сколько времени должно занять выполнение этого вычисления?

Заранее спасибо всем, кто хотел бы взглянуть.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...