##fftnとifftnだけ覚えておけば良さそう
numpy.fftには、いくつか関数があるのだが、そのうちfftnとifftnを覚えておけば良い。
fft(1次元)、fft2(2次元)などあるが、fftnが何次元でもいける。
ifftnが逆フーリエ変換。
例えば、2次元のフーリエ変換は以下のように使える。
import numpy as np
from numpy.fft import fftn
nx = 64
xmax = 2.*np.pi
xmin = 0.
x = np.linspace(xmin,xmax,nx,endpoint=False)
y = x.copy()
X, Y = np.meshgrid(x,y,indexing="ij")
z = np.sin(2.*X) + np.cos(Y)
fz = fftn(z)/nx/nx
print(np.mean(z**2))
print(np.sum(fz*np.conj(fz)))
とても簡単である。
fz = fftn(z,axes=0)
とすれば、０番目の軸に沿ってのみフーリエ変換できる。
特定の２軸以上をフーリエ変換したいときは以下のようにする。
規格化についてのオプションもあるがあとで説明する。
fz = fftn(z,axes=(0,1))
##番号と波数の関係
ここはIDLと同じ。
念のために配列の数が偶数と場合と奇数の場合に分ける。
0番目に、波数0の波が入るのは奇数でも偶数でも変わらない。
配列の個数をnx個とする。
奇数の場合
1から(nx-1)/2までに正の波数の値が入り、(nx-1)/2+1からnx-1までに負の波数の値が入る。
負の波数の方は大きい方から順番に値が入る。
例えば配列の要素が5個のときは以下のようになる。
| 配列 | 0 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| 波数 | 0 | 1 | 2 | -2 | -1 |
| 配列の要素数が奇数のときは、正の波数の数と負の波数の数は同数になる。 |
偶数の場合
1からnx/2までに正の波数の値が入り、nx/2+1からnx-1までに負の波数の値が入る。
例えば配列の要素が5個のときは以下のようになる。
| 配列 | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| 波数 | 0 | 1 | 2 | 3 | -2 | -1 |
| 配列の要素数が偶数のときは、正の波数の数が負の波数の数より１つ大きくなる。 |
波数0を配列の真ん中にもってくる
IDLには、centerというオプションがあって、波数0を配列の真ん中に持ってくるというものがあったがnumpyにはない。
これはnp.rollという配列の場所をshiftしてくれる関数を使うと良い。
import numpy as np
from numpy.fft import fftn
nx = 128
ny = 128
z = np.random.rand(nx,ny)
fz = fftn(z)
ffz = np.roll(fz,[(nx-1)//2,(ny-1)//2],axis=(0,1))
上記のように使う。負の波数の数だけシフトすればよい
奇数の時：(nx-1)/2個
偶数の時：(nx-2)/2個
(nx-1)//2と書けば整数計算で切り捨てなので、これで望む分をシフトすることができる。
シフト後は、
- 偶数の時はnx/2 - 1
- 奇数の時は(nx-1)/2
が波数０となる。
##フーリエ変換の規格化
二次元デカルト座標(x,y)のデータについて、フーリエ変換して1次元の分布を求める時の手順を記す。
データ領域の大きさをそれぞれ$L_x$, $L_y$。データ点数を$n_x$、$n_y$とする。
パーセバルの定理より、フーリエ変換の規格化は以下のようになると良い。
$
\displaystyle{\frac{\int v^2 dxdy}{L_xL_y} = \int\tilde{v}\tilde{v}^* dk_xdk_y}$
$2\pi$のファクターがあるのだが、ここは太陽業界の慣例に習ってつけないこととする。
ここで、$\tilde{v}$はフーリエ変換後の値(離散フーリエ変換)、$*$は複素共役を意味する。
一方、numpyのFFTの規格化は以下のようになっている。
$\displaystyle{\frac{\displaystyle{\sum_{i,j}}v^2}{n_xn_y}}
=\displaystyle{\frac{\displaystyle{\sum_{i,j}}\hat{v}\hat{v}^*}{n_x^2n_y^2}}
$
となるようになっている。
高速フーリエ変換をするということで、格子間隔は等幅なので下記の等式は成り立つ。
$
\langle v^2\rangle = \displaystyle{\frac{\int v^2 dxdy}{L_xL_y}} =
\displaystyle{\frac{\displaystyle{\sum_{i,j}}v^2}{n_xn_y}}
$
また、次元のある波数は以下のように表すことができる。
$
k_x = \displaystyle{\frac{2\pi}{L_x}}\hat{k}_x
$
ここで
$
\hat{k}_x
$は上記の表で示した無次元の波数である。よって
$
\displaystyle{dk_x = \frac{2\pi}{L_x}}
$
である。
$
\displaystyle{\int\tilde{v}\tilde{v}^* dk_xdk_y} =
\frac{(2\pi)^2}{L_xL_y}\sum_{i,j}\tilde{v}\tilde{v}^* =
\displaystyle{\frac{\displaystyle{\sum_{i,j}}\hat{v}\hat{v}^*}{n_x^2n_y^2}}
$
よって２次元の分布のままでは
$
\tilde{v}\tilde{v}^*=
\displaystyle{\frac{L_xL_y}{(2\pi)^2n_x^2n_y^2}}
\hat{v}\hat{v}
$
と規格化するとパーシバルの定理を満足する(最後の$v$はアステリスクがつくが、なぜか編集できないので、このままにしておく)
次に、 ２次元から1次元のスペクトルを求めたいときにどうすれば良いかを考える。
ここでは、$L_x=L_y$, $n_x=n_y$の場合を考える。
$
\displaystyle{
\int\int \tilde{v}\tilde{v}^* dk_xdk_y = \int\int \tilde{v}\tilde{v}^* k dkd\theta
}
$
とデカルト座標から極座標への変換をおこなう。
$
E_\mathrm{k}(k)=
\displaystyle{
\int \tilde{v}\tilde{v}^* k d\theta
}
$
と定義する。同じ円周上で合計することで積分のほとんどができている。$k$の次元が$2\pi/L_x$なので
このファクターをかける。
円周上に合計するときは、numpyのwhere関数を使う。ファクターを含めると以下のようにできる
nx = 64
ny = 64
ix = np.arange(-(nx-2)//2,nx//2+1)
jx = np.arange(-(ny-2)//2,ny//2+1)
ixx, jxx = np.meshgrid(ix,jx,indexing="ij")
nrr = np.sqrt(ix**2 + jx**2)
fz0 = np.zeros(nx//2,dtype=np.complex128)
for i in range(0,nx//2):
a = np.where(((nrr >=k) & (nrr < k+1)))
fz0[k] = fz[a[0],a[1]].sum(axis=0)
などとする。時間軸があった場合も最後のfor文を
for i in range(0,nx//2):
a = np.where(((nrr >=k) & (nrr < k+1)))
fz0[k,:] = fz[a[0],a[1],:].sum(axis=0)
などとする。