如何模拟给定复场方程？现有Python代码报错，需生成干涉拍信号

问题描述 #

需要模拟两个复场方程，叠加后得到干涉拍信号，但编写的Python代码运行报错。

原复场方程 #

$$ E_1 = \sum_{j=0}^{J-1} a_j e^{-i(j\Delta+w_o )t} $$
$$ E_2 = e^{-iw_ot} \sum_{m=J}^{J+M} e^{-im\Delta t}(a_m+b_m e^{-iw_ot}) $$

报错代码 #

%matplotlib inline

import random
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import ifft, fftshift

tstart = -10e-9
tstop = 10e-9
delta = 31.6e6 #rep rate
wo = 5e6    #offset frequency
i = 1j

aj = np.array([1 for i in range(500)]) # no of comblines
t = np.linspace (tstart,tstop, 1000)
E1 = np.linspace(0,0,1000).astype("complex") # PSD

for s in range(len(t)):
    for k in range(len(aj)):
        E1[s]+= aj[k]**2*(np.exp(-i*(k*delta+wo)*t[s]))
        
###########################################################

am = np.array([1 for i in range(500)])
bm = np.array([1 for i in range (500)])# no of comblines
t2 = np.linspace (tstart,tstop, 1000)
E2 = np.linspace(0,0,1000, dtype = "complex_") # PSD
m = np.array([i for i in range (500,1000,1)])

for s2 in range(len(t)):
    for k2 in range(len(m)):
        E2[s2]+= (np.exp(-i*wo*t[s2]))*[np.exp(-i*(m[k2]*delta)*t[s2])*(am[k2]+bm[k2]*np.exp(-i*(wo)*t[s2]))]

# 计划叠加E1和E2

错误分析与修正 #

核心错误点 #

1. 类型不匹配：E2计算时，[np.exp(...)]是列表类型，无法直接与复数相乘，需移除方括号
2. 方程实现偏差：E1中错误使用aj[k]**2，原方程为a_j而非平方项
3. 冗余变量：t2与t完全重复，无需重新定义
4. 效率低下：嵌套循环计算量较大，建议用NumPy向量化运算替代

修正后代码 #

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftshift

# 参数设置
tstart = -10e-9
tstop = 10e-9
delta = 31.6e6  # 重复频率
wo = 5e6        # 偏移频率
i = 1j

# 生成时间轴
t = np.linspace(tstart, tstop, 1000)

# 计算E1：严格匹配原方程，移除多余平方项
J = 500
aj = np.ones(J)
j = np.arange(J)
# 向量化运算：构建所有j对应的项后求和，替代嵌套循环
E1 = np.sum(aj * np.exp(-i * (j * delta + wo) * t[:, np.newaxis]), axis=1)

# 计算E2：修正类型错误，移除冗余变量
M = 500
am = np.ones(M)
bm = np.ones(M)
m = np.arange(J, J+M)
# 先计算求和部分，再乘以e^(-iwo t)
sum_term = np.sum(np.exp(-i * m * delta * t[:, np.newaxis]) * (am + bm * np.exp(-i * wo * t[:, np.newaxis])), axis=1)
E2 = np.exp(-i * wo * t) * sum_term

# 叠加得到总场
E_total = E1 + E2

# 可视化结果：时域信号
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(t * 1e9, np.real(E_total))
plt.xlabel('时间 (ns)')
plt.ylabel('实部振幅')
plt.title('总复场时域信号')

# 频域信号（观察干涉拍信号特征）
freq = fftshift(np.fft.fftfreq(len(t), d=t[1]-t[0]))
E_fft = fftshift(np.abs(fft(E_total)))
plt.subplot(212)
plt.plot(freq / 1e6, E_fft)
plt.xlabel('频率 (MHz)')
plt.ylabel('频谱幅度')
plt.title('总复场频域特征（干涉拍信号）')
plt.tight_layout()
plt.show()

说明 #

修正后的代码通过向量化运算大幅提升计算效率，同时严格匹配原方程定义。叠加后的总场时域信号会呈现干涉拍频特征，频域图中可观察到对应拍频的峰值。

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