ガウス関数のFFTについて考察

ガウス関数は正規分布の書き方だと，平均$x_0$と分散$\sigma$を用いて

f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp \Bigl({-\frac{(x-x_0)^2}{2 \sigma^2}} \Bigr)

となります．今回はこのガウス関数のフーリエ変換について，理論式を計算した結果とFFT(高速フーリエ変換)の結果を比較します．

#目次
・理論式によるフーリエ変換
・FFTによるフーリエ変換
・両者の比較
・FFT実行のためのコード(Python)

##理論式によるフーリエ変換
計算過程は省略しますが，ガウス関数のフーリエ変換$F$はまたガウス関数になります．

\begin{eqnarray}
F(k)
&=& \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{-ikx}dx \\
&=& \sigma e^{-\frac{\sigma^2k^2}{2}}e^{-ikx_0} \\ \\
\end{eqnarray}

振幅スペクトル$|F|=\sigma e^{-\frac{\sigma^2k^2}{2}}$と，
位相(偏角の主値)${\rm Arg}(F)={\rm Arg}(e^{-ikx_0})$を$k \geq 0$で描画します．

上が振幅スペクトル，下が位相($/\pi$)です．フーリエ変換と言えば振幅スペクトルの事を指すことが多い気がします．位相は$-kx_0$と書けますが，今回は$[-\pi \ \pi]$の主値を採用しているので，端を越えると下にジャンプします．また，振幅スペクトルが落ち込んだ先の位相情報は，意味が無いので無視します．

(ところでこのグラフは横軸が波数なので，そのまま角周波数に読み替えることが出来ますが，振動数で読むためには$2\pi$で割る必要があることに注意．)

##FFTによるフーリエ変換
FFTの実行方法の説明は他に譲るとして，ここでは理論式とFFTの結果が一致するのかを検証します．自分はPythonのNumPyモジュールで実行しました(参考コードを最後に記載しておきます)．
先程と同様の関数をFFTした結果から得られる，振幅スペクトルと位相のグラフを示します．

スペクトルの最大値が先程と違いますが，概形は似たものになりました．ちなみにデータ数を多めに取らないと，きれいな結果になりませんでした．

##両者の比較
二つの結果が一致しているのか検証します．
振幅スペクトルは比を，位相は差を取って比較しました．
理論値とFFTで得た値をそれぞれ下付き添字1,2で表しています．

$k \leq 4$では比が一定，差が0という結果になりました．FFTで求まる振幅スペクトルの最大値は異なっていますが，全体が定数倍されているだけなので，定性的には問題無いと言えそうです．
というわけで，ガウス関数のフーリエ変換の理論式とFFTの結果を比較し，一致することが確認できました．

####FFT実行のためのコード(Python)
参考までに，今回FFTを実行するにあたって用いたコードを掲載しておきます．

import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl

sigma = 1.0
dx = 0.01
x0 = 5
x = np.arange(0, 1000, dx)
N = len(x)

def y(x):
    return 1 / (np.sqrt(2 * np.pi) * sigma) \
           * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

fourier = np.fft.fft(y(x))
fourier_abs = np.abs(fourier)
fourier_phase = np.arctan2(np.real(fourier), np.imag(fourier))
omega = np.linspace(0, 1/dx, N) * 2 * np.pi

fig = plt.figure()

ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.plot(k, np.abs(fourier), linewidth=1, color='black')
ax1.set_ylabel('Amplitude')

ax2.plot(k, np.arctan2(np.real(fourier), np.imag(fourier)) / np.pi, linewidth=1, color='black')
ax2.set_xlabel('k')
ax2.set_ylabel('Phase[rad]')

ax1.set_xlim([0.0, 6.0])
ax2.set_xlim([0.0, 6.0])

mpl.rcParams['axes.xmargin'] = 0
mpl.rcParams['axes.ymargin'] = 0

plt.show()