Площадь под кривой для набора данных и двух наборов базовых линий

Код ниже находит область под кривой для прямой развертки.

Каждый цикл сначала разделяется на прямую развертку напряжения и обратную развертку. Затем вы можете предоставить данные форвардного цикла в calculate_fwd_cycle_peak_area(). Эта функция находит пик, рассматривая особенности градиента сигнала. Функция также отображает график, показывающий расположение различных «ориентиров», которые используются для обозначения границ пиков. Площадь вычисляется с использованием трапециевидного интегрирования numpy, и возвращается результат.

Функция настроена с различными значениями по умолчанию, которые, как я обнаружил, хорошо работают с предоставленными вами данными (все циклы пересылки в примерах 1, 2 и 3).

Пример результата:

Использование:

#Load data
data = load_numeric_data(f'/mnt/c/dev/CV_sample1.txt')

#Define the start potentials of the forward and reverse sweeps
fwd_sweep_start_potential = -1.3
rev_sweep_start_potential = 0.7

#Find the start indices
fwd_cycle_starts = find_cycle_starts(data, fwd_sweep_start_potential)
rev_cycle_starts = find_cycle_starts(data, rev_sweep_start_potential)

#Get the cycles into a list
fwd_cycles_list, rev_cycles_list = split_on_cycles(data, fwd_cycle_starts, rev_cycle_starts)

#Supply each cycle in turn to the area function
for i, cycle in enumerate(fwd_cycles_list):
    area = calculate_fwd_cycle_peak_area(cycle)
    print('Peak area of cycle', i, 'is:', area)

Полный код приведен ниже.

import numpy as np
from matplotlib import pyplot as plt
import pandas as pd

def load_numeric_data(filepath):
    data_start_marker = "Potential/V, Current/A"
    numeric_data = []
    
    with open(filepath, 'r') as file:
        # Read through the file until the marker is found
        for line in file:
            if data_start_marker in line:
                break
        
        # Read the remaining lines and process them
        for line in file:
            parts = line.strip().split(',')
            if len(parts) == 2:
                try:
                    potential = float(parts[0].strip())
                    current = float(parts[1].strip())
                    numeric_data.append([potential, current])
                except ValueError:
                    # This handles lines that cannot be converted to floats
                    continue

    # Convert list to a NumPy array
    return np.array(numeric_data)

def find_cycle_starts(data, start_potential):
    # Finding all indices where the potential approximately matches the start_potential
    potential_close_indices = np.where(np.isclose(data[:,0], start_potential, atol=0.01))[0]
    
    # Determine the actual start of each cycle by checking the difference between consecutive matches
    cycle_starts = [potential_close_indices[0]]
    for i in range(1, len(potential_close_indices)):
        if potential_close_indices[i] - potential_close_indices[i-1] > 10:  # Ensuring a significant index gap
            cycle_starts.append(potential_close_indices[i])
    return cycle_starts

def split_on_cycles(data, fwd_cycle_starts, rev_cycle_starts):
    """Returns fwd cycles and rev cycles in two lists
    """
    fwd_cycles_list = []
    rev_cycles_list = []
    
    #fwd cycles
    for i, start_pt in enumerate(fwd_cycle_starts[:-1]): #last one is half cycle
        cycle_period = slice(start_pt, rev_cycle_starts[i])
        fwd_cycles_list.append(data[cycle_period, :])
    
    #rev cycles
    for i, start_pt in enumerate(rev_cycle_starts):
        cycle_period = slice(start_pt, fwd_cycle_starts[i + 1])
        rev_cycles_list.append(data[cycle_period, :])
    return fwd_cycles_list, rev_cycles_list


def calculate_fwd_cycle_peak_area(
    cycle,
    clip_below_v=0.2,
    nonzero_grad_thresh=1e-6,
    grad_rise_thresh=10e-6
):
    """Locates peak and returns the area. Also renders a plot.
    
    Parameters
    ----------
    cycle: Array (n_samples, 2) of [voltages, currents]
    clip_below: Ignore voltages below this value
    nonzero_grad_thresh: Threshold for clipping the flat part before the peak
    grad_rise_thresh: Used to establish the start of the peak
    
    Returns:
    area under the curve
    """
    voltage, current = [cycle[:, i] for i in range(2)]

    #Gradients
    current_dx = pd.Series(np.gradient(current)).rolling(window=10).mean().values
    abs_current_dx = np.abs(current_dx)
    # current_d2x = pd.Series(np.gradient(np.gradient(current))).rolling(window=10).mean().values

    f, ax = plt.subplots(figsize=(11, 4.5))
    ax.plot(voltage, current, color='black', linewidth=6, label='data')
    ax.set(xlabel='voltage', ylabel='current')

    grad_clr = 'navy'
    ax2 = ax.twinx()
    ax2.plot(voltage, current_dx, color=grad_clr, linewidth=1.5, label='gradient')
    ax2.axhline(0, color=grad_clr, linewidth=1, linestyle=':')
    ax2.set_ylabel('gradient')
    ax2.tick_params(axis='y', colors=grad_clr)
    ax2.yaxis.label.set_color(grad_clr)
    ax2.spines.right.set_color(grad_clr)

    cmap = plt.get_cmap('tab10')
    line_fmt = dict(linestyle='--', linewidth=2, alpha=1)

    #"rel_start" tracks the starting point for each subsequent operation
    ax.axvline(voltage[rel_start := 0], color=cmap(0), **line_fmt, label='cycle starts')

    #Hard clip. Ignore voltages below clip_below_v
    clip_idx = np.argwhere(voltage > clip_below_v)[0].item()
    ax.axvline(voltage[rel_start + clip_idx], color=cmap(1), **line_fmt, label='hard clip ends')
    rel_start += clip_idx

    #Clip the flat part
    flatpart_endidx = np.argwhere(abs_current_dx[rel_start:] > nonzero_grad_thresh)[0].item()
    ax.axvline(voltage[rel_start + flatpart_endidx], color=cmap(2), **line_fmt, label='flat part ends')
    rel_start += flatpart_endidx

    #Peak starts where the gradient begins to rise
    peak_start_idx = np.argwhere(abs_current_dx[rel_start:] > grad_rise_thresh)[0].item()
    ax.axvline(voltage[rel_start + peak_start_idx], color=cmap(3), **line_fmt, label='peak_start_idx')
    peak_x_indices = [rel_start + peak_start_idx]
    rel_start += peak_start_idx

    #Find the midpoint of the downwards slope
    down_slope_idx = np.argmin(current_dx[rel_start:]).item()
    ax.axvline(voltage[rel_start + down_slope_idx], color=cmap(5), **line_fmt, label='down_slope_idx')
    rel_start += down_slope_idx

    #The gradient point closest to 0 after the downwards slope is where the peak ends 
    peak_end_idx = np.argmin(abs_current_dx[rel_start:]).item()
    peak_x_indices.append(rel_start + peak_end_idx)
    ax.axvline(voltage[rel_start + peak_end_idx], color=cmap(6), **line_fmt, label='peak_end_idx')

    f.legend(
        fontsize=8.5, loc='upper left', ncols=10, shadow=True, framealpha=1,
        bbox_to_anchor=(0.05, 0)
    )
    ax.set_xlim(voltage[clip_idx + flatpart_endidx - 50], voltage.max() * 0.9)
    ax.set_ylim(-0.005, 0.0075)
    ax2.set_ylim(ax2.get_ylim()[0] * 0.83, ax2.get_ylim()[1] * 0.8)

    #Points denoting line
    ax.scatter(
        voltage[peak_x_indices], current[peak_x_indices],
        marker='o', color='lime', s=40, zorder=3,
    )

    peak_indices_all = np.arange(*peak_x_indices)
    peak_voltages = voltage[peak_indices_all]
    peak_currents = current[peak_indices_all]

    line_voltages = peak_voltages[[0, -1]]
    line_currents = peak_currents[[0, -1]]
    line_poly = np.polynomial.Polynomial.fit(line_voltages, line_currents, deg=1)
    ax.plot(line_voltages, line_poly(line_voltages), linewidth=3, color='lime')
    ax.plot(peak_voltages, peak_currents, color='lime', linewidth=1.8)
    ax.fill_between(
        x=peak_voltages, y1=line_poly(peak_voltages), y2=peak_currents,
        hatch=None, facecolor='lime', alpha=0.25
    )

    peak_area = np.trapz(peak_currents, peak_voltages)
    triangle_area = 0.5 * np.ptp(line_voltages) * line_currents.max()
    peak_area -= triangle_area #subtract the triangle under the line
    ax.set_title(f'shaded area: {peak_area:>.2e}', weight='bold')
    plt.show()
    
    return peak_area

#Load data
data = load_numeric_data(f'/mnt/c/dev/CV_sample1.txt')

#Define the start potentials of the forward and reverse sweeps
fwd_sweep_start_potential = -1.3
rev_sweep_start_potential = 0.7

#Find the start indices
fwd_cycle_starts = find_cycle_starts(data, fwd_sweep_start_potential)
rev_cycle_starts = find_cycle_starts(data, rev_sweep_start_potential)

#Get the cycles into a list
fwd_cycles_list, rev_cycles_list = split_on_cycles(data, fwd_cycle_starts, rev_cycle_starts)

for i, cycle in enumerate(fwd_cycles_list):
    area = calculate_fwd_cycle_peak_area(cycle)
    print('Peak area of cycle', i, 'is:', area)

#For debugging: view the data splitting results
if False:
    plt.plot(data[:, 0])
    [plt.gca().axvline(x, linewidth=1, linestyle=':', color='tab:red') for x in fwd_cycle_starts]
    [plt.gca().axvline(x, linewidth=1, linestyle=':', color='black') for x in rev_cycle_starts]
    plt.gcf().set_size_inches(8, 2)
    plt.gca().set(
        xlabel='index', ylabel='V',
        title=f'{len(fwd_cycles_list)} fwd cycles/{len(rev_cycles_list)} rev cycles'
    )
    plt.show()

    for cycles_list in [fwd_cycles_list, rev_cycles_list]:
        #View split cycles
        for c in cycles_list:
            plt.plot(c[:, 0], c[:, 1])
        plt.gca().set(xlim=(0.25, 0.65), ylim=(-0.0075, 0.0075))
        plt.gca().set(xlabel='V', ylabel='I')
        plt.gcf().set_size_inches(8, 3)
        plt.show()