POINT
- 離散フーリエ変換(DFT)に関するまとめ.
- 計算機では有限個の離散データしか扱えない.そこで,有限個の離散データを周期的に拡張して扱う.
- 非周期的なデータは扱えず,周期的なデータとなる.
- フーリエ変換は離散フーリエ変換(DFT)として扱い,高速フーリエ変換(FFT)というアルゴリズムで計算される.
フーリエ変換と離散フーリエ変換(DFT)の関係
フーリエ変換は連続変数の関数を,無限の積分区間で計算する必要がある.しかし,現実的には離散時間の有限個のデータをサンプリングすることしかできない.サンプリングで得られた離散的で有限個のデータ列から,フーリエ変換(スペクトル)を推定する方法が離散フーリエ変換(DFT)である.どうしてDFTでフーリエ変換が推定できるか,ということについては次の記事を参照.
フーリエ変換とDFTのつながり - Notes_JP
DFTは高速フーリエ変換(FFT)というアルゴリズムで計算することが多い.以下の記事では,三角関数のDFTを「手計算」し,FFTの結果と比較した.
正弦波のFFT (numpy.fft) - Notes_JP
離散フーリエ変換(DFT)
$N$個の離散データ列$\{g(k) \,|\, k=0,1,,...,N-1\}$が与えられたとき,$g$を$g(k+N)=g(k)$を満たす周期$N$の関数と考える.このとき,$g(k)$の離散フーリエ変換(DFT)は離散フーリエ変換($N$点DFT)
\begin{aligned}
G(n) = \sum_{k=0}^{N-1} g(k) e^{-i 2\pi nk/N}
\end{aligned}
G(n) = \sum_{k=0}^{N-1} g(k) e^{-i 2\pi nk/N}
\end{aligned}
このとき,
\begin{aligned}
G(n+N)
&= \sum_{k=0}^{N-1} g(k) e^{-i (2\pi nk/N + 2\pi k)} \\
&=G(n)
\end{aligned}
なので,$G$も周期$N$の関数である.G(n+N)
&= \sum_{k=0}^{N-1} g(k) e^{-i (2\pi nk/N + 2\pi k)} \\
&=G(n)
\end{aligned}
また,
\begin{aligned}
G(0)
&=\sum_{k=0}^{N-1} g(k)
\end{aligned}
なので,スペクトルの直流成分は時間データの和になる.G(0)
&=\sum_{k=0}^{N-1} g(k)
\end{aligned}
逆離散フーリエ変換(IDFT)
逆変換は次で与えられる.$g=\mathcal{F}^{-1}[G]$と表す.逆離散フーリエ変換($N$点IDFT)
\begin{aligned}
g(k) = \frac{1}{N} \sum_{n=0}^{N-1} G(n) e^{i 2\pi nk/N}
\end{aligned}
g(k) = \frac{1}{N} \sum_{n=0}^{N-1} G(n) e^{i 2\pi nk/N}
\end{aligned}
DFTの両辺に$e^{i 2\pi nl/N}$をかけて$\displaystyle \sum_{n=0}^{N-1}$とすると
\begin{aligned}
&\sum_{n=0}^{N-1} G(n) e^{i 2\pi nl/N} \\
&=\sum_{k=0}^{N-1}
\Biggl[g(k) \cdot \sum_{n=0}^{N-1} e^{-i 2\pi n(k-l) /N} \Biggr]
\end{aligned}
となる.ここで,&\sum_{n=0}^{N-1} G(n) e^{i 2\pi nl/N} \\
&=\sum_{k=0}^{N-1}
\Biggl[g(k) \cdot \sum_{n=0}^{N-1} e^{-i 2\pi n(k-l) /N} \Biggr]
\end{aligned}
\begin{aligned}
&\sum_{n=0}^{N-1} e^{-i 2\pi n(k-l) /N} \\
&=
\begin{cases}
\displaystyle \frac{1-e^{-i 2\pi (k-l)}}{1-e^{-i 2\pi (k-l) /N}}=0 & (k\neq l) \\
\, N & (k = l)
\end{cases} \\
&=N\delta_{k,l}
\end{aligned}
であるから&\sum_{n=0}^{N-1} e^{-i 2\pi n(k-l) /N} \\
&=
\begin{cases}
\displaystyle \frac{1-e^{-i 2\pi (k-l)}}{1-e^{-i 2\pi (k-l) /N}}=0 & (k\neq l) \\
\, N & (k = l)
\end{cases} \\
&=N\delta_{k,l}
\end{aligned}
\begin{aligned}
&\sum_{n=0}^{N-1} G(n) e^{i 2\pi nl/N} \\
&= \sum_{k=0}^{N-1} g(k) \cdot N\delta_{k,l} \\
&=N g(l)
\end{aligned}
が成り立つ.//&\sum_{n=0}^{N-1} G(n) e^{i 2\pi nl/N} \\
&= \sum_{k=0}^{N-1} g(k) \cdot N\delta_{k,l} \\
&=N g(l)
\end{aligned}
ここで,
\begin{aligned}
g(0) = \frac{1}{N} \sum_{n=0}^{N-1} G(n)
\end{aligned}
なので,スペクトルを平均すると時刻ゼロのデータが得られる.g(0) = \frac{1}{N} \sum_{n=0}^{N-1} G(n)
\end{aligned}
性質
演算時の注意
$N$点DFTの演算を行う際は,演算に使う2つの信号の長さを$N$に揃える必要がある.長さが異なる場合は,短い方にゼロを加えて同じ長さにする.線形性
線形性
$G_{1}=\mathcal{F}[g_{1}], G_{2}=\mathcal{F}[g_{2}]$と定数$a_{1}, a_{2}$に対して
【証明】\begin{aligned}
\mathcal{F}[a_{1}g_{1} + a_{2}g_{2}]
=a_{1}G_{1} + a_{2}G_{2}
\end{aligned}
\mathcal{F}[a_{1}g_{1} + a_{2}g_{2}]
=a_{1}G_{1} + a_{2}G_{2}
\end{aligned}
定義から簡単に示せる.//
対称性・反対称性
$G=\mathcal{F}[g]$であるとする.また,$A$の複素共役を$A^{*}$で表す.対称性
\begin{aligned}
G^{*}(n)
&=
\begin{cases}
\, G(N-n) & (g\text{が実関数}) \\
\, -G(N-n) & (g(k)\text{が純虚数})
\end{cases}
\end{aligned}
が成り立つ.G^{*}(n)
&=
\begin{cases}
\, G(N-n) & (g\text{が実関数}) \\
\, -G(N-n) & (g(k)\text{が純虚数})
\end{cases}
\end{aligned}
注:前者を「周波数スペクトルが$n=0$に対して共役対称 (Hermitian symmetry) である」と表現し,後者を「周波数スペクトルが$n=0$に対して反共役対称 (skew Hermitian symmetry) である」と表現する.
【証明】
DFTの複素共役をとると
\begin{aligned}
G^{*}(n)
&=\pm \sum_{k=0}^{N-1} g(k) e^{i 2\pi nk/N} \\
&=\pm \sum_{k=0}^{N-1} g(k) e^{-i 2\pi (N-n)k/N} \\
&=\pm G(N-n)
\end{aligned}
となる($g(k)$が実数の場合と,純虚数の場合を併記した).//G^{*}(n)
&=\pm \sum_{k=0}^{N-1} g(k) e^{i 2\pi nk/N} \\
&=\pm \sum_{k=0}^{N-1} g(k) e^{-i 2\pi (N-n)k/N} \\
&=\pm G(N-n)
\end{aligned}
巡回時間シフト
巡回周波数シフト
$G=\mathcal{F}[g]$であるとき
【証明】\begin{aligned}
\mathcal{F} [g(k-k_{0})](n)
=G(n) e^{-i 2\pi nk_{0}/N}
\end{aligned}
注:$\mathcal{F}[g\textcolor{red}{(k)}]$のような書き方は良くないですが,公式を簡単に表現でき,意味も伝わるため,採用しました.\mathcal{F} [g(k-k_{0})](n)
=G(n) e^{-i 2\pi nk_{0}/N}
\end{aligned}
\begin{aligned}
&\mathcal{F} [g(k-k_{0})](n) \\
&= \sum_{k=0}^{N-1} g(k-k_{0}) e^{-i 2\pi nk/N} \\
&= \underbrace{\sum_{k=0}^{N-1} g(k-k_{0}) e^{-i 2\pi n(k-k_{0})/N}}_{=\sum_{l=0}^{N-1}g(l)e^{-i 2\pi nl/N}=G(n)} \cdot e^{-i 2\pi nk_{0}/N} \\
&=G(n)e^{-i 2\pi nk_{0}/N}
\end{aligned}
&\mathcal{F} [g(k-k_{0})](n) \\
&= \sum_{k=0}^{N-1} g(k-k_{0}) e^{-i 2\pi nk/N} \\
&= \underbrace{\sum_{k=0}^{N-1} g(k-k_{0}) e^{-i 2\pi n(k-k_{0})/N}}_{=\sum_{l=0}^{N-1}g(l)e^{-i 2\pi nl/N}=G(n)} \cdot e^{-i 2\pi nk_{0}/N} \\
&=G(n)e^{-i 2\pi nk_{0}/N}
\end{aligned}
巡回周波数シフト
巡回周波数シフト
$G=\mathcal{F}[g]$であるとき
【証明】\begin{aligned}
\mathcal{F}^{-1}[G(n-n_{0})](k)
&= g(k) e^{i 2\pi n_{0}k/N}
\end{aligned}
注:$\mathcal{F}^{-1}[G\textcolor{red}{(n-n_{0})}]$のような書き方は良くないですが,公式を簡単に表現でき,意味も伝わるため,採用しました.\mathcal{F}^{-1}[G(n-n_{0})](k)
&= g(k) e^{i 2\pi n_{0}k/N}
\end{aligned}
\begin{aligned}
&\mathcal{F}^{-1}[G(n-n_{0})](k) \\
&=\frac{1}{N} \sum_{n=0}^{N-1} G(n-n_{0}) e^{i 2\pi n k/N} \\
&=\underbrace{\frac{1}{N} \sum_{n=0}^{N-1} G(n-n_{0}) e^{i 2\pi (n-n_{0}) k/N}}_{=\frac{1}{N} \sum_{m=0}^{N-1} G(m) e^{i 2\pi m k/N}=g(k)} \cdot e^{i 2\pi n_{0} k/N} \\
&=g(k) e^{i 2\pi n_{0}k/N}
\end{aligned}
&\mathcal{F}^{-1}[G(n-n_{0})](k) \\
&=\frac{1}{N} \sum_{n=0}^{N-1} G(n-n_{0}) e^{i 2\pi n k/N} \\
&=\underbrace{\frac{1}{N} \sum_{n=0}^{N-1} G(n-n_{0}) e^{i 2\pi (n-n_{0}) k/N}}_{=\frac{1}{N} \sum_{m=0}^{N-1} G(m) e^{i 2\pi m k/N}=g(k)} \cdot e^{i 2\pi n_{0} k/N} \\
&=g(k) e^{i 2\pi n_{0}k/N}
\end{aligned}
巡回畳み込み
$g_{1}$と$g_{2}$の巡回畳み込みは巡回畳み込み
\begin{aligned}
(g_{1}\otimes g_{2})(\textcolor{red}{k})
&=\sum_{i=0}^{N-1} g_{1}(i) g_{2}(\textcolor{red}{k}-i) \\
&=\sum_{\substack{i+j=\textcolor{red}{k} \\ 0 \leq i,j \leq N-1}} g_{1}(i) g_{2}(j)
\end{aligned}
(g_{1}\otimes g_{2})(\textcolor{red}{k})
&=\sum_{i=0}^{N-1} g_{1}(i) g_{2}(\textcolor{red}{k}-i) \\
&=\sum_{\substack{i+j=\textcolor{red}{k} \\ 0 \leq i,j \leq N-1}} g_{1}(i) g_{2}(j)
\end{aligned}
このとき,次が成り立つ.
巡回畳み込みのDFT
$G_{1}=\mathcal{F}[g_{1}], G_{2}=\mathcal{F}[g_{2}]$とすると,
【証明】\begin{aligned}
\mathcal{F}[g_{1}\otimes g_{2}] (n)
&=G_{1}(n)\cdot G_{2}(n)
\end{aligned}
が成り立つ.\mathcal{F}[g_{1}\otimes g_{2}] (n)
&=G_{1}(n)\cdot G_{2}(n)
\end{aligned}
\begin{aligned}
&\mathcal{F}[g_{1}\otimes g_{2}] (n) \\
&= \sum_{k=0}^{N-1} \Biggl[ \sum_{i=0}^{N-1} g_{1}(i) g_{2}(k-i) \Biggr] e^{-i 2\pi nk/N} \\
&=\sum_{i=0}^{N-1} g_{1}(i)e^{-i 2\pi n i/N}
\underbrace{\sum_{k=0}^{N-1} g_{2}(k-i) e^{-i 2\pi n(k-i)/N}}_{=\sum_{l=0}^{N-1} g_{2}(l) e^{-i 2\pi nl/N}=G_{2}(n)} \\
&=G_{1}(n)\cdot G_{2}(n)
\end{aligned}
&\mathcal{F}[g_{1}\otimes g_{2}] (n) \\
&= \sum_{k=0}^{N-1} \Biggl[ \sum_{i=0}^{N-1} g_{1}(i) g_{2}(k-i) \Biggr] e^{-i 2\pi nk/N} \\
&=\sum_{i=0}^{N-1} g_{1}(i)e^{-i 2\pi n i/N}
\underbrace{\sum_{k=0}^{N-1} g_{2}(k-i) e^{-i 2\pi n(k-i)/N}}_{=\sum_{l=0}^{N-1} g_{2}(l) e^{-i 2\pi nl/N}=G_{2}(n)} \\
&=G_{1}(n)\cdot G_{2}(n)
\end{aligned}
積
$1/N$がつくことに注意.積のDFT
$G_{1}=\mathcal{F}[g_{1}], G_{2}=\mathcal{F}[g_{2}]$とすると,
【証明】\begin{aligned}
\mathcal{F}[g_{1}\cdot g_{2}] (n)
&=\frac{1}{N} (G_{1}\otimes G_{2}) (n)
\end{aligned}
が成り立つ.\mathcal{F}[g_{1}\cdot g_{2}] (n)
&=\frac{1}{N} (G_{1}\otimes G_{2}) (n)
\end{aligned}
\begin{aligned}
&\mathcal{F}[g_{1}\cdot g_{2}] (n) \\
&=\mathcal{F}[\mathcal{F}^{-1}(G_{1})\cdot g_{2}] (n) \\
&= \sum_{k=0}^{N-1} \Biggl[\frac{1}{N} \sum_{m=0}^{N-1} G_{1}(m)e^{i 2\pi mk/N} \cdot g_{2}(k)\Biggr] e^{-i 2\pi nk/N} \\
&=\frac{1}{N}\sum_{m=0}^{N-1} G_{1}(m) \underbrace{\sum_{k=0}^{N-1} g_{2}(k) e^{-i 2\pi (n-m) k/N}}_{G(n-m)} \\
&=\frac{1}{N} (G_{1}\otimes G_{2}) (n)
\end{aligned}
&\mathcal{F}[g_{1}\cdot g_{2}] (n) \\
&=\mathcal{F}[\mathcal{F}^{-1}(G_{1})\cdot g_{2}] (n) \\
&= \sum_{k=0}^{N-1} \Biggl[\frac{1}{N} \sum_{m=0}^{N-1} G_{1}(m)e^{i 2\pi mk/N} \cdot g_{2}(k)\Biggr] e^{-i 2\pi nk/N} \\
&=\frac{1}{N}\sum_{m=0}^{N-1} G_{1}(m) \underbrace{\sum_{k=0}^{N-1} g_{2}(k) e^{-i 2\pi (n-m) k/N}}_{G(n-m)} \\
&=\frac{1}{N} (G_{1}\otimes G_{2}) (n)
\end{aligned}
パーシバルの定理
パーシバルの定理
$G=\mathcal{F}[g]$であるとき
【証明】\begin{aligned}
\sum_{k=0}^{N-1} |g(k)|^{2}
&=\frac{1}{N}\sum_{n=0}^{N-1} |G(n)|^{2}
\end{aligned}
が成り立つ.\sum_{k=0}^{N-1} |g(k)|^{2}
&=\frac{1}{N}\sum_{n=0}^{N-1} |G(n)|^{2}
\end{aligned}
積のDFT公式から
\begin{aligned}
&\sum_{k=0}^{N-1} g_{1}(k)\cdot g_{2}(k) e^{-i 2\pi nk/N} \\
&=\frac{1}{N}\sum_{m=0}^{N-1} G_{1}(m)G_{2}(n-m)
\end{aligned}
である.この式で,$g_{1}=g, g_{2}=g^{*},n=0$とすれば&\sum_{k=0}^{N-1} g_{1}(k)\cdot g_{2}(k) e^{-i 2\pi nk/N} \\
&=\frac{1}{N}\sum_{m=0}^{N-1} G_{1}(m)G_{2}(n-m)
\end{aligned}
\begin{aligned}
\sum_{k=0}^{N-1} |g(k)|^{2}
&=\frac{1}{N}\sum_{m=0}^{N-1} G(m)
\cdot \underbrace{\sum_{k=0}^{N-1}g^{*}(k)e^{i 2\pi mk/N}}_{\mathrlap{=G^{*}(m)}} \\
&=\frac{1}{N}\sum_{n=0}^{N-1} |G(n)|^{2}
\end{aligned}
がわかる.//\sum_{k=0}^{N-1} |g(k)|^{2}
&=\frac{1}{N}\sum_{m=0}^{N-1} G(m)
\cdot \underbrace{\sum_{k=0}^{N-1}g^{*}(k)e^{i 2\pi mk/N}}_{\mathrlap{=G^{*}(m)}} \\
&=\frac{1}{N}\sum_{n=0}^{N-1} |G(n)|^{2}
\end{aligned}
畳み込みをプログラムで計算するには
$0\leq k \leq l_{1}-1$でのみゼロでない値をとるデータ列$\{g_{1}(k)\}_{k=-\infty}^{\infty}$と,$0\leq k \leq l_{2}-1$でのみゼロでない値をとるデータ列$\{g_{2}(k)\}_{k=-\infty}^{\infty}$が与えられたとする.このとき,これらの線形畳み込み演算\begin{aligned}
y(k)
&=\sum_{i=-\infty}^{\infty} g_{1}(i) g_{2}(k-i) \\
&=\sum_{i=0}^{\textcolor{red}{k}} g_{1}(i) g_{2}(k-i)
\end{aligned}
をコンピュータで計算する方法を考える.y(k)
&=\sum_{i=-\infty}^{\infty} g_{1}(i) g_{2}(k-i) \\
&=\sum_{i=0}^{\textcolor{red}{k}} g_{1}(i) g_{2}(k-i)
\end{aligned}
コンピュータで利用できるのは$N$点DFTであるから,
- $N (\geq l_{1}, l_{2})$を適切に定め,
- ($N$周期)データ列$\{g_{1}(k)\}_{k=0}^{N-1}$,$\{g_{2}(k)\}_{k=0}^{N-1}$をつくり,
- 巡回畳み込み\begin{aligned}によって,上の線形畳み込みを計算する
g_{1}\otimes g_{2}(k)
&=\sum_{i=1}^{N-1} g_{1}(i) g_{2}(k-i)
\end{aligned}
ここで,線形畳み込み
\begin{aligned}
y(k)
&=\sum_{i=0}^{k} g_{1}(i) g_{2}(k-i)
\end{aligned}
がゼロでない値を取り得るのは,y(k)
&=\sum_{i=0}^{k} g_{1}(i) g_{2}(k-i)
\end{aligned}
\begin{aligned}
\begin{cases}
\, 0\leq i \leq l_{1}-1 \\
\, 0\leq k -i \leq l_{2}-1
\end{cases}
\end{aligned}
が両立する範囲だから,$0 \leq k \leq (l_{1}-1) + (l_{2} - 1)$となる.\begin{cases}
\, 0\leq i \leq l_{1}-1 \\
\, 0\leq k -i \leq l_{2}-1
\end{cases}
\end{aligned}
以上より,
まとめ
データ列$\{g_{1}(k)\}_{k=0}^{l_{1}-1}, \{g_{2}(k)\}_{k=0}^{l_{2}-1}$に対し,$N \geq l_{1} + l_{2} - 1$を満たす$N$をとれば,巡回畳み込みの式
\begin{aligned}
g_{1}\otimes g_{2}(k)
&=\sum_{i=1}^{N-1} g_{1}(i) g_{2}(k-i)
\end{aligned}
は$0 \leq k \leq l_{1} + l_{2} - 2$の範囲で線形畳み込みg_{1}\otimes g_{2}(k)
&=\sum_{i=1}^{N-1} g_{1}(i) g_{2}(k-i)
\end{aligned}
\begin{aligned}
y(k)
&=\sum_{i=-\infty}^{\infty} g_{1}(i) g_{2}(k-i) \\
&=\sum_{i=0}^{\textcolor{red}{k}} g_{1}(i) g_{2}(k-i)
\end{aligned}
と一致する.y(k)
&=\sum_{i=-\infty}^{\infty} g_{1}(i) g_{2}(k-i) \\
&=\sum_{i=0}^{\textcolor{red}{k}} g_{1}(i) g_{2}(k-i)
\end{aligned}
注:
- 与えられていない添字に関するデータは,線形畳み込みでは$0$として定義し,巡回畳み込みでは周期的に定義するものとする.
- FFTでは$N=2^{n}$になるように選ぶ.
参考文献
- [1]信号とシステム (斉藤洋一):7. 離散フーリエ変換
(この記事を作ったきっかけ:F - Substring 2 | AtCoder Beginner Contest 196)