1 離散フーリエ変換

ここでの離散フーリエ変換の説明については、参考文献 [1]を大いに 参考にした。

1.1 区間 ${\boldsymbol [0,\,2\pi]}$ の場合

数学で学習するフーリエ変換は、連続的な関数を対象にする。しかし実際の測定量、例え ば電圧などは離散データであることの方が圧倒的に多い。そこで、このような離散的なデー タを計算機でフーリエ変換することを考える。これを離散フーリエ変換と言い、コンピュー ターが得意とする分野である。フーリエ変換という名前が付いているが、実際には フーリエ級数に近い。これからの説明は、フーリエ級数の展開係数を求めたことを思 い出しながら読めば理解しやすいであろう。

フーリエ級数では周期が重要であったが、離散フーリエ変換ではデータの区間がそれ に対応する。任意の区間については後ほど述べることとし、まずは話 を簡単にするためにそれが$2\pi$の場合を考える。0〜$2\pi$秒の間で電圧を測定した ようなものである。離散フーリエ変換の結果をこの区間を越えて拡張する場合²、区間 $2\pi$ と同じ 周期が現れる。

この区間で $N$ 個の等間隔でデータが得られた考える。等間隔という ことは重要で、それを

$$x_j =\frac{2\pi}{N}j,\qquad(j=0,\,1,\,2,\,\cdots,\,N-1)$$ (1)

と表す。一方、従属変数の方は、

$$f_j=f(x_j)$$ (2)

と表現することにする。実験データでは、$x_j$ は時刻、$f_j$ が電圧というようになる。 オシロスコープやAD変換器を用いてデータを取得したような状況を考えればよい。

ここでは、測定量 $(x_j,\,f_j)$ の組は $N$ 個ある。この $N$ 個のデータを操作して、 得られる1次独立なデータ数は最大 $N$ 個である。このことから、フーリエ係数は $N$ 個が限界 であることが分かる。従って、展開の基底も $N$ 個となり、それを

$$\left\{1,\,e^{ix},\,e^{2ix},\,e^{3ix},\,\cdots,\,e^{(N-1)ix}\right\}$$ (3)

のように選ぶ。基底の線形和として、元の関数 $f(x)$ は、

$$f(x) =c_0+c_1e^{ix}+c_2e^{2ix}+\cdots+c_{N-1}e^{(N-1)ix}$$

$$=\sum_{k=0}^{N-1}c_k e^{ikx}$$ (4)

として表すことができる。このように三角関数の和で表すことを、フーリエ級数と言うの である。ここで残された問題は、測定量の $x_j$ と $f_j$ から $c_k$ を決めるにはど のようにするかである。これは比較的簡単で、

$$f_j =\sum_{k=0}^{N-1}c_k e^{ikx_j}\qquad(j=0,1,2,\cdots,N-1)$$ (5)

の連立方程式を解けばよろしい。この式の形が分かりにくい。それではもう少 し分かり易く書いてみよう。

$$\begin{bmatrix}f_0\\ f_1\\ f_2\\ \vdots\\ f_{N-1} \end{bmatrix} =... ...} \end{bmatrix} \begin{bmatrix}c_0\\ c_1\\ c_2\\ \vdots\\ c_{N-1} \end{bmatrix}$$ (6)

$x_j$ は既知なので、この $N \times N$ の行列は計算可能である。左辺も既知である。 従って、この連立方程式から、未知数 $c_k$ を決めれば良い。ちまたでよく使われる高 速フーリエ変換(First Fourier Transform 略してFFT)は、この $c_k$ を高速に計算して いる。

「なるほど、連立方程式ならば計算機は得意なので、式(6)を解 けば話は終わり」と思ってしまうかもしれないが、そんなに世の中は甘くない。$N=1000$ 位になるとこの式を計算するのにかなりの時間がかかり実用に適さない。実際、FFTを使 いたい場面ではリアルタイム(少なくとも1秒くらい)で計測が終わる必要があるので、計 算の高速化を図らなくてはならない。そこで、計算量を減らすことを考える。

そのために、$N$ を整数、$\ell$ も整数であるがその範囲は $0,\,1,\,2,\,\cdots,\,N-1$ とした場合の

$$\frac{1}{N}\sum_{j=0}^{N-1}e^{i \ell x_j}= \begin{cases}0 & \ell\neq 0\text{のとき}\\ 1 & \ell=0\text{のとき}\\ \end{cases}$$ (7)

を使う。この式の計算については、3.1節を見よ。さて、ここから が本題で、式(5)の両辺に

$$\frac{1}{N}e^{-ik^\prime x_j}, \quad k^\prime=0,1,2,\cdots,N-1$$ (8)

を乗じて、$j$ について和を取る。すると、

$$\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-ik^\prime x_j} =\frac{1}{N}\sum_{j=0}^{N-1}\left(\sum_{k=0}^{N-1}c_k e^{ikx_j}\right)e^{-ik^\prime x_j}$$

$$=\frac{1}{N}\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}c_k e^{ix_j(k-k^\prime)}$$

$$k=k^\prime$$

$$=c_{k^\prime}$$ (9)

となる。これで、$c_k$ を求める式が得られた。もう少しきれいにまとめると、

$$c_k=\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-ikx_j}$$ (10)

である。この式が離散フーリエ変換(DFT)と呼ばれるものである。もう少し分 かりやすく書くと、

$$\begin{bmatrix}c_0\\ c_1\\ c_2\\ \vdots\\ c_{N-1} \end{bmatrix} =... ...} \end{bmatrix} \begin{bmatrix}f_0\\ f_1\\ f_2\\ \vdots\\ f_{N-1} \end{bmatrix}$$ (11)

である。

以上をまとめると、区間 $[0,\,2\pi]$ の場合の離散逆フーリエ変換が式 (4)、離散フーリエ変換が式(10)である。

1.2 区間 ${\boldsymbol [0,T]}$の場合

これまでは区間 $[0,\,2\pi]$ での離散フーリエ変換について説明してきたが、これからは もっと一般的な問題を取り扱うために、区間 $[0,\,T]$ について述べる。区間幅 $T$ で $[t_0,t_0+T]$ の方がより一般的と考えるかもしれないが、これは $t^\prime=t-t_0$ のよ うにシフトすれば $[0,\,T]$ となるので同じである。

独立変数を $t$ として、区間でのデータの取得は

$$t_j =\frac{T}{N}j$$

$$=\Delta t j,\qquad(j=0,\,1,\,2,\,\cdots,\,N-1)$$ (12)

とする。展開の基底を、

$$\left\{1,\,e^{i\frac{2\pi}{T}t},\,e^{2i\frac{2\pi}{T}t}, \,e^{3i\frac{2\pi}{T}t},\,\cdots,\,e^{(N-1)i\frac{2\pi}{T}t}\right\}$$ (13)

のように選ぶ。ここで、 $\omega=\frac{2\pi}{T}$ とする。もし、測定間隔が時間の場合、 この $\omega$ は角振動数となる。このようにすると、基底は

$$\left\{1,\,e^{i\omega t},\,e^{2i\omega t}, \,e^{3i\omega t},\,\cdots,\,e^{(N-1)i\omega t}\right\}$$ (14)

と書き改められる。ここで、$1$ は直流成分を表し、あとは等間隔に周波数成分が現れる ことになる。データ点の場所では、直流成分と $N\omega$ は区別がつかないことに注意が 必要である³。同様に、角振動数が $N\omega$ シフトした 場合も区別が付かない。

基底の線形和として、元の関数$f(t)$は、

$$f(t) =c_0+c_1e^{i\omega t}+c_2e^{2i\omega t}+\cdots+c_{N-1}e^{(N-1)i\omega t}$$

$$=\sum_{n=0}^{N-1}c_n e^{in\omega t}$$ (15)

となる。これは、フーリエ逆変換である。今は離散データを問題としているので、 $f(\Delta t j)$ を $f_j$ と書くのが良いだろう。すると、離散フーリエ逆変換は

$$f_j=\sum_{n=0}^{N-1}c_n e^{in\omega t_j}$$ (16)

と書き表せる。この $C_n$を求めれば良いのである。これを前節と全く同じ方法で求める。 $\omega t_j=\frac{2\pi}{N}j$ なので前節の $x_j$ と同一である。従って、

$$\frac{1}{N}\sum_{j=0}^{N-1}e^{i \ell \omega t_j}= \begin{cases}0 & \ell\neq 0\text{のとき}\\ 1 & \ell=0\text{のとき}\\ \end{cases}$$ (17)

である。式(16)の両辺に $\frac{1}{N}e^{-in^\prime\omega t_j}$ を乗じて、$j$ で和をとると、

$$\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-in^\prime\omega t_j} =\frac{1}{N}\sum_{j=0}^{N-1}\left( \sum_{n=0}^{N-1}c_n e^{in\omega t_j} \right)e^{-in^\prime\omega t_j}$$ (18)

$$=\frac{1}{N}\sum_{n=0}^{N-1}\sum_{j=0}^{N-1} c_ne^{i(n-n^\prime)\omega t_j}$$

$$=c_{n^\prime}$$

となる。従って、この場合の離散フーリエ変換は

$$c_n=\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-in\omega t_j}$$ (19)

となる。

1.3 まとめと実際の計算

これまで示した離散フーリエ変換の結果をまとめて、整理しておく。区間 $[0,\,2\pi]$の場合は、

$$c_k=\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-ikx_j} 離散フーリエ変換$$ (20)

$$f_j=\sum_{k=0}^{N-1}c_k e^{ikx_j} 逆離散フーリエ変換$$ (21)

である。区間 $[0,\,T]$の場合は

$$c_n=\frac{1}{N}\sum_{j=0}^{N-1}f_je^{-in\omega t_j} 離散フーリエ変換$$ (22)

$$f_j=\sum_{n=0}^{N-1}c_n e^{in\omega t_j} 逆離散フーリエ変換$$ (23)

である。ただし、 $\omega=2\pi/T$である。

データが $\left[\{0,\,f_0\},\,\{\Delta t,\,f_1\},\,\{2\Delta t\,f_2\},\,\cdots,\,\{(N-1)\Delta t,\,f_{N-1}\}\right]$ の場合の離散フーリエ変換 のプログラムを作成する式を示しておくのが良いだろう。そこで、区間 $[0,\,T]$ の 式を $\omega=2\pi/T$ と $t_j=\frac{T}{N}j$ を使って、

$$c_n=\frac{1}{N}\sum_{j=0}^{N-1}f_j e^{-2\pi ijn/N} 離散フーリエ変換$$ (24)

$$f_j=\sum_{n=0}^{N-1}c_n e^{2\pi ijn/N} 逆離散フーリエ変換$$ (25)

と書き改める。コンピューターで離散フーリエ変換を計算する場合は、これらの式を計算すべき である。 $c_n$ は $n\omega$ の振幅、 $f_j$ は$j\Delta t$ の時の値を表している。
ホームページ: Yamamoto's laboratory
著者: 山本昌志