大规模数据下Scipy信号相关性直接法：高效计算局部滞后范围

import numpy as np

def custom_lagged_correlation(x1, x2, max_lag):
    """
    计算两个一维数组在指定滞后范围内的互相关。

    该函数通过直接法迭代计算从 -max_lag 到 +max_lag 的所有滞后，
    适用于输入数组非常大，但仅需局部滞后结果的场景，
    尤其当 FFT 方法因数据特性（如稀疏性）不适用时。

    参数:
    x1 (array_like): 第一个输入数组（或可转换为 NumPy 数组的对象）。
    x2 (array_like): 第二个输入数组（或可转换为 NumPy 数组的对象）。
    max_lag (int): 正整数，定义计算的滞后范围为 [-max_lag, +max_lag]。

    返回:
    np.ndarray: 包含对应滞后值的互相关结果。
                结果数组的索引 `max_lag + i` 对应滞后 `i`。
                例如，`results[max_lag]` 对应滞后 0。
    """
    # 确保输入是NumPy数组，避免在切片操作时产生不必要的副本，
    # 尤其对于大规模数据，内存效率至关重要。
    x1 = np.asarray(x1)
    x2 = np.asarray(x2)

    # 初始化结果数组。其长度为 2*max_lag + 1，
    # 用于存储从 -max_lag 到 +max_lag 的所有滞后结果。
    # 索引 `max_lag` 处存储滞后为 0 的结果。
    correlation_results = np.zeros(2 * max_lag + 1)

    # 遍历从 -max_lag 到 +max_lag 的每一个滞后值
    for lag_i in range(-max_lag, max_lag + 1):
        # 根据当前的滞后值 `lag_i` 获取 `x1` 和 `x2` 的重叠切片
        if lag_i < 0:
            # 如果滞后为负（x2相对于x1向左移动），x1 从头开始，x2 从 `-lag_i` 处开始
            slice1, slice2 = x1, x2[-lag_i:]
        else:
            # 如果滞后为正（x1相对于x2向左移动），x1 从 `lag_i` 处开始，x2 从头开始
            slice1, slice2 = x1[lag_i:], x2

        # 裁剪两个切片，使其长度相同。
        # 这是为了确保点积操作在有效的重叠区域进行。
        common_length = min(len(slice1), len(slice2))
        slice1 = slice1[:common_length]
        slice2 = slice2[:common_length]

        # 计算重叠部分的点积。
        # 点积操作 `np.dot(slice1, slice2)` 实际上计算了这两个切片的互相关值。
        # 将结果存储到 `correlation_results` 数组中，
        # 索引 `max_lag + lag_i` 确保了结果与实际滞后值的一一对应。
        correlation_results[max_lag + lag_i] = np.dot(slice1, slice2)

    return correlation_results

# 示例用法
# 模拟实际的大规模数据场景，但为演示方便，使用较小的数据长度
data_length = 10000 # 假设实际数据长度可能为 240_000_000
base_noise = np.random.randn(data_length) * 0.1 # 模拟背景噪声

# 创建一个小的信号模式
signal_pattern = np.sin(np.linspace(0, 4 * np.pi, 20)) # 一个正弦波模式

# 将信号模式嵌入到两个数组中，并引入一个已知滞后
x1_large = base_noise.copy()
x2_large = base_noise.copy()

# 在 x1 中嵌入信号
x1_start_idx = data_length // 2
x1_large[x1_start_idx : x1_start_idx + len(signal_pattern)] += signal_pattern

# 在 x2 中嵌入相同信号，但滞后 5 个单位
known_lag = 5
x2_start_idx = x1_start_idx + known_lag
x2_large[x2_start_idx : x2_start_idx + len(signal_pattern)] += signal_pattern

# 定义我们感兴趣的最大滞后范围
max_lag_to_compute = 20 # 我们只关心从 -20 到 +20 的滞后

print(f"模拟数据长度: {data_length}")
print(f"将计算的滞后范围: [{-max_lag_to_compute}, {max_lag_to_compute}]")

# 调用自定义函数计算局部滞后互相关
results = custom_lagged_correlation(x1_large, x2_large, max_lag_to_compute)

# 分析结果
# 找到最大相关性的滞后索引
peak_lag_index = np.argmax(results)
# 将索引转换为实际的滞后值
actual_peak_lag = peak_lag_index - max_lag_to_compute

print(f"互相关结果数组长度: {len(results)}")
print(f"最大相关性值: {results[peak_lag_index]:.4f}")
print(f"最大相关性发生在滞后: {actual_peak_lag}")

# 预期结果：actual_peak_lag 应该接近 known_lag (即 5)