Python(numpy)の高速フーリエ変換の使い方と規格化について

##fftnとifftnだけ覚えておけば良さそう

numpy.fftには、いくつか関数があるのだが、そのうちfftnとifftnを覚えておけば良い。
fft(1次元)、fft2(2次元)などあるが、fftnが何次元でもいける。
ifftnが逆フーリエ変換。

例えば、2次元のフーリエ変換は以下のように使える。

```python3
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)))

```

とても簡単である。

```python3
fz = fftn(z,axes=0)
```
とすれば、０番目の軸に沿ってのみフーリエ変換できる。

特定の２軸以上をフーリエ変換したいときは以下のようにする。
規格化についてのオプションもあるがあとで説明する。

```python3
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してくれる関数を使うと良い。

```python3
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$とする。

[パーセバルの定理](https://ja.wikipedia.org/wiki/%E3%83%91%E3%83%BC%E3%82%BB%E3%83%90%E3%83%AB%E3%81%AE%E5%AE%9A%E7%90%86)より、フーリエ変換の規格化は以下のようになると良い。

$
\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関数を使う。ファクターを含めると以下のようにできる

```python3

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文を

```python3
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)
```
などとする。