复杂函数傅里叶变换代码求助：功率谱不随参数fb变化

我尝试对特定形式的激光场自相关函数做傅里叶变换，预期得到含双边带的功率谱。参考相关论文编写Python代码后，发现调整参数fb时，功率谱完全没有变化，怀疑是傅里叶变换的实现有问题，附上代码求助：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sici
from scipy.fft import fft, ifft, fftfreq, fftshift

def plot_autocorrelation_and_spectrum(x, p, ha, hb, fb, w0, wb):
    a0 = np.exp(w0 * x * 1j)
    a1 = np.exp(-(ha - hb / fb) * (wb * x * sici(wb * x)[0] - 2 * np.sin(wb * x / 2)**2))
    a2 = np.exp(-hb * np.pi**2 * np.abs(x))
    
    autocorr = (p * a0 * a1 * a2) / np.max(a0 * a1 * a2)
    
    t = np.linspace(min(x), max(x), len(autocorr))
    
    plt.subplot(3, 1, 1)
    plt.plot(t, autocorr)
    plt.title('Autocorrelation Plot')
    plt.xlabel('Time')
    plt.ylabel('Autocorrelation')
    
    fft_result = np.fft.fft(autocorr)
    power_spectrum = np.abs(fft_result) ** 2
    power_spectrum = np.fft.fftshift(power_spectrum)/max(power_spectrum)
    
    frequency = np.fft.fftfreq(len(t), d=t[1] - t[0])
    frequency = np.fft.fftshift(frequency)
    
    plt.subplot(3, 1, 2)
    plt.plot(frequency, power_spectrum)
    plt.title('Original Power Spectrum')
    plt.xlabel('Frequency')
    plt.ylabel('Power Spectrum')

x_values = np.linspace(0, 10, 1000) 
plot_autocorrelation_and_spectrum(x_values, p=1.0, ha=100.0, hb=1000.0, fb=100.0, w0=1, wb=2 * np.pi * fb)

问题根源与修正方案 #

核心问题

1. 时间轴非对称：自相关函数本质是关于$\tau=0$对称的，但当前代码仅取了$x>0$的范围，丢失负时间部分信息，导致傅里叶变换无法生成双边带特征，且参数变化的影响被掩盖。
2. 复数归一化错误：直接对复数数组取max会默认比较实部，丢失虚部信息，导致归一化失真。
3. 冗余时间变量：重新生成的t数组与输入x完全一致，属于冗余操作，且未解决时间轴不对称问题。

修正后的代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sici
from scipy.fft import fft, fftfreq, fftshift

def plot_autocorrelation_and_spectrum(x, p, ha, hb, fb, w0, wb):
    # 计算自相关函数各分量
    a0 = np.exp(w0 * x * 1j)
    sici_val = sici(wb * x)[0]
    a1 = np.exp(-(ha - hb / fb) * (wb * x * sici_val - 2 * np.sin(wb * x / 2)**2))
    a2 = np.exp(-hb * np.pi**2 * np.abs(x))
    
    autocorr = p * a0 * a1 * a2
    # 基于复数模的最大值归一化，避免实部比较的错误
    autocorr = autocorr / np.max(np.abs(autocorr))
    
    # 绘制自相关函数的实部与虚部
    plt.subplot(3, 1, 1)
    plt.plot(x, np.real(autocorr), label='实部')
    plt.plot(x, np.imag(autocorr), label='虚部')
    plt.title('自相关函数')
    plt.xlabel('时间 τ')
    plt.ylabel('自相关值')
    plt.legend()
    
    # 傅里叶变换与功率谱计算
    fft_result = fft(autocorr)
    power_spectrum = np.abs(fft_result) ** 2
    power_spectrum = fftshift(power_spectrum) / np.max(power_spectrum)
    
    frequency = fftfreq(len(x), d=x[1] - x[0])
    frequency = fftshift(frequency)
    
    plt.subplot(3, 1, 2)
    plt.plot(frequency, power_spectrum)
    plt.title('功率谱')
    plt.xlabel('频率')
    plt.ylabel('归一化功率')
    plt.xlim(-500, 500)  # 根据参数范围调整显示区间

# 生成对称时间轴，包含正负时间
x_values = np.linspace(-10, 10, 2000) 
# 测试不同fb值，观察功率谱变化
for fb in [50, 100, 150]:
    plt.figure(figsize=(10,8))
    plot_autocorrelation_and_spectrum(x_values, p=1.0, ha=100.0, hb=1000.0, fb=fb, w0=1, wb=2 * np.pi * fb)
    plt.suptitle(f'参数fb = {fb} 时的结果')
    plt.tight_layout()
    plt.show()

关键修正点

• 对称时间轴：生成从-10到10的对称数组，还原自相关函数的对称性，确保傅里叶变换能正确生成双边带。
• 正确归一化：使用np.max(np.abs(autocorr))对复数自相关函数归一化，保留完整的复振幅信息。
• 可视化优化：绘制自相关函数的实部和虚部，便于观察参数变化的影响；添加频率范围限制，突出双边带特征。

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