我正在尝试从我正在阅读的实验数据集中制作功率谱，然后将其拟合到理论曲线。现在一切正常，我没有出错，除了我的曲线与数据相差 1000 倍，我完全不知道问题可能是什么。我问了几个人，但无济于事。（希望大家能帮忙）
无论如何，我很确定它不是单位，因为我和其他 2 人对它们进行了三次检查。基本上，我需要使用最小二乘法将功率谱拟合到方程中。我不能发布整个代码，因为它相当长而且有点混乱，但这是傅立叶部分，我为所有未在代码中声明的数组和变量添加了注释）
#Calculate stuff
Nm = 10**-6 #micro to meter
KbT = 4.10E-21 #Joule
T = 297. #K
l = zvalue*Nm #meter
meany = np.mean(cleandatay*Nm) #meter (cleandata is the array that I read in from a cvs at the start.)
SDy = sum((cleandatay*Nm - meany)**2)/len(cleandatay) #meter^2
FmArray[0][i] = ((KbT*l)/SDy) #N
#print FmArray[0][i]
print float((i*100/len(filelist)))#how many % done?
#fourier
dt = cleant[1]-cleant[0] #timestep
N = len(cleandatay) #Same for cleant, its the corresponding time to cleandatay
这是傅立叶部分开始的地方，我将 fft 转换为功率谱。然后我用数组freqs计算相应的频率步长
fouriery = np.fft.fft((cleandatay*(10**-6)))
fourierpower = (np.abs(fouriery))**2
fourierpower = fourierpower[1:N/2] #remove 0th datapoint and /2 (remove negative freqs)
fourierpower = fourierpower*dt #*dt to account for steps
freqs = (1.+np.arange((N/2)-1.))/50.
#Least squares method
eta = 8.9E-4 #pa*s
Rbead = 0.5E-6#meter
constant = 2*KbT/(3*eta*pi*Rbead)
omega = 2*pi*freqs #rad/s
Wcarray = 2.*pi*np.arange(0,30, 0.02003) #0.02 = 30/len(freqs)
ChiSq = np.zeros(len(Wcarray))
for k in range(0, len(Wcarray)):
Py = (constant / (Wcarray[k]**2 + omega**2))
ChiSq[k] = sum((fourierpower - Py)**2)
pylab.loglog(omega, Py)
print k*100/len(Wcarray)
index = np.where(ChiSq == min(ChiSq))
cutoffw = Wcarray[index]
Pygoed = (constant / (Wcarray[index]**2 + omega**2))
print cutoffw
print constant
print min(ChiSq)
pylab.loglog(omega,ChiSq)
所以我不知道可能出了什么问题，我认为它是 fft，因为没有其他东西真的会出错。下面是我根据光谱绘制所有拟合线时得到的图片，你可以看到它偏离了大约 1000（实际上正好是 1000，因为这留下了 10^-22 的最小二乘残差，但我不能只是随机相乘而不知道为什么）
只是为了详细说明图片。绿点是 fft 频谱，线条是拟合，红点是它认为的截止频率，蓝线是卡方拟合，寻找最低值。