修正離散コサイン変換をFFTで表して高速化するには
この記事は、このブログで以前に書いた記事に追加して書こうと思いましたが
記事の容量限界を超えたらしくエラーがでたので、新規記事として書きます。
まずは、修正離散コサイン変換をFFTで表す。から
\begin{equation} S(r)=\sum_{k=0}^{N-1} s(k){\rm cos}(\frac {2 \pi} {N} (k+ \frac{1}{2} + \frac{N}{4})(r+\frac{1}{2})) ただし r=0 ～ \frac{N}{2}-1 \end{equation}
\begin{equation} = \frac{1}{2} \sum_{k=N/4+1/2}^{3N/4-1/2}{ \rm exp}(i \frac{ \pi k}{N})(s(k- \frac{N}{4}- \frac{1}{2})-s( \frac{3N}{4}- \frac{1}{2}-k)){\rm exp}(i \frac {2 \pi}{N}kr ) \end{equation}
\begin{equation} + \frac{1}{2} \sum_{k=3N/4+1/2}^{N-1/2}{ \rm exp}(i \frac{ \pi k}{N})(s(k- \frac{N}{4}- \frac{1}{2})+s( \frac{7N}{4}- \frac{1}{2}-k)){ \rm exp}(i \frac {2 \pi}{N}kr) \end{equation}
\begin{equation} - \frac{1}{2} \sum_{k=1/2}^{N/4-1/2}{ \rm exp}(i \frac{ \pi k}{N})(s(k+ \frac{3N}{4}- \frac{1}{2})+s( \frac{3N}{4}- \frac{1}{2}-k)){ \rm exp}(i \frac {2 \pi}{N}kr) \end{equation}
ここまでは、前の記事の結果です。
さて、しかし、これではまだ不完全で、フーリエ変換の範囲が $ k=1/2～N-1/2 $
になっているので、 $ k=0～N-1 $ になおします。
\begin{equation} S(r)= \frac{1}{2} \sum_{k=N/4}^{3N/4-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) (s(k- \frac{N}{4})-s( \frac{3N}{4}-1-k)) exp(i \frac{2π}{N} (k+ \frac{1}{2})r) \end{equation}
\begin{equation} + \frac{1}{2} \sum_{k=3N/4}^{N-1} exp(iπ \frac{k+ \frac{1}{2}}{N})(s(k- \frac{N}{4})+s( \frac{7N}{4}-1-k)) exp(i \frac{2π}{N} (k+ \frac{1}{2})r) \end{equation}
\begin{equation} - \frac{1}{2} \sum_{k=0}^{N/4-1} exp(iπ \frac{k+ \frac{1}{2}}{N})(s(k+ \frac{3N}{4})+s( \frac{3N}{4}-1-k)) exp(i \frac{2π}{N} (k+ \frac{1}{2})r) \end{equation}
となります。また、フーリエ変換のかたちではなくなったので、 $ k+ \frac{1}{2} $ の
$ \frac{1}{2} $ を外に出します。この項は、 $ k $ を含まないので $ \sum $ の前に
だせます。
\begin{equation} S(r)= \frac{1}{2} exp(i \frac{πr}{N}) \sum_{k=N/4}^{3N/4-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) (s(k- \frac{N}{4})-s( \frac{3N}{4}-1-k)) exp(i \frac{2π}{N} kr) \end{equation}
\begin{equation} + \frac{1}{2} exp(i \frac{πr}{N}) \sum_{k=3N/4}^{N-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) (s(k- \frac{N}{4})+s( \frac{7N}{4}-1-k)) exp(i \frac{2π}{N} kr) \end{equation}
\begin{equation} - \frac{1}{2} exp(i \frac{πr}{N}) \sum_{k=0}^{N/4-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) (s(k+ \frac{3N}{4})+s( \frac{3N}{4}-1-k)) exp(i \frac{2π}{N} kr) \end{equation}
これで終了ですが、慣れてない方にさらに説明すると、上の式の3行は、実は次の
ように1行にまとめられます。
\begin{equation} S(r)= \frac{1}{2} exp(i \frac{πr}{N}) \sum_{k=0}^{N-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) p(k) exp(i \frac{2π}{N} kr) \end{equation}
$ p(k) $ は、$ \sum $の範囲でそれぞれ定義される関数をつなげて1個にした関数だから
です。この $ p(k) $ は、前半の奇対称数列部と後半の偶対称数列部をつなげたものを
4等分して、おしりの1/4を左にもってきて、左にもってきたものだけ符号をかえた
ものになります。奇対称数列になります。
フーリエ変換対象関数は、これを $ 0～N-1 $の区間で $ Re-Im-t $ 空間
で螺旋状に半回転させたものになります。
また $ r=0～N-1 $ で値がもとまりますが、実際は半分の $ r=0～N/2-1 $
しか計算しません。
奇関数のフーリエ変換はサイン変換になり半分の計算ですむからです。
残りの半分計算しても、同じものを計算することになるからです。
さて、逆変換について説明します。
$ p(k) $ は、奇対称数列でした。前半奇対称数列、後半偶対称数列になる逆変換結果
を求めるには、$ p(k) $ を $ S(r) $ から求めて、数列を4等分して、一番前の $ 1/4 $
を一番後ろの後に符号をかえてくっつければ、完成となります。
そこで $ p(k) $ を $ S(r) $ の逆フーリエ変換であらわせばよいことになります。
ここで $ S(r) $ は、修正離散コサイン変換結果で $ r=0～N/2-1 $ の範囲の数列でした。
しかし、実際は、 $ S(r) $ を奇関数展開して、$ r=0～N-1 $ の範囲にしたものでないと
逆フーリエ変換できません。そこで、奇関数展開した $ S(r) $ を $ SS(r) $ とおくと
\begin{equation} SS(r)= \frac{1}{2} exp(i \frac{πr}{N}) \sum_{k=0}^{N-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) p(k) exp(i \frac{2π}{N} kr) \end{equation}
\begin{equation} 2 exp(-i \frac{πr}{N})SS(r) = \sum_{k=0}^{N-1} exp(iπ \frac{k+ \frac{1}{2}}{N}) p(k) exp(i \frac{2π}{N} kr) \end{equation}
\begin{equation} \sum_{r=0}^{N-1} 2 exp(-i \frac{πr}{N})SS(r)exp(-i \frac{2π}{N} kr) =exp(iπ \frac{k+ \frac{1}{2}}{N}) p(k) \end{equation}
\begin{equation} p(k) = exp(-iπ \frac{k+ \frac{1}{2}}{N}) \sum_{r=0}^{N-1} 2 exp(-i \frac{πr}{N})SS(r)exp(-i \frac{2π}{N} kr) \end{equation}
となり $ p(k) $ が求まります。もとまった $ p(k) $を先ほど述べた操作を加えて
前半奇対称数列、後半偶対称数列の逆変換結果が得られるわけです。
C言語でのプログラム例のソースは
zuruyasumineko2002.blog.fc2.com
に載せています。
以上です。
