nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

ガウス関数のフーリエ変換 (scipy.fftpack)

ガウス関数フーリエ変換しても、数式上は(幅は変わるが)またガウス関数になる。
高速フーリエ変換FFT)とフーリエ変換との関係を以前まとめたため、FFTを用いてガウス関数フーリエ変換する。
koideforest.hatenadiary.com

ここではscipy.fftpackを用いてFFTを行うが、(少なくともこの場合)出力される配列データの並びが「(正の振動数のデータ)、(負の振動数のデータ)」となっており、そのままプロットするとガウス関数に見えない。

試しにやってみる。

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

a = 1
f = lambda x: np.exp( - a * x**2 )  # Gauss function

max_x  = 10.
min_x  = -10.
N      = 2**7
x = np.linspace( min_x, max_x, N )
data   = f( x )
delta_x  = x[1] - x[0]  # ( max_x - min_x ) / ( N - 1 )

# Fast Fourier Transformation 
tilde_f = fft( data ) * delta_x

delta_k = 1. / delta_x / ( N - 1 )
k = np.arange( 0, 1./delta_x + delta_k, delta_k ) * 2. * np.pi

plt.plot( k,  abs( tilde_f ), "o" )
plt.show()

f:id:koideforest:20180916224622p:plain

データの並びに関しては、以下のサイトが参考になる。
fft - memoring
SciPy FFTpack
Fourier Transforms (scipy.fftpack) — SciPy v1.1.0 Reference Guide

もう少し説明すると、FFTにおけるアルゴリズムの性質上、fft で生成した配列は、直感的な[ -N/2, ..., N/2 ]という順番では無く、

[ 0, 1, ,2, ... (N/2 - 1), N/2, -N/2, (-N/2 + 1), ..., -1 ]

という風になっている。
つまり、

[ (正の振動数:昇順)、(負の振動:降順)]

である。

fftpackに内蔵されている fftfreq を使うと、fft と同様な並びの横軸を作成してくれる。

from scipy.fftpack import fftfreq 

freq = fftfreq( N, delta_x ) * 2. * np.pi

plt.plot( freq, abs( tilde_f ), "-o" )
plt.show()

f:id:koideforest:20180916224702p:plain

点でプロットするならば問題無いが、線を追加すると N/2 → -N/2 に移るときの線が現れてしまう。
そこで、配列の並びを直感的に戻してくれる関数 fftshift が用意されている。

from scipy.fftpack import fftshift

tilde_f_shift = fftshift( tilde_f )
freq_shift = fftshift( freq )

gauss_k = lambda k: np.sqrt( np.pi / a ) * np.exp( - k**2 / ( 4. * a ) )

plt.plot( freq_shift, abs( tilde_f_shift ), "-o" )
put.plot( freq_shift, gauss_k( freq_shift ), "x" )
plt.show()

f:id:koideforest:20180916224802p:plain
解析的に求めた解も一緒にプロットした。

\displaystyle
f( x ) = e^{ - a x^2 }
\\
\displaystyle
f( k ) = \sqrt{ \frac{ \pi }{ a } } e^{ - \frac{ k^2 }{ 4 a }}


解析解の導出は以下のサイトが分かり易い。ただし、指数の係数 aは省略されている。
ガウス関数のフーリエ変換 | ポップラーン

しかし、FFTを用いたガウス関数フーリエ変換で最も問題なのは、実は数値解が振動している点である。
虚部は実部に比べて小さいため、実部のみプロットすると、

plt.plot( freq_shift, tilde_f_shift.real, "-o" )
plt.plot( freq_shift, abs( tilde_f_shift ), "-" )
plt.show()

f:id:koideforest:20180916230025p:plain

何故、振動が現れるのかは、今の所パッと思い付かないが、フーリエ変換はなかなか侮れないことを教えてくれていると思う。