随机共振（Stochastic Resonance）下含噪信号与周期驱动力的相移计算问题

问题描述 #

计算含噪信号与周期驱动力之间的相移时，互相关方法得到的结果极小且不符合预期，尝试滤波也未解决问题，求准确计算相移的可行方法。

原仿真代码及结果：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.signal as signal

simulation_time = 100.0
sigma = 0.3
A = 10
omega = 0.1
dt = 0.001
sqdt = np.sqrt(dt) #Precompute
time = np.arange(0, simulation_time, dt)
x = np.empty(len(time))
x2 = np.empty(len(time))
#Initial conditions 
x[0] = 1.0

i = 1
for t in time[1:]:
    fx = x[i-1] - x[i-1]*x[i-1]*x[i-1] + A*np.cos(omega*t)
    x2[i] = A*np.cos(omega*t)
    gx = sigma
    x[i] = x[i-1] + dt*fx + sqdt*gx*np.random.standard_normal()
    i += 1

plt.figure()
plt.plot(time, x)
plt.plot(time, x2)
plt.xlabel('time')
plt.ylabel('signal')
plt.show()

#cross corelation method 
output = np.correlate(x,x2*0.15,mode='same')
lags   = time-50
plt.plot(lags,output)
maximum = np.max(output)
phase_shift = lags[output==maximum][0]
phase_shift

# 输出结果：0.05799999999999983

原方法问题分析 #

1. 你的系统是非线性系统（更新方程含x³项），输出信号除驱动力基频外还包含谐波成分，直接互相关会受谐波和噪声干扰，无法准确捕捉基频相位差。
2. 原互相关处理存在错误：x2*0.15的无依据缩放、mode='same'的滞后计算逻辑不严谨，导致结果偏差。

可行解决方案 #

方法1：频率域相位提取（FFT） #

利用FFT提取驱动力基频对应的相位，直接计算两者的相位差，不受谐波和噪声的干扰（只要基频成分足够显著）。

# FFT方法计算相移
n = len(time)
freq = np.fft.fftfreq(n, dt)
# 转换驱动力角频率为Hz，找到对应频谱索引
target_freq = omega/(2*np.pi)
target_idx = np.argmin(np.abs(freq - target_freq))

# 计算FFT并提取相位
fft_x = np.fft.fft(x)
fft_x2 = np.fft.fft(x2)
phase_x = np.angle(fft_x[target_idx])
phase_x2 = np.angle(fft_x2[target_idx])

# 计算相位差并转换为时间相移，调整到[-T/2, T/2]范围
T_drive = 2*np.pi/omega
phase_diff = phase_x - phase_x2
time_shift = phase_diff / omega
if time_shift > T_drive/2:
    time_shift -= T_drive
elif time_shift < -T_drive/2:
    time_shift += T_drive

print(f"FFT方法得到的时间相移：{time_shift}")

方法2：希尔伯特变换求瞬时相位 #

通过希尔伯特变换得到解析信号，提取瞬时相位后与驱动力相位做平均，抵消噪声和暂态的影响。

# 希尔伯特变换方法
from scipy.signal import hilbert

# 计算解析信号并提取解缠绕后的瞬时相位
analytic_x = hilbert(x)
inst_phase_x = np.unwrap(np.angle(analytic_x))

# 驱动力的瞬时相位
inst_phase_x2 = omega * time

# 去掉初始暂态（前10%数据），计算平均相位差
transient_idx = int(0.1 * len(time))
phase_diff_avg = np.mean(inst_phase_x[transient_idx:] - inst_phase_x2[transient_idx:])

# 转换为时间相移并调整范围
time_shift_hilbert = phase_diff_avg / omega
if time_shift_hilbert > T_drive/2:
    time_shift_hilbert -= T_drive
elif time_shift_hilbert < -T_drive/2:
    time_shift_hilbert += T_drive

print(f"希尔伯特变换方法得到的时间相移：{time_shift_hilbert}")

方法3：改进的互相关方法 #

先通过带通滤波提取信号的基频成分，再做互相关计算，避免谐波干扰。

# 改进的互相关方法
fs = 1/dt
center_freq = omega/(2*np.pi)
bandwidth = 0.01  # 窄带宽滤波，只保留基频

# 设计零相位带通滤波器
b, a = signal.butter(4, [center_freq - bandwidth/2, center_freq + bandwidth/2], 
                     btype='bandpass', fs=fs)
x_filtered = signal.filtfilt(b, a, x)

# 计算全范围互相关并获取对应滞后
corr = np.correlate(x_filtered, x2, mode='full')
lags_corr = np.arange(-len(time)+1, len(time)) * dt

# 找到互相关最大值对应的滞后，调整到周期范围内
max_idx = np.argmax(corr)
time_shift_corr = lags_corr[max_idx]
if time_shift_corr > T_drive/2:
    time_shift_corr -= T_drive
elif time_shift_corr < -T_drive/2:
    time_shift_corr += T_drive

print(f"改进互相关方法得到的时间相移：{time_shift_corr}")

关键注意事项 #

• 非线性系统输出存在谐波，必须先提取基频成分再计算相移，否则结果会受谐波干扰。
• 噪声的影响可通过滤波、时间平均（去掉暂态）等方式降低。
• 最终相移需调整到驱动力周期的[-T/2, T/2]范围内，避免出现超过半周期的不合理结果。

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