next up previous contents index
Next: 信号処理によるノイズ除去 Up: 計測のためのバックグランドとしての数学 Previous: 自己相関関数   目次   索引


フーリエ変換

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

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

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


\begin{wrapfigure}[4]{r}{.38\textwidth}
\begin{center}\vskip-1em
\includegraphi...
...}
\mfcaption{周期中不連続点を含む関数}
\end{center}\end{wrapfigure}
この係数は、

$\displaystyle a_0$ $\textstyle =$ $\displaystyle \frac{2}{p}\int_{d-p}^{d+p}f(t)dt$ (5.22)
$\displaystyle a_n$ $\textstyle =$ $\displaystyle \frac{1}{p}\int_{d-p}^{d+p}f(t)\cos\frac{n\pi t}{p}dt$ (5.23)
$\displaystyle b_n$ $\textstyle =$ $\displaystyle \frac{1}{p}\int_{d-p}^{d+p}f(t)\sin\frac{n\pi t}{p}dt$ (5.24)

または
$\displaystyle c_0$ $\textstyle =$ $\displaystyle \frac{1}{p}\int^p_{-p}f(t)dt$ (5.25)
$\displaystyle c_n$ $\textstyle =$ $\displaystyle \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)}$  
  $\textstyle =$ $\displaystyle \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 のアルゴリズムがはじめてと言われている。 高速化のために、 データ数に2の基数を用いることから4や8の基数を用いる場合もある。 また、離散サイン変換(DST)や離散コサイン変換(DCT)という FFT よりも高速なアルゴリズムも作られている。 ここでは一般のデジタルフーリエ変換に関する注意を書く。

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

\begin{displaymath}
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}}
\end{displaymath}

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

\begin{displaymath}
n\Delta\omega < \omega < (n+1)\Delta\omega
\end{displaymath}

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

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

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

\begin{displaymath}
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
\end{displaymath}

であり、 $\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$ であることを用いて、上の式は

\begin{displaymath}
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\}
\end{displaymath}

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

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

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

$\displaystyle f(t)$ $\textstyle =$ $\displaystyle \int_{-\infty}^{\infty}g(\omega)e^{j \omega t}d\omega$ (5.27)
$\displaystyle g(\omega)$ $\textstyle =$ $\displaystyle \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.2で あり、 $g(\omega)=a(\omega)+jb(\omega)$として、
$\displaystyle \mbox{パワー or 電力}$ $\textstyle P(\omega)$ $\displaystyle = g(\omega)g^*(\omega) = a^2(\omega) + b^2(\omega)$ (5.29)
$\displaystyle \mbox{位相}$ $\textstyle \angle(\omega)$ $\displaystyle = \tan^{-1}a(\omega)/b(\omega)$ (5.30)

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

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

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


next up previous contents index
Next: 信号処理によるノイズ除去 Up: 計測のためのバックグランドとしての数学 Previous: 自己相関関数   目次   索引
Ken Kishimoto 2014-06-02