math - numpy による離散化された連続フーリエ変換

Question

関数 f(t) を考えてみましょう。連続フーリエ変換 g(w) を計算してプロットするにはどうすればよいですか (numpy と matplotlib を使用)。

フーリエ積分の解析解が存在しない場合、この問題または逆問題 (g(w) が与えられ、f(t) のプロットが不明) が発生します。

Answer 1

そのためにnumpy FFT モジュールを使用できますが、追加の作業を行う必要があります。まず、フーリエ積分を見て離散化します。 ここで、k,m は整数で、N は f(t) のデータ ポイントの数です。この離散化を使用して、

最後の式の合計は、numpy が使用する離散フーリエ変換 (DFT) とまったく同じです ( numpy FFT モジュールの「実装の詳細」セクションを参照してください)。この知識があれば、次の python スクリプトを書くことができます。

import numpy as np
import matplotlib.pyplot as pl

#Consider function f(t)=1/(t^2+1)
#We want to compute the Fourier transform g(w)

#Discretize time t
t0=-100.
dt=0.001
t=np.arange(t0,-t0,dt)
#Define function
f=1./(t**2+1.)

#Compute Fourier transform by numpy's FFT function
g=np.fft.fft(f)
#frequency normalization factor is 2*np.pi/dt
w = np.fft.fftfreq(f.size)*2*np.pi/dt


#In order to get a discretisation of the continuous Fourier transform
#we need to multiply g by a phase factor
g*=dt*np.exp(-complex(0,1)*w*t0)/(np.sqrt(2*np.pi))

#Plot Result
pl.scatter(w,g,color="r")
#For comparison we plot the analytical solution
pl.plot(w,np.exp(-np.abs(w))*np.sqrt(np.pi/2),color="g")

pl.gca().set_xlim(-10,10)
pl.show()
pl.close()

結果のプロットは、スクリプトが機能することを示しています