3.5 Fourier[] による離散Fourier変換

(考えてみたら、周期を一般の $ T$ とした式を書くべきだな…)

Mathematica には、 離散フーリエ変換をするための Fourier[ ] がある。

周期 $ 2\pi$ の周期関数 $ u\colon\R\to\C$ があるとき、 区間 $ [0,2\pi]$ $ N$ 等分点

(1) $\displaystyle t_j=\frac{2\pi(j-1)}{N}$   $\displaystyle \mbox{($j=1,2,\cdots,N$)}$

における $ u$ の値

(2) $\displaystyle u_j:=u(t_j)$   $\displaystyle \mbox{($j=1,2,\cdots,N$)}$

が得られたとする。このとき、 $ u$ の複素 Fourier 係数

(3) $\displaystyle c_k:=\frac{1}{2\pi}\int_0^{2\pi} u(t)e^{-i k t}\;\Dt$   $\displaystyle \mbox{($k\in\Z$)}$

を台形則で近似したものを $ C_k$ とすると、

(4) $\displaystyle C_k=\frac{1}{2\pi}\sum_{j=1}^{N} u(t_j)e^{-i k t_j}\cdot\frac{2\pi}{N} =\frac{1}{N} \sum_{j=1}^{N}u_j\exp\left(-\frac{2\pi i(j-1)k}{N}\right)$   $\displaystyle \mbox{($k\in\Z$)}$$\displaystyle .$

InverseFourier[] という逆変換が用意されている。

簡単に証明できる良く知られた事実
  • 普通の Fourier 係数

    (5) $\displaystyle a_k:=\frac{1}{2\pi}\int_{0}^{2\pi} u(t)\cos kt\;\Dt,\quad b_k:=\frac{1}{2\pi}\int_{0}^{2\pi} u(t)\sin kt\;\Dt$   $\displaystyle \mbox{($k\in\Z$, $k\ge 0$)}$

    と複素 Fourier 係数との間に、

    (6) $\displaystyle c_0=\frac{a_0}{2},\quad c_k=\left\{ \begin{array}{ll} \dfrac{a_k-...
...$)} [1ex] \dfrac{a_{-k}+i b_{-k}}{2} & \mbox{($k\le -1$)} \end{array} \right.$

    という関係がある ($ \{a_k\}$ , $ \{b_k\}$ から $ \{c_k\}$ を求めるのに使える)。 言い替えると (こちらは $ \{c_k\}$ から $ \{a_j\}$ , $ \{b_k\}$ を求めるのに使える)

    (7) $\displaystyle a_0=2c_0,\quad a_k=c_k+c_{-k},\quad b_k=i(c_k-c_{-k}) \quad ($$\displaystyle \mbox{$k\in\N$}$$\displaystyle ).$

    特に $ u$ が実数値関数の場合、$ \{a_k\}$ $ \{b_k\}$ は実数列であるから、 $ c_{-k}=\overline{c_k}$ が成り立ち (特に $ c_0$ は実数)、

    (8) $\displaystyle a_0=2 c_0=2\MyRe c_0,\quad a_k=2\MyRe c_k,\quad b_k=-2\MyIm c_k$   $\displaystyle \mbox{($k\in\N$)}$$\displaystyle .$

  • $ \{C_k\}$ は周期 $ N$ の周期数列である: $ C_{k+N}=C_k$ ($ k\in\Z$ ).
  • $ \vert k\vert\ll N/2$ のとき、$ C_k$ $ c_k$ の ``良い近似'' である。 $ 0\le k\ll N/2$ のとき

    (9) $\displaystyle c_k\simeq C_k,\quad c_{-k}\simeq C_{-k}=C_{N-k}.$

    ($ a_k$ , $ b_k$ の近似が必要な場合) $ 1\le k\ll N/2$ のとき、

    (10) $\displaystyle a_0\simeq 2 C_0,\quad a_k\simeq C_k+C_{-k}=C_k+C_{N-k},\quad b_k\simeq i(C_k-C_{-k})=i(C_k-C_{N-k}).$

Mathematica では、 長さ $ N$ のリスト u={ $ u_1,\dots,u_N$ } があるとき、
v=Fourier[u,FourierParameters->{$ a$ ,$ b$ }]
とすると、

$\displaystyle v_s:=\frac{1}{N^{(1-a)/2}}\sum_{r=1}^N u_r
\exp\left(\frac{2\pi i b (r-1)(s-1)}{N}\right)$   $\displaystyle \mbox{($s=1,2,\dots,N$)}$

で定まる $ v_1$ , $ \dots$ , $ v_{N}$ を並べたリスト v が得られる。 \fbox{\texttt{,FourierParameters->\{$a$,$b$\}}} を省略して 単に v=Fourier[u] とすると、 $ (a,b)=(0,1)$ と指定したのと同じになる (いわゆる ``default'')。

特に
  v=Fourier[u,FourierParameters->{-1,-1}]
とすると、

$\displaystyle v_s:=\frac{1}{N}\sum_{r=1}^N u_r
\exp\left(-\frac{2\pi i(r-1)(s-1)}{N}\right)=C_{s-1}$   $\displaystyle \mbox{($s=1,2,\dots,N$)}$$\displaystyle .$

すなわち

   v$\displaystyle =\{C_0,\cdots,C_{N-1}\}.
$

三角級数の係数を再生
a0 = 1

a = {1, 2, 3, 4, 5}

b = {6, 7, 8, 9, 10}

u[t_]:=a0/2+a . Table[Cos[n t],{n,Length[a]}]+b . Table[Sin[n t],{n,Length[b]}]

n = 12

ts = Table[2 Pi (j - 1)/n, {j, n}];

us = u[ts];

cs = Fourier[us, FourierParameters -> {-1, -1}]

A0 = 2 cs[[1]]

A = Table[cs[[k + 1]] + cs[[n - k + 1]], {k, 5}]

B = I Table[cs[[k + 1]] - cs[[n - k + 1]], {k, 5}]

桂田 祐史
2016-11-16