2017-02-13 3 views
2

私はいくつかのデータの導関数を計算しようとしていますが、私は有限差分とスペクトル法出力の出力を比較しようとしていました。しかし、結果は非常に異なり、私は正確に理由を把握することはできません。`numpy.diff`と` scipy.fftpack.diff`は微分時に異なる結果を出します

import numpy as np 
from scipy import fftpack as sp 
from matplotlib import pyplot as plt 
x = np.arange(-100,100,1) 
y = np.sin(x) 

plt.plot(np.diff(y)/np.diff(x)) 
plt.plot(sp.diff(y)) 

plt.show() 

以下のサンプルコードは、これは、以下の結果 enter image description here

オレンジ出力fftpack出力される出力することを検討してください。微妙なことを心配しないでください。これは単なる例のためです。

なぜ、それらはどう違うのですか?彼らは(ほぼ)同じではないでしょうか?

fftpack.diffの周期キーワードでさまざまな振幅を補正できることは間違いありませんが、どの周期が正しいかはわかりません(私はperiod=1だと思っていましたが動作しません)。

また、numpyを使用して独自のスペクトル差別化を行うにはどうすればよいですか?

答えて

4

関数scipy.fftpack.diffは、導関数を計算しますが、入力は周期的であるとみなします。 period引数は、入力シーケンスの期間(すなわち、x間隔の合計長さ)を与えます。

あなたの場合、これはlen(x)*dxです。dx = x[1] - x[0]です。

単純な(中央の)有限差分(青色)と引数(赤色)を使用したdiffの結果をプロットするコードを示します。変数xyは、あなたのコード内で使用したものと同じである:あなたの入力は、実際に定期的でない場合は、diffによって計算された誘導体は、の端部付近が不正確になること

In [115]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b') 
Out[115]: [<matplotlib.lines.Line2D at 0x1188d01d0>] 

In [116]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r') 
Out[116]: [<matplotlib.lines.Line2D at 0x1188fc9d0>] 

In [117]: plt.xlabel('x') 
Out[117]: <matplotlib.text.Text at 0x1157425d0> 

plot

注意間隔。ここで

は、区間内の正弦関数のひとつの完全な期間が含まれているより短い配列を使用して、もう一つの例だ[0、1]:

In [149]: x = np.linspace(0, 1, 20, endpoint=False) 

In [150]: y = np.sin(2*np.pi*x) 

In [151]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b') 
Out[151]: [<matplotlib.lines.Line2D at 0x119872d90>] 

In [152]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r') 
Out[152]: [<matplotlib.lines.Line2D at 0x119c49090>] 

In [153]: plt.xlabel('x') 
Out[153]: <matplotlib.text.Text at 0x1197823d0> 

plot2

+0

私は単純にその 'period'考えていましたある測定値と他の測定値との間の期間を参照している(したがって、私の場合は「1」)。しかし、ええ、これは明らかに多くの意味があります。私が理解したところで、私はnumpyでfftに 'i * 2 * pi'を掛けて再現することもできます。 – TomCho

0

1ラジアンは、差分近似のためのかなり粗いストライドで、あなたはかなりよく、データセット内の周期の整数で

x = np.arange(-200,200,1) 
y = np.sin(np.pi/50*x) 

plt.plot(np.diff(y)/np.diff(x)) 
plt.plot(sp.diff(y,order=1, period=400)) 

マッチしたいはずです - しかし、私はのための正確な合理的を知りませんfftルーチンの期間/正規化

関連する問題