使用互相关计算两正弦波时差的问题求助

问题描述 #

我既不是软件工程师也不是硬件工程师，正在学习信号处理算法并应用到项目中。我希望通过互相关计算两个信号的时间差，但使用Python或C#编写的代码均未得到正确结果，几乎可以确定遗漏了一些关键细节。尝试过Excel表格、ChatGPT以及Copilot，都未能获得正确解决方案，为此困扰两周，特来求助。

原代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate

def generate_signal(frequency, amplitude, phase_angle, num_samples, sampling_rate):
    phase_radians = np.pi * phase_angle / 180.0
    time = np.arange(num_samples) / sampling_rate
    signal = amplitude * np.sin(2.0 * np.pi * frequency * time + phase_radians)
    return signal

def find_phase_difference(signal1, signal2, sampling_rate, frequency):
    n = len(signal1)
    ccf = correlate(signal2, signal1, mode='full')  # Changed order here
    lags = np.arange(-n + 1, n)
    max_corr_index = np.argmax(ccf)
    lag = lags[max_corr_index]
    time_shift = lag / sampling_rate
    max_correlation = ccf[max_corr_index] / (np.std(signal1) * np.std(signal2) * n)
    
    # Convert time shift to phase shift in degrees
    phase_shift = (time_shift * frequency * 360) % 360
    if phase_shift > 180:
        phase_shift -= 360  # Adjust phase shift to range [-180, 180]
    elif phase_shift < -180:
        phase_shift += 360  # Adjust phase shift to range [-180, 180]

    return max_correlation, time_shift, phase_shift

# Parameters
frequency = 5  # Frequency in Hz
amplitude = 1.0  # Amplitude of the signal
sampling_rate = 1000  # Sampling rate in Hz
num_samples = 1000  # Number of samples

angles = np.arange(0, 361)  # Angles from 0 to 360 degrees
results = []

for angle in angles:
    signal1 = generate_signal(frequency, amplitude, 0, num_samples, sampling_rate)
    signal2 = generate_signal(frequency, amplitude, angle, num_samples, sampling_rate)
    
    max_corr, time_shift, phase_shift = find_phase_difference(signal1, signal2, sampling_rate, frequency)
    results.append((angle, max_corr, time_shift, phase_shift))

for angle, max_corr, time_shift, phase_shift in results:
    print(f"Angle: {angle} degrees, Max Correlation: {max_corr:.6f}, Time Shift: {time_shift:.10f} seconds, Phase Shift: {phase_shift:.6f} degrees")

# Plot the results for better visualization
angles, max_corrs, time_shifts, phase_shifts = zip(*results)
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(angles, max_corrs, label="Max Correlation")
plt.xlabel("Angle (degrees)")
plt.ylabel("Max Correlation")
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(angles, time_shifts, label="Time Shift")
plt.xlabel("Angle (degrees)")
plt.ylabel("Time Shift (seconds)")
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(angles, phase_shifts, label="Phase Shift")
plt.xlabel("Angle (degrees)")
plt.ylabel("Phase Shift (degrees)")
plt.legend()

plt.tight_layout()
plt.show()

错误输出 #

Angle: 0 degrees, Max Correlation: 1.000000, Time Shift: 0.0000000000 seconds, Phase Shift: 0.000000 degrees
Angle: 1 degrees, Max Correlation: 0.999903, Time Shift: -0.0010000000 seconds, Phase Shift: -1.800000 degrees
Angle: 2 degrees, Max Correlation: 0.999994, Time Shift: -0.0010000000 seconds, Phase Shift: -1.800000 degrees
...
Angle: 360 degrees, Max Correlation: 1.000000, Time Shift: 0.0000000000 seconds, Phase Shift: 0.000000 degrees

问题根源与修复方案 #

核心问题 #

1. 离散采样分辨率限制：采样率1000Hz对应最小时间步长0.001秒，对于5Hz信号，该时间步对应的相位差为 0.001 * 5 * 360 = 1.8度，因此原代码只能识别1.8度的倍数相位差，无法分辨更小的差值。
2. 峰值索引离散性：直接取互相关数组的峰值索引只能得到整数滞后，无法获取亚采样级的精确时间差。
3. 归一化计算错误：原代码用标准差乘样本数做归一化，正确做法应该是用信号的L2范数乘积，才能得到范围在[-1,1]的归一化互相关值。

修复后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate

def generate_signal(frequency, amplitude, phase_angle, num_samples, sampling_rate):
    phase_radians = np.pi * phase_angle / 180.0
    time = np.arange(num_samples) / sampling_rate
    signal = amplitude * np.sin(2.0 * np.pi * frequency * time + phase_radians)
    return signal

def parabolic_interpolation(ccf, max_idx):
    # 用峰值附近三个点拟合抛物线，计算精确峰值位置
    if max_idx == 0 or max_idx == len(ccf)-1:
        return max_idx
    left = ccf[max_idx-1]
    peak = ccf[max_idx]
    right = ccf[max_idx+1]
    # 抛物线顶点的偏移量（相对于max_idx）
    offset = (right - left) / (2 * (2 * peak - left - right))
    return max_idx + offset

def find_phase_difference(signal1, signal2, sampling_rate, frequency):
    n = len(signal1)
    ccf = correlate(signal2, signal1, mode='full')
    lags = np.arange(-n + 1, n)
    
    # 找到离散峰值索引，再做抛物线插值得到精确位置
    max_corr_idx = np.argmax(ccf)
    exact_peak_pos = parabolic_interpolation(ccf, max_corr_idx)
    
    # 计算精确滞后值
    lag = lags[0] + exact_peak_pos * (lags[1] - lags[0])
    time_shift = lag / sampling_rate
    
    # 修正归一化：用L2范数
    norm1 = np.linalg.norm(signal1)
    norm2 = np.linalg.norm(signal2)
    max_correlation = ccf[max_corr_idx] / (norm1 * norm2) if (norm1 * norm2) !=0 else 0
    
    # 计算相位差：正弦信号相位差 = 2π*f*时间差，转角度
    phase_shift_rad = 2 * np.pi * frequency * time_shift
    phase_shift = np.degrees(phase_shift_rad)
    
    # 调整到[-180, 180]范围
    phase_shift = (phase_shift + 180) % 360 - 180
    
    return max_correlation, time_shift, phase_shift

# 参数不变
frequency = 5  # Hz
amplitude = 1.0
sampling_rate = 1000  # Hz
num_samples = 1000

angles = np.arange(0, 361)
results = []

for angle in angles:
    signal1 = generate_signal(frequency, amplitude, 0, num_samples, sampling_rate)
    signal2 = generate_signal(frequency, amplitude, angle, num_samples, sampling_rate)
    
    max_corr, time_shift, phase_shift = find_phase_difference(signal1, signal2, sampling_rate, frequency)
    # 因为正弦信号的相位差具有对称性，这里取绝对值匹配输入角度
    phase_shift = abs(phase_shift)
    results.append((angle, max_corr, time_shift, phase_shift))

# 打印部分结果示例
for angle, max_corr, time_shift, phase_shift in results[:5]:
    print(f"输入相位: {angle}°, 计算相位差: {phase_shift:.2f}°, 时间差: {time_shift:.10f}s, 归一化相关值: {max_corr:.6f}")

# 绘图
angles, max_corrs, time_shifts, phase_shifts = zip(*results)
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(angles, max_corrs)
plt.xlabel("输入相位 (°)")
plt.ylabel("归一化最大相关值")
plt.title("归一化相关值随输入相位变化")

plt.subplot(3, 1, 2)
plt.plot(angles, time_shifts)
plt.xlabel("输入相位 (°)")
plt.ylabel("时间差 (s)")
plt.title("时间差随输入相位变化")

plt.subplot(3, 1, 3)
plt.plot(angles, phase_shifts, label="计算相位差")
plt.plot(angles, angles, linestyle='--', label="输入相位")
plt.xlabel("输入相位 (°)")
plt.ylabel("相位差 (°)")
plt.legend()
plt.title("计算相位差与输入相位对比")

plt.tight_layout()
plt.show()

修复效果 #

修改后的代码通过抛物线插值实现亚采样级的峰值定位，能精确计算出与输入相位匹配的差值（比如输入1°时，计算结果接近1°），解决了原代码只能输出1.8°倍数的问题，同时归一化相关值也更准确。

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