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

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

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

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

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

における $u$ の値

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

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

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

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

$$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) \mbox{($k\in\Z$)} .$$ (4)

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

簡単に証明できる良く知られた事実

• 普通の Fourier 係数

  $$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 \mbox{($k\in\Z$, $k\ge 0$)}$$ (5)

  
  
  と複素 Fourier 係数との間に、

  $$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.$$ (6)

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

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

  
  

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

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

  
  

• $\{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$ のとき

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

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

  $$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}).$$ (10)

  
  

Mathematica では、 長さ $N$ のリスト u={ $u_1,\dots,u_N$} があるとき、

v=Fourier[u,FourierParameters->{$a$,$b$}]

とすると、

$$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)$$   $$\mbox{($s=1,2,\dots,N$)}$$

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

特に

  v=Fourier[u,FourierParameters->{-1,-1}]

とすると、

$$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}$$   $$\mbox{($s=1,2,\dots,N$)}$$$$.$$

すなわち

   v$$=\{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}]