デジタル信号処理:離散フーリエ変換
フーリエ変換 で説明したように,可積分関数 \(f\) のフーリエ変換は
\begin{align} F(t) = \int_{-\infty}^{\infty} f(x) e^{-itx} dx \label{eq:fourier-transform} \end{align} と定義される。また,フーリエ逆変換
\begin{align} f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(t) e^{itx} dt \label{eq:fourier-inverse-transform} \end{align}
により,関数 \(F\) から元の関数 \(f\) を得ることができる。しかし,フーリエ変換・フーリエ逆変換は連続関数に対する変換処理であるため,このままではコンピュータによるデジタル信号処理で使うことができない。
デジタル信号処理では,アナログ信号に対し標本化(サンプリング)を行うことでデジタル信号に変換する。つまり,デジタル信号は離散信号である。なお,この標本化の最大周波数の2倍よりも高い周波数(周期)で標本化をすることで元の波形に戻すことができるというシャノン・染谷の標本化定理により,アナログ信号のうち処理対象の最大周波数の2倍よりも高い周波数で標本化する。例えば,CD 音源では人の可聴域が 20 kHz 以下であるため 20 kHz を最大周波数としており,この2倍よりも大きい 44.1 kHz で標本化を行っている。デジタル信号処理では,このようにアナログ信号をデジタル信号に変換したものを処理している。
上述のフーリエ変換は可積分関数,つまり連続信号に対する変換であるため,このままの形では離散信号であるデジタル信号に対して適用できない。離散信号 \(x_n = \{ x_0, x_1, \dots, x_{N-1} \}\) に対するフーリエ変換(時間領域の有限個の離散信号から周波数領域の信号への変換)は離散フーリエ変換 (DFT: Discrete Fourier Transform) として以下のように定義される。
\begin{align} X_n = \sum_{k=0}^{N-1} x_k e^{-\frac{2\pi}{N}ikn} \label{eq:dft} \end{align}
また,周波数領域から時間領域への変換である逆離散フーリエ変換は
\begin{align} x_k = \frac{1}{N} \sum_{n=0}^{N-1} X_n e^{\frac{2\pi}{N}ikn} \label{eq:idft} \end{align}
と表される。
導出
畳み込みの定理
離散フーリエ変換を導出するにあたり,後述の離散時間フーリエ変換の周期性の説明で畳み込みの定理(Convolution Theorem)を用いるため,まず畳み込みの定理について説明する。畳み込みとは,2つの関数 \(f(x)\) と \(g(x)\) を組み合わせたもので,\(f(x)\) と \(g(x)\) の畳み込みを \((f * g)(x)\) と表すと,以下のように定義される。
\begin{align*} (f * g)(x) = \int_{-\infty}^{\infty} f(y) g(x - y) dy \end{align*}
コンピュータサイエンスに慣れ親しんだ方は * が積でないこと(つまり \(f(x) g(x)\) ではないこと)に注意してください。畳み込みには,
- 交換則
- 結合則
- 分配則
が成り立つ。
フーリエ変換における畳み込みの定理とは,関数 \(\cdot\) に対するフーリエ変換を \(\mathcal{F}\{\cdot\}\) と表したとき,
\begin{align*} \mathcal{F}\{f * g\} = \mathcal{F}\{f\} \mathcal{F}\{g\} \end{align*}
が成り立つことを言う。畳み込みの定理の説明では使用しないが,この流儀に従ってフーリエ逆変換を \(\mathcal{F^{-1}}\{\cdot\}\) と表すこともある。
ここではまず,畳み込みの定理を導出する。
\begin{align*} \mathcal{F}\{f * g\} &= \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} f(y) g(x - y) dy \right) e^{-itx} dx \\ &= \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} f(y) g(x - y) e^{-itx} dy \right) dx \\ &= \int_{-\infty}^{\infty} f(y) \left( \int_{-\infty}^{\infty} g(x - y) e^{-itx} dx \right) dy \end{align*}
と書ける。ここで \(z = x - y\) とおくと,
\begin{align*} \mathcal{F}\{f * g\} &= \int_{-\infty}^{\infty} f(y) \left( \int_{-\infty}^{\infty} g(z) e^{-it(y+z)} dz \right) dy \\ &= \int_{-\infty}^{\infty} f(y) e^{-ity} \left( \int_{-\infty}^{\infty} g(z) e^{-itz} dz \right) dy \\ &= \mathcal{F}\{f\} \mathcal{F}\{g\} \end{align*}
となり,畳み込みの定理が導かれる。
また,関数 \(F, G\) をそれぞれ, 関数 \(f, g\) のフーリエ変換とすると,関数 \(f, g\) の積 \((fg)\) のフーリエ変換は,
\begin{align} \mathcal{F}\{fg\} &= \int_{-\infty}^{\infty} f(x) g(x) e^{-itx} dx \notag \\ &= \int_{-\infty}^{\infty} f(x) \left( \frac{1}{2\pi} \int_{-\infty}^{\infty} G(z) e^{ixz} dz \right) e^{-itx} dx \notag \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty} G(z) \left( \int_{-\infty}^{\infty} f(x) e^{-i(t-z)x} dx \right) dz \notag \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty} G(z) F(t-z) dz \notag \\ &= \frac{1}{2\pi} (G * F)(t) \notag \\ &= \frac{1}{2\pi} \left( \mathcal{F}\{f\} * \mathcal{F}\{g\} \right) \label{eq:convolution-fg} \end{align}
と変形できる。
デルタ関数とデルタ関数列
連続信号の離散化にあたり,ディラックのデルタ関数が広く用いられている。デルタ関数とは
\begin{align*} \delta(x) = \begin{cases} \infty && (x = 0) \\ 0 && (x \neq 0) \end{cases} \end{align*}
となり,積分すると \( 1 \)になる関数である。
関数 \( f(x) \) とデルタ関数の積を積分すると
\begin{align*} \int_{-\infty}^{\infty} f(x) \delta(x) dx = f(0) \end{align*}
となる性質がある。
デルタ関数のこの性質を用い,連続信号に対し周期的にデルタ関数をかけることで標本化(サンプリング)できる。つまり,デルタ関数を周期 \(T\) で並べた関数(デルタ関数列)を
$$ \delta_T(x) = \sum_{k=-\infty}^{\infty} \delta(x - kT) $$
と定義すると,連続信号 \(f(x) \)に対する周期 \(T\) でサンプリングした関数は
$$ f(x) \delta_T(x) = \sum_{k=-\infty}^{\infty} f(kT) \delta(x - kt) $$
と表せる。
ここでデルタ関数列 \( \delta_T \) のフーリエ変換を考える。\( \delta_T \) は周期関数であるため,フーリエ級数展開により
\begin{align} \delta_T(x) &= \sum_{k=-\infty}^{\infty} \left( \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} \delta(x) e^{-\frac{2\pi}{T}ikx} dx \right) e^{\frac{2\pi}{T}ikx} \notag \\ &= \frac{1}{T} \sum_{k=-\infty}^{\infty} e^{\frac{2\pi}{T}ikx} \label{eq:delta0} \end{align}
と書ける。これを用い,
\begin{align} \mathcal{F}\{\delta_T(x)\} &= \mathcal{F}\left\{ \frac{1}{T} \sum_{k=-\infty}^{\infty} e^{\frac{2\pi}{T}ikx} \right\} \notag \\ &= \frac{1}{T} \sum_{k=-\infty}^{\infty} \mathcal{F}\left\{ e^{\frac{2\pi}{T}ikx} \right\} \label{eq:delta1} \end{align}
\begin{align*} \mathcal{F^{-1}}\{\delta(t - p)\} = \frac{1}{2\pi} \int_{-\infty}^{\infty} \delta(t - p) e^{itx} dt = \frac{e^{ipx}}{2\pi} \end{align*}
より,
\begin{align*} \mathcal{F}\left\{ e^{\frac{2\pi}{T}ikx} \right\} &= 2\pi \mathcal{F}\left\{ \mathcal{F^{-1}}\left\{ \delta(t - \frac{2\pi k}{T}) \right\} \right\} \\ &= 2\pi \delta(t - \frac{2\pi k}{T}) \end{align*}
となるので,式\ref{eq:delta1}にこれを代入して,
\begin{align*} \mathcal{F}\{\delta_T(x)\} &= \frac{2\pi}{T} \sum_{k=-\infty}^{\infty} \delta(t - \frac{2\pi k}{T}) \end{align*}
となる。\(\omega = \frac{2\pi}{T} \) とおくと,
\begin{align} \mathcal{F}\{\delta_T(x)\} &= \omega \sum_{k=-\infty}^{\infty} \delta(t - k\omega) \notag \\ &= \omega \delta_\omega(t) \label{eq:fourier-delta-t} \end{align}
となる。
離散時間フーリエ変換
ここまでにデルタ関数列のフーリエ変換まで扱った。次に時間領域を離散化した信号のフーリエ変換を考える。前述したディラックのデルタ関数列を用い,連続関数 \( f \) をサンプリングし,離散信号を得ることができる。関数 \(f\) に対し,サンプリング周期 \(T\) でサンプリングした関数は \( f(x) \delta_T(x) \) と連続信号に対するインパルス応答の和の関数として表せる。この離散信号のフーリエ変換 \( \mathcal{F}\{f(x)\delta_T(x)\} \) は,式\ref{eq:fourier-transform} より,
\begin{align*} \mathcal{F}\{f(x)\delta_T(x)\} &= \int_{-\infty}^{\infty} \left( \sum_{k=-\infty}^{\infty} f(kT) \delta(x - kT) \right) e^{-itx} dx \\ &= \sum_{k=-\infty}^{\infty} f(kT) \int_{-\infty}^{\infty} \delta(x - kT) e^{-itx} dx \\ &= \sum_{k=-\infty}^{\infty} f(kT) e^{-ikTt} \end{align*}
となる。ここで,\( f[n] = f(nT) \) とすると,
\begin{align} \mathcal{F}\{f(x)\delta_T(x)\} &= \sum_{k=-\infty}^{\infty} f[k] e^{-ikTt} \label{eq:dtft} \end{align}
となる。式\ref{eq:dtft} を離散時間フーリエ変換 (DTFT: Discrete-Time Fourier Transform) という。式\ref{eq:dtft} において \(T\) は定数なので \( Tt \) を1つの変数として表す書き方もできる。
離散時間フーリエ変換は,連続関数 \( f \) のフーリエ変換を \( F(t) = \mathcal{F}\{f(x)\} \) とすると,この \( F(t) \) を使って表すことができる。式\ref{eq:convolution-fg}を用い,\( F(t) = \mathcal{F}\{f(x)\} \) とすると,
\begin{align*} \mathcal{F}\{f(x)\delta_T(x)\} &= \frac{1}{2\pi} \left( \mathcal{F}\{f(x)\} * \mathcal{F}\{\delta_T(x)\} \right) \\ &= \frac{1}{2\pi} \left(F(t) * \omega \delta_\omega(t) \right) \quad (\because \text{式\ref{eq:fourier-delta-t}}) \\ &= \frac{1}{T} \left(F(t) * \delta_\omega(t) \right) \end{align*}
と元の連続関数 \( f \) のフーリエ変換と周期 \( \omega (=\frac{2\pi}{T}) \) のデルタ関数列のフーリエ変換の畳み込みの \( \frac{1}{T} \) 倍と表せる。
これをさらに変形すると,
\begin{align*} \mathcal{F}\{f(x)\delta_T(x)\} &= \frac{1}{T} \int_{-\infty}^{\infty} F(x) \left( \sum_{n=-\infty}^{\infty} \delta(t - n \omega - x) \right) dx \\ &= \frac{1}{T} \sum_{n=-\infty}^{\infty} F(t - n \omega) \end{align*}
と表すことができる。つまり,離散時間フーリエ変換で得られるスペクトルは,元の信号をフーリエ変換して得られるスペクトルを周期 \(\omega\) (の整数倍)ずらしたものを足し合わせて \( \frac{1}{T} \) 倍したものである。つまり,\( \omega \) が元の信号をフーリエ変換して得られるスペクトラムの最大周波数の2倍(\(F(t)\) の幅の倍)以上であれば,このずらしたスペクトラムが重ならないため,離散時間フーリエ変換をして得られるスペクトラムから元のスペクトラムを復元できる。このようにして,前述のシャノン・染谷の標本化定理が証明できる。
離散フーリエ変換
前述した離散時間フーリエ変換は,区間 \([-\infty, \infty]\) の信号が必要であり,また式\ref{eq:dtft} からも分かる通り,フーリエ変換後の関数(つまり,周波数領域において)は連続関数となる。コンピュータで現実の信号を扱えるようにするには,有限区間かつ周波数領域でも離散値としたい。
上述の通り,離散時間フーリエ変換をして得られるスペクトラムは周期性を持つ。フーリエ級数が離散的なスペクトラムの和で表されたように,入力信号が周期性を持つとき周波数領域も離散化できる。つまり,周期関数をフーリエ(逆)変換すると離散値となり,離散化した関数をフーリエ(逆)変換すると周期性を持つ。
そこで,\( f \) を区間 \( [0, NT] \) の周期関数とすることで,有限区間の信号として扱う。つまり,サンプリング周期 \( T \) で \( N \) 個の信号をサンプリングした離散信号を周波数領域に変換する。
まず,式\ref{eq:dtft} で \( \mathcal{F}\{f(x)\delta_T(x)\} = F(t),\ Tt = \frac{2\pi}{N} s \) とおくと,
\begin{align} \mathcal{F}\{f(x)\delta_T(x)\} = F\left(\frac{2\pi}{NT} s \right) = \sum_{k=-\infty}^{\infty} f[k] e^{-\frac{2\pi}{N}iks} \label{eq:dft01} \end{align}
となる。\( f(x) \) が周期 \( NT \) の周期関数であるため,
$$ \forall m,n \in \mathbb{Z}, \; f[n - mNT] = f[n] $$
とかける。よって,式\ref{eq:dft01} を以下のように変形できる。
\begin{align} F\left(\frac{2\pi}{NT}s\right) &= \sum_{m=-\infty}^{\infty} \sum_{l=0}^{N-1} f[l - mN] e^{-\frac{2\pi}{N}i (l - mN) s} \notag \\ &= \sum_{l=0}^{N-1} f[l] e^{-\frac{2\pi}{N}i l s} \sum_{m=-\infty}^{\infty} e^{2\pi i m s} \label{eq:dft02} \end{align}
ここで,式\ref{eq:delta0} より,
$$ \sum_{k=-\infty}^{\infty} e^{2 \pi ikx} = \delta_1 (x) $$
であるため,式\ref{eq:dft02} は,
\begin{align*} F\left(\frac{2\pi}{NT}s\right) &= \sum_{l=0}^{N-1} f[l] e^{-\frac{2\pi}{N}i l s} \delta_1(s) \\ &= \sum_{l=0}^{N-1} f[l] e^{-\frac{2\pi}{N}i l s} \left(\sum_{k=-\infty}^{\infty} \delta(s - k) \right) \end{align*}
と離散化できる。つまり,\(F\) は \( s \) が整数のときだけ値を持ち,それ以外のときは \( 0 \) となる。ここで,\( F\left(\frac{2\pi}{NT}s\right) = F[s] \) とおくと,
\begin{align} F[n] &= \sum_{k=0}^{N-1} f[k] e^{-\frac{2\pi}{N}i k n} \label{eq:dft-proof} \end{align}
となり,離散時間フーリエ変換から離散フーリエ変換が導けた。
離散フーリエ変換は,離散時間フーリエ変換ではなくフーリエ級数から導くこともできる。\(f(x)\delta_T(x)\) のフーリエ係数 \( c_n \) は \( \epsilon \) を正の実数とする(つまり区間 \([0, NT)\) の定積分で考える)と,
\begin{align*} c_n &= \lim_{\epsilon \rightarrow 0} \frac{1}{NT} \int_{0}^{NT-\epsilon} \left( f(x) \sum_{k=0}^{N-1} \delta(x - kT) \right) e^{-\frac{2 \pi}{NT}inx} dx \\ &= \lim_{\epsilon \rightarrow 0} \frac{1}{NT} \sum_{k=0}^{N-1} f(kT) \int_{0}^{NT-\epsilon} \delta(x - kT) e^{-\frac{2 \pi}{NT}inx} dx \\ &= \frac{1}{NT} \sum_{k=0}^{N-1} f[k] e^{-\frac{2\pi}{N}ikn} \end{align*}
となる。ここで,
$$ F[n] = \sum_{k=0}^{N-1} e^{-\frac{2\pi}{N}ikn} $$
は,式\ref{eq:dft-proof} であり,離散フーリエ変換の式と一致する。フーリエ級数展開より,
\begin{align*} f(x)\delta_T(x) &= \sum_{n=-\infty}^{\infty} c_n e^{\frac{2\pi}{NT}inx} \\ &= \frac{1}{NT} \sum_{n=-\infty}^{\infty} F[n] e^{\frac{2\pi}{NT}inx} \\ &= \frac{1}{NT} \sum_{m=-\infty}^{\infty} \sum_{l=0}^{N-1} F[l+mN] e^{\frac{2\pi}{NT}i(l+mN)x} \quad (l,m \in \mathbb{Z}) \\ &= \frac{1}{NT} \sum_{l=0}^{N-1} F[l] e^{\frac{2\pi}{NT}ilx} \sum_{m=-\infty}^{\infty} e^{\frac{2\pi}{N}imx} \\ &= \frac{1}{N} \sum_{l=0}^{N-1} F[l] e^{\frac{2\pi}{NT}ilx} \delta_T(x) \quad (\because \text{式\ref{eq:delta0}}) \end{align*}
よって,\(x = kT \) を代入して,
\begin{align} f[k] = \frac{1}{N} \sum_{n=0}^{N-1} F[n] e^{\frac{2\pi}{N}ikn} \label{eq:idft-proof} \end{align}
と書ける。式\ref{eq:idft-proof} は式\ref{eq:idft} の逆離散フーリエ変換と一致しており,フーリエ級数展開から逆離散フーリエ変換の式が得られる。このように,フーリエ級数展開から離散フーリエ変換と逆離散フーリエ変換を導ける。