基于双段基线线性回归的数据集目标曲线下面积计算求助

我需要计算电位-电流曲线的特定区域面积，但曲线起始和终点未归零。具体需求如下：

• 分别选取预基线和后基线的数据点进行线性回归
• 计算原始曲线与这两条基线所围成区域的面积

原始数据曲线：

目标面积示意图：

更具体的操作逻辑：选择预基线、后基线的数据点分别做线性回归，基于预基线回归的最后一个点与后基线回归的第一个点，计算曲线与两条基线间限定区域的面积。

我尝试用ChatGPT生成实现代码，但代码未正确理解需求——它计算的并非我需要的双基线与曲线之间的面积，现有代码如下：

from scipy.integrate import trapz
from scipy.stats import linregress

def baseline_correction(data, indices):
    # Calculate a unified baseline using linear regression on selected indices
    slope, intercept, _, _, _ = linregress(data[indices, 0], data[indices, 1])
    return slope, intercept

def calculate_area(data, start_index, end_index, slope, intercept):
    # Correct the current by subtracting the baseline
    corrected_current = data[start_index:end_index, 1] - (data[start_index:end_index, 0] * slope + intercept)
    # Calculate the area (charge) under the corrected current curve
    charge = trapz(corrected_current, x=data[start_index:end_index, 0])
    return charge, corrected_current

def process_cycles(data, cycle_starts):
    charges = []
    for i in range(len(cycle_starts) - 1):
        start_index = cycle_starts[i]
        end_index = cycle_starts[i + 1]

        # Selecting only the increasing potential part of the cycle
        increasing_potential_indices = np.where(np.diff(data[start_index:end_index, 0]) > 0)[0] + start_index
        
        # Combine pre-peak and post-peak data points for baseline
        pre_peak_indices = increasing_potential_indices[np.where((data[increasing_potential_indices, 0] >= 0.35) & (data[increasing_potential_indices, 0] <= 0.4))]
        post_peak_indices = increasing_potential_indices[np.where((data[increasing_potential_indices, 0] >= 0.525) & (data[increasing_potential_indices, 0] <= 0.55))]
        combined_indices = np.concatenate((pre_peak_indices, post_peak_indices))
        
        # Calculate unified baseline
        slope, intercept = baseline_correction(data, combined_indices)
        
        # Define peak region (between 0.35V to 0.4V)
        peak_start_idx = np.min(np.where(data[:, 0] >= 0.4)[0])
        peak_end_idx = np.max(np.where(data[:, 0] <= 0.525)[0])
        
        # Calculate the charge
        charge, corrected_current = calculate_area(data, peak_start_idx, peak_end_idx, slope, intercept)
        charges.append(charge)
        
        # Visualization
        plt.figure(figsize=(10, 6))
        plt.plot(data[start_index:end_index, 0], data[start_index:end_index, 1], label='Original Data', color='gray')
        plt.plot(data[combined_indices, 0], data[combined_indices, 1], 'o', label='Baseline Points', color='red')
        baseline_curve = data[peak_start_idx:peak_end_idx, 0] * slope + intercept
        plt.plot(data[peak_start_idx:peak_end_idx, 0], baseline_curve, label='Baseline', color='green')
        plt.plot(data[peak_start_idx:peak_end_idx, 0], corrected_current, label='Corrected Current', color='purple')
        plt.fill_between(data[peak_start_idx:peak_end_idx, 0], 0, corrected_current, color='purple', alpha=0.3, label='Area Under Curve')
        plt.xlabel('Potential (V)')
        plt.ylabel('Current (A)')
        plt.title(f'Cycle {i+1} Charge Calculation')
        plt.legend()
        plt.grid(True)
        plt.show()
    
    return charges

补充信息 #

数据加载与预处理代码：

import numpy as np
import matplotlib.pyplot as plt

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 plot_cycle_data(data, cycle_starts, cycle_number):
    if cycle_number > len(cycle_starts):
        print("Cycle number exceeds the available cycles.")
        return
    
    start_index = cycle_starts[cycle_number - 1]
    end_index = cycle_starts[cycle_number] if cycle_number < len(cycle_starts) else len(data)
    
    plt.figure(figsize=(10, 5))
    plt.plot(data[start_index:end_index, 0], data[start_index:end_index, 1], label=f'Cycle {cycle_number}')
    plt.xlabel('Potential (V)')
    plt.ylabel('Current (A)')
    plt.title(f'Cycle {cycle_number} from Potential {data[start_index, 0]}V to {data[end_index-1, 0]}V')
    plt.legend()
    plt.grid(True)
    plt.show()

# Usage
file_path = 'path_to_your_file.txt'
data = load_numeric_data(file_path)
cycle_starts = find_cycle_starts(data, -1.3)
plot_cycle_data(data, cycle_starts, 1)  # Plot the first cycle
charges = process_cycles(data, cycle_starts)
print(charges)

内容的提问来源于stack exchange，提问作者Mortenfre96