8

I try to validate my understanding of Numpy's FFT with an example: the Fourier transform of exp(-pi*t^2) should be exp(-pi*f^2) when no scaling is applied on the direct transform.

However, I find that to obtain this result I need to multiply the result of FFT by a factor dt, which is the time interval between two sample points on my function. I don't understand why. Can anybody help ?

Here is a sample code:

# create data
N = 4097
T = 100.0
t = linspace(-T/2,T/2,N)
f = exp(-pi*t**2)

# perform FT and multiply by dt
dt = t[1]-t[0]
ft = fft(f)  * dt      
freq = fftfreq( N, dt )
freq = freq[:N/2+1]

# plot results
plot(freq,abs(ft[:N/2+1]),'o')
plot(freq,exp(-pi*freq**2),'r')
legend(('numpy fft * dt', 'exact solution'),loc='upper right')
xlabel('f')
ylabel('amplitude')
xlim(0,1.4)
4

1 に答える 1

17

連続時間フーリエ変換を計算していないことに注意してください。コンピューターは離散データを処理します。Numpy もそうです。numpy.fft.fftドキュメントを参照すると、次のように書かれています。

numpy.fft.fft(a、n=なし、軸=-1)[ソース]

1 次元の離散フーリエ変換を計算します。

この関数は、効率的な高速フーリエ変換 (FFT) アルゴリズムを使用して、1 次元の n 点離散フーリエ変換 (DFT) を計算します。

これは、次の式で定義される DFT を計算していることを意味します。

ここに画像の説明を入力

連続時間フーリエ変換は次のように定義されます。

ここに画像の説明を入力

そして、それらの間の関係を探すために計算を行うと、次のようになります。

ここに画像の説明を入力

ご覧の1/Nとおり、正確にスケール値である定数係数がありますdt(x[n] - x[n-1]ここで、n は [0,T] 間隔にあり、 に相当し1/Nます)。

コードにコメントするだけです。from numpy import *代わりにすべてをインポートすることはお勧めできません:

import numpy as np
import matplotlib.pyplot as plt

# create data
N = 4097
T = 100.0
t = np.linspace(-T/2,T/2,N)
f = np.exp(-np.pi*t**2)

# perform FT and multiply by dt
dt = t[1]-t[0]
ft = np.fft.fft(f) * dt      
freq = np.fft.fftfreq(N, dt)
freq = freq[:N/2+1]

# plot results
plt.plot(freq, np.abs(ft[:N/2+1]),'o')
plt.plot(freq, np.exp(-np.pi * freq**2),'r')
plt.legend(('numpy fft * dt', 'exact solution'), loc='upper right')
plt.xlabel('f')
plt.ylabel('amplitude')
plt.xlim([0, 1.4])
plt.show()

ここに画像の説明を入力

于 2013-11-14T11:12:05.810 に答える