フーリエ変換

ある種の物理現象で用いられる線形の微分方程式では $f(t)=A_1e^{j\omega t}$ という解の形を想定し、得られた特性方程式を代数的に解くことで解けることが 多い。この解の形式$e^{j\omega x}$は、三角関数でも記述でき物理現象のうち 周期 (-p,p)G の時系列データf(t ) を以下のようにフーリエ級数で表現すること ができる。

$$f(t) = \frac{a_0}{2}+a_1 \cos\frac{\pi t}{p} +a_1 \cos\frac{\pi t}{p} +a_2 \cos\frac{2\pi t}{p}+\cdots +a_n \cos\frac{n\pi t}{p}+\cdots$$

$$+a_1 \sin\frac{\pi t}{p} +a_2 \sin\frac{2\pi t}{p}+\cdots +a_n \sin\frac{n\pi t}{p}+\cdots$$ (5.20)

$$= \sum_{n=-\infty}^{\infty}c_ne^{nj\pi t/p}$$ f(t ) (5.21)

この式における a_n, b_n, c_n をオイラー・フーリエの係数と呼ぶ。

Dirichletの条件
不連続な点を含む有界な周期関数f(t ) がある時、フーリエ級数は連続な部分で はf(t ) に収束し、不連続な点では不連続点の左右の極限値の平均値に収束する。

周期中不連続点を含む関数

この係数は、

$$= \frac{2}{p}\int_{d-p}^{d+p}f(t)dt$$ (5.22)

$$= \frac{1}{p}\int_{d-p}^{d+p}f(t)\cos\frac{n\pi t}{p}dt$$ (5.23)

$$= \frac{1}{p}\int_{d-p}^{d+p}f(t)\sin\frac{n\pi t}{p}dt$$ (5.24)

または

$$= \frac{1}{p}\int^p_{-p}f(t)dt$$ c₀ (5.25)

$$= \frac{a_n - jb_n}{2}\hspace{2ex}\mbox{(at $n<$\ 0)},\ \ = \frac{a_n + jb_n}{2}\hspace{2ex}\mbox{(at $n\ge$\ 0)}$$ c_n

$$= \frac{1}{2p}\int^p_{-p}f(t)e^{-nj\pi t/p}dt$$ (5.26)

となる。

同時に周期が無限大、すなわち$p=\infty$という考え方を持ち込むことで、周期 的でない関数にもこのフーリエ展開を適用することが可能になる。現象を記述す る関数に対してこの概念を入れると、 $\sin(n\pi/p)/n\pi/p$を用いて展開する ことも可能になる。これをフーリエ積分という。

デジタルフーリエ変換では、 対象とする有限個数の時系列データは$p=\infty$の近似である。 高速化するアルゴリズムにはいくつかあるが、 Cookey と Tukey のアルゴリズムがはじめてと言われている。 高速化のために、 データ数に２の基数を用いることから４や８の基数を用いる場合もある。 また、離散サイン変換(DST)や離散コサイン変換(DCT)という FFT よりも高速なアルゴリズムも作られている。 ここでは一般のデジタルフーリエ変換に関する注意を書く。

周期を NΔt とする 式(5.21)のフーリエ変換は

$$f_n = \cdots + c_k W_n^k + c_{k+1} W_n^{k+1} + \cdots,\qquad W_n = e^{j\frac{n\pi}{N}}$$

である。 $W_n^k=e^{j\frac{n\pi}{N}k}$では角速度 $\omega=n\Delta\omega = n\frac{\pi}{N\Delta t}$である。 つまり、角速度は分解能 Δω の整数倍 n となっている。 このとき、 ω という周波数成分が含まれている信号を デジタルフーリエ変換することを考える。

$$n\Delta\omega < \omega < (n+1)\Delta\omega$$

となる分解能の整数倍に一致しない角速度ω (周波数を fとすると$\omega=2\pi f$) を含む信号では、

$$c_n = \frac{1}{N} \sum_{k=0}^{N-1} f(t) W_n^{-k}$$

となるが、このf(t ) に含まれるωを考える。 $W_n^k=e^{j\omega_n k\Delta t}$とすると、

$$C=\frac{1}{N} \sum_{k=0}^{N-1} Ce^{jwt} W_n^{-k} = \frac{1}{... ...pprox \frac{1}{2p} \int_{-p}^{p} Ce^{j(\omega - \omega_n) t}dt$$

であり、 $\omega=\omega_n$とするとこの式の値はCとなる。 そして、 $\omega\ne \omega_n$のときには式の値は 0 となる必要がある。

そこで、 $\omega_n<\omega$として $\delta=\omega - \omega_n$と表し、 $e^\delta=1+\delta +\frac{\delta^2}{2}+\frac{\delta^3}{6}\cdots$ であることを用いて、上の式は

$$C^{*}= \frac{C}{2pj\delta} (e^{j\delta p}-e^{-j\delta p}) = ... ...dots\right\} =C\left\{1-\frac{(\delta p)^2}{6}+\cdots \right\}$$

となる。pは時系列の標本範囲であるので、ゼロではないので、 δがゼロでない場合には、 係数Cが正しく求めることができないことを示している。 積分区間 pが長い場合には $C^{*}$は大きくCからは外れることになる。 例として、 分解能が1Hzの処理系で、 49.3Hzと49.7Hzをそれぞれ1と2の大きさを持つ波形を デジタルフーリエ変換した結果を 図に示した。

多くの周波数成分を含む原信号を ディジタル化された 周波数を用いた変換した場合、 その周波数とその周波数近傍の成分も重畳してくるので 決して正しい係数を導出してはくれないので、 DFT(FFT)を用いる時には注意が必要となる。

これらの関係から、 $\pi/p\rightarrow 0$として、 $\sum\rightarrow\int$の変 換を行うと一対のフーリエ変換を定義できる。

$$= \int_{-\infty}^{\infty}g(\omega)e^{j \omega t}d\omega$$ f(t ) (5.27)

$$g(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}f(s)e^{-j \omega s}ds$$ (5.28)

このように、フーリエ変換は時系列のデータf(t ) を周波数空間データ $g(\omega)$へ写像するものであり、$g(\omega)$からf(t ) への変換は逆 フーリエ変換である。$g(\omega)$は複素数^5.3で あり、 $g(\omega)=a(\omega)+jb(\omega)$として、

$$\mbox{パワー or 電力} P(\omega) = g(\omega)g^*(\omega) = a^2(\omega) + b^2(\omega)$$ (5.29)

$$\angle(\omega) = \tan^{-1}a(\omega)/b(\omega)$$ 位相 (5.30)

として周波数 $\omega=2\pi f$のもつエネルギと位相を表すことができる。

以下に FFT の Sande-Tukey アルゴリズムを C で示す。

prog/fft/fft.prog

このアルゴリズムではバタフライ演算とビット逆順を用いている。この詳細につ いては解説書に譲る。