对称自相关函数实现及Python代码正确性评估请求

需求说明 #

我需要计算两个复值自相关核函数：
$$R_z[n,t] = z[n+t]z^[n-t]$$
$$R_z[n,t'] = z[n+t']z^[n-t']$$
其中：

• $*$ 表示复共轭
• $z[n] = \exp\left(1i \cdot 2\pi \cdot \frac{\text{cumsum}(wav_file_signal[n])}{wav_fs}\right)$
• $t' = t + \text{int}(wav_fs/4)$
• $N$ 为信号长度，$t$ 在指定范围内取值

待评估代码 #

import numpy as np

# Example usage
fs = 100  # Sample rate
Ts = 1 / fs  # Time step
N = 1000  # Length of the signal
t = np.arange(0, N) * Ts  # Time vector
f0 = 5  # Fundamental frequency of the signal

# Generate a noisy signal for testing
z = np.sin(2 * np.pi * f0 * t) + 0.5 * np.random.randn(N)

tau_lim = 30
taus = np.arange(-tau_lim, tau_lim + 1)
alphas = round(fs / (4 * np.max(np.abs(z))))
tau_delta = round(fs / (4 * f0))

Rx = np.zeros((len(z), len(taus)))
R_1_result = np.zeros(len(z))
R_2_result = np.zeros(len(z))

z = np.exp(1j * 2 * np.pi * alphas * np.cumsum(z) / fs)


# Circular shift function
def circular_shift(signal, shift_amount):
    return np.roll(signal, shift_amount)


for i, tau in enumerate(taus):
    # R_1_backward_shift = z[:len(z) - abs(tau)]
    # R_1_forward_shift = z[abs(tau):]

    R_1_backward_shift = circular_shift(z, - abs(tau))
    R_1_forward_shift = circular_shift(z, abs(tau))

    R_1_length = len(R_1_forward_shift)  # R_1 and R_2 have the same length

    R_1_result = R_1_forward_shift * np.conj(R_1_backward_shift)

    # R_2_backward_shift = z[:len(z) - abs(tau) - tau_delta]
    # R_2_forward_shift = z[abs(tau) + tau_delta:]

    R_2_backward_shift = circular_shift(z, - abs(tau) - tau_delta)
    R_2_forward_shift = circular_shift(z, abs(tau) + tau_delta)

    R_2_length = len(R_2_forward_shift)

    R_2_result = R_2_forward_shift * np.conj(R_2_backward_shift)

代码问题分析 #

1. 结果存储逻辑完全错误 #

循环中每次迭代都会覆盖R_1_result和R_2_result，最终仅保留最后一个tau值的计算结果，完全丢失了其他t（即tau）对应的结果，不符合“t在指定范围内取值”的需求。

2. 边界处理错误：循环移位不符合公式定义 #

公式中的$z[n+t]$和$z[n-t]$要求：当$n+t$或$n-t$超出信号索引范围（0到N-1）时，该位置的$R_z$值无效，应舍弃或设为0。但代码用np.roll实现循环移位，会将信号首尾样本拼接复用，完全违背了原公式的边界约束。

3. t'的偏移量计算错误 #

公式中$t' = t + \text{int}(fs/4)$，但代码中tau_delta = round(fs/(4*f0))引入了无关的基频$f0$，属于公式映射的明显错误，应改为tau_delta = int(fs/4)。

4. 测试阶段变量混淆 #

测试信号生成时，最初的z是实值正弦加噪声，后续重新赋值为复值信号逻辑没问题，但alphas的计算基于实值z的最大值，和原公式中基于wav信号累积和的定义不符，属于测试代码的瑕疵。

修正示例片段 #

# 初始化二维数组存储所有tau对应的结果
R1 = np.zeros((N, len(taus)), dtype=np.complex128)
R2 = np.zeros((N, len(taus)), dtype=np.complex128)
tau_delta = int(fs / 4)  # 修正t'的偏移量计算

for i, tau in enumerate(taus):
    abs_tau = abs(tau)
    # 计算R_z[n,tau]，仅保留有效索引范围
    valid_n = np.arange(abs_tau, N - abs_tau)
    R1[valid_n, i] = z[valid_n + abs_tau] * np.conj(z[valid_n - abs_tau])
    
    # 计算R_z[n,tau']，tau'=tau+tau_delta
    tau_prime = tau + tau_delta
    abs_tau_prime = abs(tau_prime)
    valid_n_prime = np.arange(abs_tau_prime, N - abs_tau_prime)
    R2[valid_n_prime, i] = z[valid_n_prime + abs_tau_prime] * np.conj(z[valid_n_prime - abs_tau_prime])

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