本記事は窓関数の役割(端点処理)や離散フーリエ変換と連続関数版フーリエ変換の間の関係性だったりを主眼においた説明をしようとしたが、教科書に書いてあるような内容の記事になってしまった。より正確にはフーリエ変換の専門書を読んだほうが分かりやすいかもしれない。とはいえ折角記述してしまったので、記事にしておこう。
フーリエ変換と畳み込み
周期関数のフーリエ変換
周期Tの関数$f(x)=f(x+T)$は、
$$f(x)=\sum_{k=-\infty}^{\infty}\mathcal c_k\exp(i\frac{2\pi k}{T}x)$$
と、基底周波数の倍数となる周波数の振動関数を重ね合わせることで表現できる。このとき、
$$c_k=\frac{1}{T}\int_{-T/2}^{T/2}f(x)\exp(-i\frac{2\pi k}{T}x)dx$$
と計算できる。任意の周期関数は三角関数の重ね合わせで表現することができ、このスペクトルを求める操作をフーリエ変換と呼ぶ。
フーリエ変換した関数は周波数の依存したスペクトルとなるが、これを周波数空間上の関数とよぶ。
$c_k$が上記計算で出せることは、$k$を整数として
$$\int_{-T/2}^{T/2}\exp(ikx)dx=\begin{cases}
T \quad k\neq 0\
0 \quad k=0
\end{cases}$$
であることからわかる。これはクロネッカーのデルタ$\delta_{ij}$を用いて
$$\int_{-T/2}^{T/2}\exp(ikx)dx=T\delta_{k0}$$
とかける。
(証明)
$$\frac{1}{T}\int_{-T/2}^{T/2}\mathcal f(x) \exp(-i\frac{2\pi k}{T}x)dx=\frac{1}{T}\int\sum_{k'} c_{k'}\exp(i\frac{2\pi k'}{T}x)\exp(-i\frac{2\pi k}{T}x)dx$$
$$\therefore=\frac{1}{T}\sum_{k'} c_{k'}\int\exp(i\frac{2\pi (k'-k)}{T}x)dx=\frac{1}{T}\sum_{k'} c_{k'}T\delta_{kk'}=c_k$$
となることより、成り立つ。
非周期関数
なお、周期関数以外の関数においても、周期無限
の関数だと考えれば、
$$f(x)=\sum_k c_k \exp(i\frac{2\pi k}{T}x)$$
$$c_k=\frac{1}{T}\int_{-T/2}^{T/2}f(x)\exp(-i\frac{2\pi k}{T}x)dx$$
とできる。が、このままでは$T\rightarrow\infty$のときに$c_k$
がどうなるかがわかりにくいため、ちょっと式変形する。
$$w=\frac{2\pi k}{T}$$
$$c_k=\frac{2\pi}{T}F(w)$$
とすると、以下のように変形できる。
$$f(x)=\lim_{a\rightarrow \infty,T\rightarrow \infty}\sum_{k=-aT}^{aT}\frac{2\pi}{T}F(\frac{2\pi k}{T})\exp(i\frac{2\pi k}{T})$$
$w=\frac{2\pi k}{T}$として、区分求積法をすれば、
$$f(x)=\int_{-\infty}^{\infty}F(w)\exp(iwx)dw$$
となる。$v=\frac{k}{T}$とした場合も、
$$f(x)=\int_{-\infty}^\infty 2\pi F(2\pi v)\exp(2\pi ivx)dv$$
$dv=2\pi dw$として変換すれば、結局おなじである。
同様に、
$$c_k=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}f(x)\exp(-i\frac{2\pi k}{T}x)dx$$
$$\therefore F(w)=\lim_{T\rightarrow\infty}\frac{1}{2\pi}\int_{-T/2}^{T/2}f(x)\exp(-iwx)dx=\frac{1}{2\pi}\int_{-\infty}^{\infty}f(x)\exp(-iwx)dx$$
まとめると、非周期関数のフーリエ変換は
$$F(w)=\frac{1}{2\pi}\int f(x)\exp(-iwx)dw$$
$$f(x)=\int F(w)\exp(iwx)dw$$
なお、これは$F'(w)=2\pi F(w)$と置き直せば
$$F'(w)=\int f(x)\exp(-iwx)dx$$
$$f(x)=\frac{1}{2\pi}\int F'(w)\exp(iwx)dw$$
と書くことができる。同様に、$F'(w)=\sqrt{2\pi}F(w)$とすれば
$$F'(w)=\frac{1}{\sqrt{2\pi}}\int f(x)\exp(-iwx)dx$$
$$f(x)=\frac{1}{\sqrt{2\pi}}\int F'(w)\exp(iwx)dw$$
と対称性のよい定義にすることもできる。もっといえば、
$$F(w)=\int_{-\infty}^{\infty}f(x)\exp(iwx)dx$$
$$f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(w)\exp(-iwx)dx$$
という定義もある($w\rightarrow-w$と置き直した形式で、工学分野に多い)
このフーリエ変換と逆変換操作を$F(w)=\mathcal F(f(x)),f(x)=\mathcal F^\dagger(F(w))$と記述する。$\dagger$はエルミート共役である。
どれもスペクトルの強度基準をどこにおくかの違いであるが、どれも定義として利用され得るため変換と逆変換のペアは対応させて使わないと変換逆変換で元の波形の定数倍だけズレた結果となることがある。
ようは、得られたスペクトルはどの定義式で変換器されたものなのかを意識する必要がある。単にライブラリの逆フーリエ変換をテキトーに使っただけでは必ずしも望む結果にならない可能性がある点は、注意である。
同じようなことはDFTやFFTについてもいえる。(FFTはwikiの記事まんまになってしまうので、本記事ではDFTまでの紹介とする)
畳み込み、フィルター計算
信号処理の過程でフーリエ変換を挟む必要があるとすれば、その主な理由の一つになり得る操作が畳み込み計算だろう。フィルタ計算とも呼ばれることがある。
畳み込み計算の定義式は以下である。
$$f\ast h=\int_{-\infty}^{\infty}f(x-x')h(x')dx'$$
あるいは、
$$f\ast h=\int_{-\infty}^{\infty}f(x+x')h(x')dx'$$
ここで$f(x)$を入力関数、$h(x)$はカーネル関数またはフィルタカーネルと呼ぶ。
グラフでわかりやすくする便宜上、後者の定義で進める(カーネル関数を前後逆に作用させるかどうかの違いであって、大きな差はない)。
有名なフィルターの効果としては平滑化や
エッジ抽出などがあり、以下のような関数をカーネルに用いたときに効果が得られる。
たとえば、有名な平滑化フィルターとしてガウシアンカーネルというものがある。その関数形状は以下のように定義される。
$$h(x)=\frac{1}{Z}\exp(-x^2/\sigma^2)$$
※細かい係数が異なることがあるが、議論の大局には影響しないのでこのままとする。
なお$\sigma$は分散とか半径とか呼ばれることがあり、大きいほど遠くの画素を巻き込んだ平滑化の作用が得られる。言い換えると、大きいほど強い平滑化がかかる。
これの半径を小さくした上で三角関数と同期させたものは、微分フィルタとして機能する。
微分は以下の計算で瞬間的勾配を出すのであった。
$$\frac{df}{dx}=\lim_{\epsilon\rightarrow 0}\frac{f(x+\epsilon/2)-f(x-\epsilon/2)}{\epsilon}$$
これを$f\ast h=\int f(x+x')h(x)dx'$と表現したとき、hは
$$h(x)=\frac{\delta(x-e)-\delta(x+e)}{2e}$$
とかき、$e\rightarrow 0$としたときに該当する。
ところで、
$$\delta(x)=\lim_{\epsilon\rightarrow 0}\frac{1}{\epsilon\sqrt\pi}\exp(-x^2/\epsilon^2)$$
であるので、
$$h(x)=\frac{1}{2\epsilon^2\sqrt\pi}\lbrace\exp(-(x-\epsilon)^2/\epsilon^2)-\exp(-(x+\epsilon)^2/\epsilon^2)\rbrace$$
として、十分小さな$\epsilon$を設定すれば、フィルタ結果は勾配となるはずだ。やってみよう。
一番上の段が入力波形、二段目からはかけたフィルターを左図、入力波形と畳み込み後の波形を重ねたものが右図である。
一番下以外はガウシアンフィルタでσを変えたものであり、一番下の段は微分フィルタの結果である。効果をわかりやすくするため、最後のフィルタのみノイズを平滑化した区系
入力信号はガウシアンフィルタの径が小さいうちはノイズを平滑化する効果があるが、径が大きいと大きな構造まで平滑化してしまうことがわかる。
一方、微分フィルタの結果は勾配の向きに応じた値変化をもたらす意味で意図通りの結果になっているといえる。
周波数空間上での畳み込み処理
さきほどの畳み込み処理
$$f\ast h=\int f(x+x')h(x')dx'$$
とフーリエ変換の関係を調べよう。結論をいうと、
$$F(w)H^\ast(x)=\int (f\ast h)(x)\exp(iwx)dx=\mathcal F(f\ast h)$$
となることが知られている。導出は以下である:
(下記の各積分範囲は実数全範囲であるが省略する)
$$\int \int f(x+x')h(x')dx'\exp(iwx)dx=\frac{1}{4\pi^2}\int\int \int F(w_1)\exp(-iw_1(x+x'))dw_1\int H(w_2)\exp(-iw_2x')dw_2dx'dx$$
$$=\frac{1}{4\pi^2}\int\int\int \exp(i(w-w_1)x)dxF(w_1)\exp(-iw_1x')dw_1\int H(w_2)\exp(-iw_2x)dw_2dx'$$
$$=\frac{1}{2\pi}\int F(w)\int\exp(-i(w+w_2)x')dx'H(w_2)dw_2=F(w)H(-w)=F(w)H^\ast(w)$$
ここで、
$$\int_{-\infty}^{\infty} \exp(\pm i(w-w')x)dx=2\pi\delta(w-w')$$
$$\int \delta(w-w')F(w')dw'=F(w')$$
を用いて計算した(定数=1を出す関数のフーリエ変換は$w\neq0$ですべて0であり、フーリエ変換結果を逆フーリエ変換すると定数1となることからデルタ関数の定義を満たす)。
同様に、周波数空間上で畳み込みを定義することができ、
$$f(x)h^\ast(x)=\mathcal F^\dagger(F(w)\ast H(w))$$
である。つまり、畳み込みを
$$f(q) \ast h(q)=\int f(q+q')h(q)dq$$
と定義したとき、周波数空間上の畳み込みは実空間上では成分ごとの複素共役積になり、実空間上の畳み込みは周波数空間上での複素共役積となる。なお、
$$f(q) \ast h(q)=\int f(q-q')h(q)dq$$
と定義したときは復素共役積ではなく、通常の成分積になる。両者が一致するとき、畳み込むフィルタは偶関数である必要があり、畳み込むフィルターのフーリエ変換結果は実関数となる。
$$f(x)h(x)=\mathcal F(F(w)\ast H(w))$$
$$F(w)H(w)=\mathcal F(f(x)\ast h(x))$$
上記式を逆に、関数同士の成分積で作られた関数はフーリエ変換すると関数同士の畳み込み結果となることがいえる。この事実は次の節で重要となる。
平滑化フィルタと微分フィルターの周波数特性
平滑化フィルター$h_L$と微分フィルター$h_H$はそれぞれフーリエ変換すると以下のようになる。
$$\mathcal F(h_L)=\int_{-\infty}^{\infty}\exp(-x^2/\sigma^2)\exp(iwx)dx=\int\exp(-\frac{(x-i\sigma^2w/2)^2}{\sigma^2}-\frac{\sigma^2w^2}{4})dx\propto\exp(-\frac{\sigma^2w^2}{4})=C\exp(-Aw^2)$$
$$\mathcal F(h_H)=C\exp(-A(w\epsilon)^2)(\exp(iw\epsilon)-\exp(-iw\epsilon))=C'\exp(-A(w\epsilon)^2)\sin(w\epsilon)$$
(定積分は、ガウス積分と無限範囲の複素経路積分で計算できるが省略)
それぞれをグラフにすると、以下のようになる。
平滑化フィルターは$w=0$付近のみ0でない値をもち、微分フィルターは$w$が大きいところほど大きな値をもつ。
実空間上の畳み込み(フィルター)処理は周波数空間上でのかけ算になるため、カーネル関数のフーリエ変換において出力が0となる周波数に着目すると、畳み込み関数の出力も0となる。つまり畳み込みによって特定周波数のスペクトルがカットされる。平滑化フィルターのように周波数の低いところのみ情報が残り高周波がカットされるのでローパスフィルター、微分フィルターのように低い周波数をカットし高い周波数を残すフィルターをハイパスフィルターと呼ぶ。
波形の一部分の周波数解析と窓関数
音声などで得られる波形について、短時間の区間で周波数解析したいことがある。これは音声に含まれる音色や和音を解析する際に有用なてがかりになるはずだ。この方法について議論しよう。
- 音声の周期性を調べて、ちょうど一周期になるところを積分範囲として周期関数のフーリエ変換
- 音声データ全区間を非周期関数とみなして全域一括フーリエ変換
- 音声データを一定苦寒でフーリエ変換
ちょうど一周期の波形を特定する？
1については、ちょうど一周期の範囲を特定するのが難しい。音声は一般的に様々な音色の重ね合わせであり、たとえば以下の二つの振動の重ね合わせで生まれる波形の合成周波数はいくらでも長くなる。
$$y=\cos k_1x+\cos k_2x=\cos\frac{k_1+k_2}{2}x\cos\frac{k_1-k_2}{2}x$$
いま、$k_1=k-\epsilon,k_2=k+\epsilon$とすると
$$y=\cos kx\cos\epsilon x$$
つまり、二波長で構成される波形の合成周期は波長の差分に等しい長さになる。これをうなりという。
厳密には重ね合わせた波形の周波数は$k_1=\frac{n_1}{T},k_2=\frac{n_2}{T}$となるため、二つの周波数が無理数な波形の場合はそもそも周期が定義不能である。
ゆえに、音声の解析方法として本方法は不適だろう。
非周期関数として全域一括フーリエ変換？
これはある程度うまくいくが、全部の時間内に含まれる音をすべて一気に周波数化するので
特定の時間内しか鳴っていないような音色を識別することに向かない。遊びでフーリエ変換して、三角関数の重ね合わせで音声を復元できる不思議な事象をデモしたりする際に使われることがある。
一方、音声解析のモチベーションは普通、いつどんな音がなっているかを区別したい興味をもつことが多いため音声全体をフーリエ変換することも不適である。
一定区間でフーリエ変換？
これが音声解析のアプローチとしては一番ストレートだといえるだろう。が、これにも問題がある。
まずは以下の図をみよう。単なるsin/cos波を特定区域を切ってフーリエ変換したものだ。
ちょうど整数周期で切った結果は計算誤差を除いて想定通りの結果になっている。一方、少し周波数をずらしてみると、高周波(や低周波)に謎の持ち上げが出てくる。この理由は、何回かこのパターンを繰り返せば見える。
ようは、区間を切ってフーリエ変換するということは、切った区間が周期的に続く関数をフーリエ変換していることと同じ意味を持つ。冒頭の周期関数のフーリエ変換の処理そのものであるためである。
一つ目のグラフでは、両端が不連続なところが急なエッジだと解釈されるため、高周波成分としてカウントされてしまうのである。なぜエッジが高周波成分を含んでいるかの定性的な意味は、すでに微分フィルターが高周波抜き出しとして機能すると説明したところから類推するとよい。
二つ目のグラフでは、両端はつながっているが微分が連続ではないところが高周波としてカウントされる。低周波の波形は総じて振幅の割に一階微分が小さく、重ね合わせでこの波形を表現するには高周波な三角関数を使わないと表現できないからである(一つ目の理由も、じつはこっちで類推したほうがいいかもしれない)。
なお低周波成分が載っているのは、波形が上に偏っているから、が主な理由である(定数は周期0の波)。
窓関数の出番
さて、課題が見えたので窓関数の説明をする。
音声解析などで時間範囲を切って周波数解析をしたいが、適当に切るのは存在しない成分がのるし、周期を絞り込むことも難しい。適当に切るのがダメなのは不連続点と微分不連続性が理由。
では、適当に切った波形に対して不連続点と、微分不連続性が起きないように整形してからフーリエ変換してやればどうだろうか？この目的で用いるのが窓関数である。
窓関数のフーリエ変換応答
一定の範囲決めうちで信号の周波数を特定する際、端点の不連続性などに起因する偽信号を除去するには窓関数を利用する。
窓関数のかけかたは単純で、実空間で
$$g(x)=f(x)h(x)$$
と、成分ごとにかけ算するだけである。
これをフーリエ変換して得られる周波数は、さきほど導出したかけ算と畳み込みの関係により、
$$\mathcal F(g)=\mathcal F(f)\ast\mathcal F(h)$$
とわかっている。すなわち、窓関数後のフーリエ変換は元信号の周波数特性に対して、窓関数のフーリエ変換関数で周波数空間上でフィルター処理した結果となる。
この特性を理解するのは、フーリエ変換前後で関数が変わらないガウス関数を考えるとわかりやすい。
ガウス窓前後の応答
今回は効果が分かりやすいように、周期がわかっている関数に対してガウス窓をかけ、フーリエ変換器してみよう。
以下の図は、ガウス窓とフーリエ変換波形(上段)、元信号である三個の周波数の合成波形とそのフーリエ変換波形(中段)、ガウス窓の波形とフーリエ変換波形(下段)である。
上段の通り、ガウス関数は平滑化フィルタとして作用する特性をもつため、**中段の周波数スペクトルを平滑化処理したようなスペクトルが下段の周波数特性になる(周波数が負の値は正の値が折り返されているとみる)。
これは、実空間上でかけ算をしたことによる影響である。
このように、窓関数を使うと元の周波数の情報はほかの周波数の信号と混ざって出てくる。下段の図は元々三本の周波数があった波形が、周波数的に太い二本のピークで表現されている。つまり、一本は区別(分解)できていないことが起きる。
この、ピーク同士を区別できる性能のことを分解能とよぶ。
窓関数の種類については他にも色々ありwiki参照であるが、ここでは次にハン窓について語ろう。
ハン窓
窓サイズTのハン窓の考え方は単純なcosの基本振動数で定義される。
$$h(x)=\cos^2 \frac{\pi x}{T}=\frac{1+\cos\frac{2\pi x}{T}}{2}, \quad x\in[-T/2,T/2]$$
波形と周波数応答(および0付近の拡大図)は以下。
ここで注意したいのは、一定区間範囲のフーリエ変換で定義される周波数は離散的であるということである。冒頭の周期関数の定義を再掲しよう。
$$f(x)=\sum_{k=-\infty}^{\infty}\mathcal c_k\exp(i\frac{2\pi k}{T}x)$$
ただしハン窓は偶関数なので負の周波数成分は同じであるため省略している。位相成分(波形の開始位置のずれ)を考えないとき、実関数のフーリエ変換で得られた$c_k$の絶対値は正負で同じになるため、省略して表記している。
先ほどの波形にハン窓をかけてフーリエ変換しよう。
先ほどより分解能がいいように見えるが、これは半値幅のパラメーターの問題が支配的なので主な効果ではない。ハン窓のメリットは両端が微分も含めて0に落ちるため端点起因で高周波成分がのらないことにある。
この要件を満たし、単なる初等関数で表現でき、かつ半値幅も端の条件を満たす中では広い波形であることも人気の理由である。
もうひとつ重要な特性として、ハン窓による周波数空間上の畳み込みは、だいたい隣接周波数までしか影響しないことである。
厳密には計算範囲の倍数以外の周波数の波が信号のなかに入ってしまっている場合に、近似的に離散的な周波数のどれかにカウントされてしまうことがあるが、その影響は、計算範囲の長さを解析したい周波数の3倍から数倍程度の長さ確保していれば十分少ない。
波長の非整数倍の周期が紛れていた場合の結果の例を以下に示す。
低周波側の、周期の倍数でない周波数はそれを挟む二つの倍数となる周波数に補間する形でカウントされる。結果、実際の周波数とは若干ずれることはあるが、そのずれは最低限といえる。
このズレをなるべく小さくしたいならば窓の幅を大きくとるべきで、周波数解析の時間分解をよくしたいならば窓の幅を小さくし、代わりに周波数の分解能を諦める必要がある。
余談
量子力学の不確定性原理
量子力学の不確定性原理に時間とエネルギーの両者を同時に正確に求めることはできない、とある。エネルギーは波(物質波含む)の振動数に比例するため、これも今回の現象と類似した話である。少しだけ紹介しよう。
ある粒子の波動関数が
$$\phi=\exp(ikx-i\nu t)$$
であるとき、この運動量$p$とエネルギー$E$は
$$p=\hbar k$$
$$E=\hbar \nu$$
であり、不確定性原理は
$$\Delta x\Delta p,\Delta E\Delta t\geq \frac{\hbar}{2}$$
を満たす、というものである。(ここで、$\hbar=\frac{h}{2\pi}$はディラック定数といい、プランク定数を$2\pi$で割った値である)
これはつまり、
$$\Delta k\Delta x,\Delta \nu \Delta t\geq\frac{1}{2}$$
となる。
ある周波数$s$の波を
$$f(s)=\exp(2i\pi st)$$
としたとき、目標となる窓の長さを$T$とすると、周波数分解能$\Delta s$は
$$\Delta s\geq \frac{1}{4\pi T}$$
が限界である、と言っているものと解釈できる。
今のハン窓の手法ではここまでの分解能は出ないことは明らかなため、かなり頑張っている定義だと言える。
画像信号や多変量値のフーリエ変換
画像信号においても高周波、低周波は定義できる。画像は二次元の空間的に広がりを持つデータであるため、画像は座標ごとに明るさを定義した関数として扱うことができる。これを$f(x,y)$としてフーリエ変換係数で展開すると、
$$f(x,y)=\sum_{n_x,n_y}c_{k_x,k_y}\exp(2\pi i\lbrace\frac{k_xx}{T_x}+\frac{k_yy}{T_y}\rbrace)$$
とかける。これを変形すると
$$f(x,y)=\sum_{k} c_k'(y)\exp(2\pi i k\frac{x}{T_x})$$
$$c_k'(y)=\sum_{k_y}c_{k,k_y}\exp(2\pi ik_y\frac{y}{T_y})$$
とできる。この2行の式は各々がフーリエ変換係数による展開と考えることができるので、
以下の順序でフーリエ変換すれば$c_{k,k_y}$は求まる。
- 画像の行ごとに列方向の断面図をとり、得られた波形をフーリエ変換して$c_k'(y)$を求める
- $c_k'(y)$を$k$を固定してフーリエ変換する
逆変換は上記手順を逆に踏めば得られる。なお、式の対称性から行からフーリエ変換しても列からフーリエ変換しても同じ結果となる。
音声では時間を横軸に周波数を定義したが、画像では空間を横軸に周波数が定義される。これを空間周波数とよぶ。
空間周波数は$(x,y)$の二方向で定義できるので二次元ベクトルである。
例として、あるヒートマップをフーリエ変換したものを示そう。
一番上が0行目の断面図(横軸を列方向として.pixelの値をプロットした)であり、中段が入力したヒートマップである。
一番下がフーリエ変換したものであり、いくつかの数えられる程度の周波数ベクトルを離散的に重ね合わせたものであることがわかる。
より一般的には、任意の多変量関数$$$f(x_n)=f(x_1,x_2,x_3,...)$$
のフーリエ変換はn回の1次元フーリエ変換の繰り返しで定義できる。周波数ベクトルの次元の数は元の関数の次元の数に一致する。
$$F(w_1,w_2,...,w_n)=\mathcal F(f(x_1,x_2,...,x_n))$$
多変量の例として動画であれば$f(x,y,t)$、3DCG動画であれば$f(x,y,z,t)$として扱うことも可能である。$f(x,y,t)$をフーリエ変換したものを時空間周波数と呼ぶ・・・のかは知らないが、名称があるなら、きっとこの名前になるだろう。
なお、画像における(x,y)のことをpixelとよび、3DCGにおいて(x,y,z)のことをvoxelと呼ぶ。ついでに紹介しておく。
離散フーリエ変換の導出
画像や音声など、実際には連続ではなく離散的に定義されたデータについては離散フーリエ変換が便利である。
離散フーリエ変換は以下の式で表現される。
$$x_n=\sum_{k=0}^{N-1} c_k\exp(i\frac{2\pi k}{N}n)$$
このとき、
$$c_k=\frac{1}{N}\sum_{n=0}^{N-1}x_n\exp(-i\frac{2\pi n}{N}k)$$
であるが、この節で証明しよう。
$$\frac{1}{N}\sum_{n=0}^{N-1}x_n\exp(-i\frac{2\pi n}{N}k)=\frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}c_m\exp(2\pi i\frac{m-k}{N}n)=\frac{1}{N}\sum_{m=0}^{N-1}c_m\sum_{n=0}^{N-1}\exp(2\pi i\frac{m-k}{N}n)$$
と変形できる。ここで、
$$\sum_{n=0}^{N-1}\exp(2\pi i\frac{m-k}{N}n)=N\delta_{mk}$$
であることを証明できれば、導出完了である。
$$z=\exp(\frac{2\pi (m-k)i}{N})$$
とおくと、
$$z^N=1$$
である。等比級数の式
$$\sum_{n=0}^{N-1}z^n=\frac{1-z^n}{1-z}$$
を使えば、$z^N=1$は以下に因数分解可能だとわかる(勝手に分母を払って大丈夫なのか心配な方は、下記式を展開すれば検算可能である)。
$$(z-1)(\sum_{n=0}^{N-1}z^n)=0$$
因数分解は恒等変形なので、
$$m\neq k\Rightarrow z\neq 1\Rightarrow\sum_{n=0}^{N-1}z^n=0$$
$$m=q\Rightarrow z=1\Rightarrow \sum_{n=0}^{N-1}z^n=N$$
となる(ここで、$m,k\in[0,N-1]$)ことから、
$$\sum_{n=0}^{N-1}\exp(2\pi i\frac{m-k}{N}n)=\sum_{n=0}^{N-1}z^n=N\delta_{mk}$$
が成り立つ。これを冒頭の式に用いれば、
$$\frac{1}{N}\sum_{n=0}^{N-1}x_n\exp(-i\frac{2\pi kn}{N}k)=\frac{1}{N}\sum_{m=0}^{N-1}c_m\sum_{n=0}^{N-1}\exp(2\pi i\frac{m-k}{N}n)=\frac{1}{N}\sum_{m=0}^{N-1}c_m N\delta_{mk}=c_k$$
以上で証明完了である。
離散フーリエ変換のユニタリ性
前述の結果を用いると、
$$v_{kn}=\frac{1}{\sqrt{N}}\exp(\frac{2\pi k}{N}n)$$
で構成される行列はユニタリ行列である。
$$VV^\dagger=v_{kn}v_{nk'}^\ast=\frac{1}{N}\sum_n \exp(2\pi i\frac{k-k'}{N}n)=\delta_{kk'}=I$$
のため、
$$V^{-1}=V^\dagger$$
すなわち、Vはユニタリ行列である。
いま、
$$V_k=\lbrace v_{kn},n=0,1,2,..N-1\rbrace$$
というベクトルを定義したとき、その内積もまた
$$V_{k'}^TV_k^\ast=\delta_{kk'}$$
となる。つまり$V_k$は正規直交基底である。
フーリエ変換と離散フーリエ変換の対応関係
周期関数f(x)のフーリエ変換は、
$$F_k=\mathcal F(f)=\frac{1}{T}\int_{-T/2}^{T/2}f(x)\exp(\frac{2\pi i k}{T}x)dx$$
である。あるいは、$f(x)=f(x+T)$であり
$$F_k=\frac{1}{T}\int_0^T f(x)\exp(\frac{2\pi ik}{T}x)dx$$
としても等価である。これを用いると、
$$f(x)=\sum_{k=-\infty}^{\infty}F_k \exp(-\frac
{2\pi ik}{T}x)$$
とできるのであった。
これと離散フーリエ変換の対応関係をみてみよう。いま、実空間上で
$$h(x)=\sum_{n=-\infty}^{\infty} \delta(x-n)$$
という関数(くし型関数という)をかけ算した関数を定義する。
$$g(x)=f(x)h(x)$$
このカーネル関数と実空間上でかけ算をしてフーリエ変換しよう。周期$T$は整数値$N$に等しいとする。
$$G_k=\mathcal F(g)=\frac{1}{N}\int_{0}^{N}f(x)h(x)\exp(\frac{2\pi i k}{N}x)dx=\frac{1}{N}\left\lbrace\sum_{n=0}^{N}f(n)\exp(\frac{2\pi ik}{N}n)-\frac{1}{2}f(0)-\frac{1}{2}f(N)\right\rbrace$$
$$(w=\frac{2\pi k}{N})$$
なお、この展開には
$$\int_0^\infty \delta(x)f(x)dx=\frac {f(0)}{2}$$
$$\int_{-\infty}^{\infty}\delta(x)f(x)dx=f(0)$$
を利用する。
いま、関数の周期性から$f(0)=f(N)$であることに注意すると、
$$\mathcal F(g)=\frac{1}{N}\sum_{n=0}^{N-1}f(n)\exp(\frac{2\pi k}{N}n)$$
$$G_k=\mathcal F(g)=\frac{1}{N}\sum_{n=0}^{N-1}f(n)\exp(\frac{2\pi ik}{N}n)$$
と変形できる。ここで、$G_k$は$G_{k+N}=G_k$を満たす周期関数であることが一つ目の重要な知見である。
前の節で離散フーリエ変換を導出した通り、
$$f(n)=\sum_{k=0}^{N-1}G_k\exp(\frac{-2\pi ikn}{N})$$
と変形できるが、冒頭における
$$f(x)=\sum_{k=-\infty}^{\infty}F_k\exp(-\frac{2\pi ikx}{N})$$
を$x=n=0,1,2,...$として
$$f(n)=\sum_{k=-\infty}^{\infty}F_k\exp(-\frac{2\pi ikn}{N})$$
とし、
$$\exp(-\frac{2\pi i(k+mN))n}{N})=\exp(-\frac{2\pi ikn}{N}) \quad(\because n\in \mathbb Z)$$
を利用して変形すると、
$$f(n)=\sum_{k=0}^{N-1}\sum_{m=-\infty}^\infty F_{k+mN}\exp(-\frac{2\pi ikn}{N})$$
となる。$f(n)$を表現する式が二つでてきて、どの$n$に対しても成り立つことを満たすため、
$$0=\sum_{k=0}^{N-1}\left(\sum_{m=-\infty}^\infty F_{k+mN}-G_k\right)\exp(-\frac{2\pi ikn}{N})\quad { }^\forall n\in N$$
$$\therefore G_k=\sum_m F_{k+mN}$$
つまり、離散フーリエ変換の係数$G_k$は連続関数のフーリエ変換の係数$F_k$をサンプル数で切りそろえて足し合わせた値と等価である。
離散フーリエ変換の式は連続関数に対して櫛形関数を掛けてフーリエ変換した式と同じである。
これらを踏まえて、用語紹介もかねてまとめよう。
- 連続関数$f(x)$に対して離散的に$f(n)$という信号列を得る操作をサンプリングといい、信号列の長さ$N$をサンプリング長という。
- サンプリングしてから離散フーリエ変換して得られた数列(DFT応答)は元の連続関数に櫛形関数をかけてフーリエ変換した関数と等しく、元の連続関数に対してサンプリング長で切り揃えて足し合わせた数列に等しい。定義上、DFT応答はサンプリング長の周期で同じ値が繰り返される。
- DFT応答で最も高周波な応答は$N/2$であり、これをナイキスト周波数とよぶ。それよりも高い$k$は低周波の応答の一部となる。これを周波数の折り返しとよぶ。
周波数の折り返しは、以下の数式で簡単に表現できる。
$$G_k=G_{k+N}$$
$$G_{N/2+k}=G_{N/2+k-N}=G_{-(N/2-k)}$$
フーリエ変換において$c_k$,$c_{-k}$は同じ周波数(共役な振動)の係数であるから、$N/2$を挟んで等距離にある周波数は同じ振動数の成分であるといえる。
なお、サンプリングした数列が実数の数列であった場合、折り返す前と後の応答は共役関係になる。
$$G_{N/2+k}=G_{N/2-k}^\ast$$
ナイキスト周波数とモアレ
連続関数で定義した波形について、離散的に信号を取り出す操作をサンプリングと呼ぶことは前の節で紹介した。
サンプリングの操作は元の信号の一部を取り出す操作であることより、元の情報を失う操作だとおえるだろう。今回、その落ちる度合いについて少し議論しよう。
例として、周期8の波形から離散的に整数画素をサンプリングする操作を考える。
すると、元の波形の位相によってサンプリングする振幅が変わる。
各行は元波形の周期を変えながらプロットした様子を示している。左の列における青線は元の連続関数の波形、赤マーカーが整数値でサンプリングした点、緑マーカーが半画素ずらしてプロットした例である。右の図は赤マーカー、緑マーカーをそれぞれ線で結んだグラフである。
上から三段目は元波形の周期がちょうど2である場合を指している。言い換えれば、2サンプルで1周期となる周波数である。このとき、赤プロットは全く振幅のない波形として見なされ、一方半画像だけずらしてプロットすると振幅が最大となる。このような周波数をナイキスト周波数といい、この付近が元の波形信号を表現できる限界だといえる。それ以上に短い周波数を元波形とするとサンプリング波形は低周波な波形として現れることがみえるが、これが先述の周波数折り返しが現れた例であり、本来存在しないはずの周波数が見えてしまったことになる。これをモアレという。