離散フーリエ変換 (DFT)

\(\)
基準座標の式 (8-77) は, 「デジタル信号処理」(Digital Signal Processing : DSP) に於ける「離散フーリエ変換」(DFT) の式と全く同じ形をしていることに気付いた.従って, 基準座標 \(Q_{\alpha}\) は変位 \(q_j\) のフーリエ変換と見做すことも出来るかも知れない.そこで参考のために, 小畑:「信号処理」の第2章から DFT に関係する部分を抜粋したものに多少の修正を施したものを示しておこう.


サンプリング定理

ディジタル信号処理に於いては, 信号 \(x(t)\) を一定時間ごとにサンプリングし, そのサンプル値

\begin{equation}
\dotsb,x(-T),x(0),x(T),x(2T),\dotsb,x(nT),\dotsb
\tag{1}
\end{equation}

を数値に変換してディジタル演算を実行する.この場合, まず時間軸 \(t\) がサンプリング周期 \(T\) によって離散化され, 次に信号値 \(x(nT)\) がアナログ・ディジタル変換の量子化単位 \(\varepsilon\) によって離散化される.通常, アナログ・ディジタル変換は実用上十分な精度のものとするので数値化したデータは信号値 \(x(nT)\) をほぼ正確に表現すると考えてよい.ここで問題にするのは, 時間軸 \(t\) を離散化するサンプリング過程, すなわち \(x(nT)\) に変換する過程である(図 1 を参照).

図 1.

信号 \(x(t)\) は連続した値をとる時間 \(t\) の関数として与えられている.そこから一定時間 \(T\) 刻みの値 \(x(nT)\) だけを取り出して使用することにすると, 元の信号 \(x(t)\) の持っていた情報の一部が失われるのではないか?, 従って, サンプリング過程が不可避であるディジタル信号処理を採用すると近似的な結果しか求めることが出来ないのではないか?, との疑問が発生する.しかし, これは実際問題としては正しくない.我々が取り扱う信号は, 通常, ある周波数以下の成分しか含んでいないので, この性質を利用すると, サンプル値 \(x(nT)\) だけ使用して \(x(t)\) のスペクトル \(X(j\omega)\) を正確に求めることが可能である.ただし,「サンプリング周期 \(T\)」は上限周波数から定まるある値以下としなければならない.しかし,「サンプリング時間 \(T\)」を無限に小さくする必要はない.

もちろん, 実際問題としては, 完全に正確なスペクトル \(X(j\omega)\) を求めるのは困難である.その原因はサンプリング周期 \(T\) が有限であるためではなく, 利用できる信号長が有限であるためである.我々は信号 \(x(t)\) を \(t=-\infty\) から \(t=\infty\) までの無限長の時間にわたって観測することは出来ず, ある有限長 \(L\) のデータ \(x(t)\) \(( -\frac{L}{2}\le t \le \frac{L}{2}) \) しか利用できないからである.

理論的取り扱いとしては, \(x(t)\) \((-\infty<t<\infty)\) 或いは, \(x(nT)\) \( (-\infty<n<\infty)\) が利用できるとすると便利である.この場合, \(x(nT)\) だけから完全に正確なスペクトル \(X(j\omega)\) が求まるが, 信号の時間領域表現 \(x(t)\) と周波数領域表現 \(X(j\omega)\) は完全に等価であるから, \(X(j\omega)\) が正確に求まることは \(x(t)\) も正確に求まることを意味する.すなわち, \(x(nT)\) から任意の時刻 \(t\) の信号値 \(x(t)\) を計算, 復元することが出来る.この事実は「サンプリング定理」として知られているもので, 染谷 勲C.E. Shannon によって独立に研究されたものである.

この信号値系列 \(x(nT)\) を, ある形の時間信号 \(x_s(t)\) として表現することが出来る.この信号 \(x_s(t)\) は \(t=nT\) 以外ではもちろんゼロの値をとる:

\begin{equation}
x_s(t)=0,\quad t\ne nT,\qquad n=0,\pm1,\pm2,\dotsb
\tag{2}
\end{equation}

従って, この信号 \(x_s(t)\) は「デルタ関数」を用いて次のように表現することが出来る:
\begin{align}
x_s(t)&=\dotsb+x(-T)\delta(t+T)+x(0)\delta(t)+x(T)\delta(t-T)+\dotsb\notag\\
&=\sum_{n=-\infty}^{\infty} x(nT)\delta(t-nT)
\tag{3}
\end{align}

このとき,
\begin{align}
\int_{-\infty}^{\infty} x_s(t)\,dt
&= \int_{-\infty}^{\infty}dt\,\left\{\sum_{n=-\infty}^{\infty} x(nT)\delta(t-nT) \right\}\notag\\
&= \sum_{n=-\infty}^{\infty}x(nT)\int_{-\infty}^{\infty}\delta(t-nT)\,dt
=\sum_{n=-\infty}^{\infty}x(nT)
\tag{4}
\end{align}

この信号 \(x_s(t)\) のフーリエ変換を考える.それを \(X_s(j\omega)\) と書くと,
\begin{align}
X_s(j\omega)&=\int_{-\infty}^{\infty} x_s(t)e^{-j\omega t}\,dt \notag\\
&= \int_{-\infty}^{\infty}\left\{\sum_{n=-\infty}^{\infty} x(nT)\delta(t-nT) \right\}\,e^{-j\omega t}\,dt\notag\\
&=\sum_{n=-\infty}^{\infty} x(nT)\int_{-\infty}^{\infty} e^{-j\omega t}\,\delta(t-nT)\,dt\notag\\
&=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j\omega nT}
\tag{5}
\end{align}

もちろん, この \(X_s(j\omega)\) のフーリエ逆変換は \(x_s(t)\) となる.

サンプリング周期 \(T\)」から「サンプリング角周波数 \(\Omega_0\)」を次とする:

\begin{equation}
\Omega_0=\frac{2\pi}{T}
\tag{6}
\end{equation}

すると, \(k\) を任意の整数とするならば, 式 (5) から次が言える:
\begin{align}
X_s[j(\omega+k\Omega_0)]&=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j(\omega+k\Omega_0) nT}\notag\\
&=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j\omega nT}\,e^{-jkn\Omega_0T}
=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j\omega nT}\,e^{-j2\pi kn}\notag\\
&=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j\omega nT}=X_s(j\omega)
\tag{7}
\end{align}

従って,「\(X_s(j\omega)\) の周期は, サンプリング角周波数 \(\Omega_0\) に等しい」ことが分かる.このことと, 式 (5) から, \(x(nT)\) は \(X_s(j\omega)\) のフーリエ級数展開の係数になっていると言える:
\begin{equation}
x(nT)=\frac{1}{\Omega_0}\int_{-\Omega_0/2}^{\Omega_0/2} X_s(j\omega)\,e^{jn\omega T}\,d\omega
=\frac{T}{2\pi}\int_{-\Omega_0/2}^{\Omega_0/2} X_s(j\omega)\,e^{jn\omega T}\,d\omega
\tag{8}
\end{equation}

それでは,「原信号 \(x(t)\) のスペクトル \(X(j\omega)\)」は「サンプル信号 \(x_s(t)\) のスペクトル \(X_s(j\omega)\)」とどのような関係を持っているのであろうか?.\(x(t)\) のフーリエ変換が \(X(j\omega)\) であるので, その逆変換から
\begin{equation}
x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty} X(j\omega)\,e^{j\omega t}\,d\omega,\ \rightarrow\
x(nT)=\frac{1}{2\pi}\int_{-\infty}^{\infty} X(j\omega)\,e^{j\omega nT}\,d\omega
\tag{9}
\end{equation}

この \(\omega\) に関する \(-\infty\) から \(\infty\) までの積分を, 次のような区間 \(\omega_k\) に分割して積分することを考える:
\[ k\Omega_0 – \frac{\Omega_0}{2} \le \omega_k \le k\Omega_0 +\frac{\Omega_0}{2},\quad k= 0,\pm1,\,\pm2,\,\dotsb \]
すると, 式 (9) の \(x(nT)\) は次のように書ける:
\begin{equation}
x(nT)=\frac{1}{2\pi}\sum_{k}\int_{\omega_k} X(j\omega)\,e^{j\omega nT}\,d\omega
\tag{10}
\end{equation}

さらに変数変換 \(\omega’=\omega-k\Omega_0\) を行うと, \(-\frac{\Omega_0}{2}\le \omega’_k\le \frac{\Omega_0}{2},\,d\omega’=d\omega,\,\Omega_0T=2\pi\) となるから,
\begin{align}
x(nT)&= \frac{1}{2\pi}\sum_{k}\int_{-\Omega_0/2}^{\Omega_0/2} X(j\omega’+jk\Omega_0)\,e^{j(\omega’+k\Omega_0) nT}\,d\omega’
\notag\\
&=\frac{1}{2\pi}\int_{-\Omega_0/2}^{\Omega_0/2} \left\{\sum_{k}X(j\omega’+jk\Omega_0)\right\}\,e^{j\omega’nT}\,d\omega’
\tag{11}
\end{align}

この結果を式 (8) と比較すると,次式が言えることが分かる:
\begin{equation}
X_s(j\omega)=\frac{1}{T}\sum_k X(j\omega+jk\Omega_0)
\tag{12}
\end{equation}

従って, 原信号 \(x(t)\) のスペクトル \(X(j\omega)\) を \(0,\,\pm\Omega_0,\,\pm 2\Omega_0,\,\dotsb\) だけ位置をずらせて加算したものが \(TX_s(j\omega)\) になると言える.よって, 原信号 \(x(t)\) のスペクトル \(X(j\omega)\) が \(\Omega_0/2\) 以上の周波数成分を含まない場合には, \(-\frac{\Omega_0}{2}\le \omega\le \frac{\Omega_0}{2}\) の範囲において \(X(j\omega)\) と \(TX_s(j\omega)\) とは完全に一致する.しかしながら, \(X(j\omega)\) が \(\Omega_0/2\) 以上の周波数成分を含む場合には, それが \(\Omega_0/2\) 以下の周波数成分に重なり合うことになる.この現象を「エイリアシング」(aliasing) と言う (図 2 を参照).

fig2-7

図 2.

このエイリアシングが発生するのを避けるためには,「サンプリング周期 \(T\) を, 原信号 \(x(t)\) の最高周波数成分の周期の半分以下に設定することが必要」である.

離散フーリエ変換

信号値系列 \(x(nT)\) の時間表現信号を \(x_s(t)\) とするとき, そのフーリエ変換 \(X_s(j\omega)\) は式 (5) に示したように次で与えられる:

\begin{equation}
X_s(j\omega)=\sum_{n=-\infty}^{\infty} x(nT)\,e^{-j\omega nT}
\tag{5′}
\end{equation}

この \(X_s(j\omega)\) の数値計算法を以下で示す.前述したように \(X_s(j\omega)\) は \(\omega\) に関して周期 \(\Omega_0=2\pi/T\) の周期関数であるから, \(-\Omega_0/2\) から \(\Omega_0/2\) までの範囲, 或いは, \(0\) から \(\Omega_0\) までの範囲の \(X_s(j\omega)\) を求めれば十分である.
我々が実際に観測できるのは, ある有限な範囲の信号値である.それを, 次としよう:
\begin{equation}
x(0),\ x(T),\ \dotsb,\ x((N-1)T)
\tag{13}
\end{equation}

この範囲以外の信号値 \(x(nT)\) は不明であるが, 簡単化のためにそれらは全てゼロと置くことにするか, または, 上記の信号値を周期 \(N\) で繰り返す周期的な信号系列と見做してもよい.すると,
\begin{equation}
X_s(j\omega)= \sum_{n=0}^{N-1} x(nT)\,e^{-jn\omega T}
\tag{14}
\end{equation}

となる.観測した信号の時間長は \(NT\) なので, 最も低い角周波数成分としては
\begin{equation}
\omega_0=\frac{2\pi}{NT}
\tag{15}
\end{equation}

のものを計算することにし, \(\omega=k\omega_0\ (k=0,\,1,\,2,\,\dotsb,\,N-1)\) に於ける \(X_s(j\omega)\) を求めることにすると,
\begin{align}
X_s(k\omega_0)&= \sum_{n=0}^{N-1} x(nT)\,\exp\left(-jnk\frac{2\pi}{NT} T\right)\notag\\
&=\sum_{n=0}^{N-1} x(nT)\,\exp\left(-j\frac{2\pi}{N}kn\right),\quad k=0,\,1,\,2,\,\dotsb,\,N-1
\tag{16}
\end{align}

となる.これは数値計算可能である.\(k=N\) に於いては,
\[ \omega=N\omega_0=N\frac{2\pi}{NT}=\frac{2\pi}{T}=\Omega_0 \]
であるから, \(\omega=0\) から \(\omega=\Omega_0\) まで, \(\omega_0\) の刻みで \(X_s(j\omega)\) が定まることになる.式 (16) を信号値系列 \(\big\{x(nT)\big\}\) の「離散フーリエ変換」(Discrete Fourier Transformation, DFT) という.次式が成り立つことは, 容易に確かめられる (前述のブログ記事を参照すべし):
\begin{equation}
\frac{1}{N}\sum_{n=0}^{N-1} \exp\left(-j\frac{2\pi}{N}kn\right)\exp\left(j\frac{2\pi}{N}ln\right)
=\begin{cases} 1, & k-l が N の整数倍のとき \\ 0, & それ以外のとき\end{cases}
\tag{17}
\end{equation}

この関係を用いると,
\begin{equation}
x(nT)= \frac{1}{N}\sum_{k=0}^{N-1} X_s(k\omega_0)\,\exp\left(j\frac{2\pi}{N}kn\right), \quad n=0,\,1,\,2,\,\dotsb,\,N-1
\tag{18}
\end{equation}

この式 (18) を「離散フーリエ逆変換」(Inverse Discrete Fourier Transformation, IDFT) という.