読み込んでいる実験データセットからパワースペクトルを作成し、それを理論曲線に当てはめようとしています。現在、すべてが正常に機能しており、エラーは発生していませんが、曲線がデータとは 1000 倍も異なっており、何が問題なのかまったくわかりません。数人に聞いてみましたがだめでした。（皆様のお役に立てれば幸いです）
とにかく、ユニットは私と他の2人によって3回チェックされたので、ユニットではないことは確かです. 基本的に、最小二乗法を使用して、powerspectrum を方程式に適合させる必要があります。かなり長くて少し面倒なので、コード全体を投稿することはできませんが、これはフーリエ部分です。コードで宣言されていないすべての配列と変数にコメントを追加しました)
#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 を使用して、対応する freq ステップを計算します。
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 スペクトル、線は適合、赤色の点はカットオフ周波数であると考えられる場所、青色の線は最小値を探すカイ 2 乗適合です。