【发布时间】：2012-12-21 12:05:41
【问题描述】：
我正在尝试从我正在读取的实验数据集中制作功率谱，然后将其拟合到理论曲线。现在一切正常，我没有收到错误，除了我的曲线与数据的差异一直为 1000 倍，我完全不知道问题可能是什么。我问了几个人，但无济于事。 （希望大家多多帮忙）
无论如何，我很确定它不是单位，因为我和其他 2 人对它们进行了三次检查。基本上，我需要使用最小二乘法将功率谱拟合到方程中。 我不能发布整个代码，因为它相当长而且有点凌乱，但这是傅立叶部分，我将 cmets 添加到所有未在代码中声明的数组和变量中）
#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 频谱，线条是拟合，红点是它认为截止频率所在的位置，蓝线是卡方拟合，寻找最低值。
【问题讨论】：
-
尝试对您在 numpy 中生成的已知信号进行分析，您可以准确地知道功率谱应该是什么样子，例如，一个简单的正弦波，然后看看结果如何。 （此外，这对于本网站上的人来说更容易调试，因为我们可以运行代码，假设有人如此倾向于。）这些示例是数据分析的单元测试，并且会给你对结果的合理信心。