はじめに
以前、離散サイン変換（DST-I）と fft との関係は確認しましたが、DST-II と DST-III のタイプの離散サイン変換も必要になったので、改めて fft 関数を使って計算する方法を確認してみました。
何通りもやり方がありますし各所にコードも転がっている（例えば File Exchange: Wavelet Scale Spectra ）んですが、やっぱり確認しておかないと・・ということで備忘録代わりに。
ちなみに DST-I については MATLAB の Partial Differential Equation Toolbox に dst という関数があります。
定義式
まずは定義式から。長さ $N$の信号 $x$ に対して離散フーリエ変換を $Y(k)$、離散コサイン変換を $y(k)$ とおきます。$\delta_{kN}$ はクロネッカーデルタです。
\begin{array}{lcl}
{{\rm DST-I}}: & y(k) & =\sqrt{\frac{2}{N+1}}\sum_{j=1}^N x(j)\sin \left(\frac{\pi }{N+1}jk\right),\\
{{\rm DST-II}}: & y(k) & =\sqrt{\frac{2}{N}}\sum_{j=1}^N x(j)\frac{1}{\sqrt{1+\delta_{kN} }}\sin \left(\frac{\pi }{2N}(2j-1)k\right),\\
{{\rm DST-III}}: & y(k) & =\sqrt{\frac{2}{N}}\sum_{j=1}^N x(j)\frac{1}{\sqrt{1+\delta_{jN} }}\sin \left(\frac{\pi }{2N}j(2k-1)\right),\\
{{\rm FFT}}: & Y(k) & =\sum_{j=1}^N x(j)\exp \left(\frac{-2\pi i}{N}(j-1)(k-1)\right).
\end{array}
注意: ご存知の通り MATLABのベクトルは $0$ から $N-1$ではなく $1$ から $N$で定義されるため、 和は $1$ から $N$ としています。コピペするときは要注意（笑）
ちなみに FFTW に掲載されている定義式は以下の通り。
\begin{array}{lrl}
{{\rm DST-I}}: & y(k) & =2\sum_{j=0}^{N-1} x(j)\sin \left(\frac{\pi }{N+1}(j+1)(k+1)\right).\\
{{\rm DST-II}}: & y(k) & =2\sum_{j=0}^{N-1} x(j)\sin \left(\frac{\pi }{N}(j+1/2)(k+1)\right).\\
{{\rm DST-III}}: & y(k) & =(-1)^k X_{n-1} +2\sum_{j=0}^{N-2} x(j)\sin \left(\frac{\pi }{N}(j+1)(k+1/2)\right).
\end{array}
DST-I の計算
これは既に 離散フーリエ変換から離散サイン変換までの式展開 で一例を記載していますが、式展開がもう少し楽になる方法を改めて記載します。
変換対象の元データを x （サイズ N )とします。
xx=[0,x(1),x(2),x(3),\cdots ,x(N),0,0,0,\cdots ,0,0]
と拡張した長さ $2N+2$ のベクトル $xx$ に対しての離散フーリエ変換から式を展開していきます。
上の FFT の定義式から
Y(k)=\sum_{j=1}^{2N+2} xx(j)\exp \left(\frac{-2\pi i}{2N+2}(j-1)(k-1)\right)=\sum_{j=2}^{N+1} xx(j)\exp \left(\frac{-2\pi i}{2N+2}(j-1)(k-1)\right)
x(0) = 0, j= N+1 以降はすべて 0 であることを反映しました。$k$ と $j$ について微調整すると DST-I の定義式に近づきます。
Y(k+1)=\sum_{j=1}^N xx(j+1)\exp \left(\frac{-\pi i}{N+1}jk\right)=\sum_{j=1}^N x(j)\exp \left(\frac{-\pi i}{N+1}jk\right)
虚部を取ると
\begin{array}{cl}
{{{\rm Im}(Y(k+1))}} & =-\sum_{j=1}^N xx(j+1)\sin \left(\frac{\pi }{N+1}jk\right)
\end{array}
DST-I の定義
y(k)=\sqrt{\frac{2}{N+1}}\sum_{j=1}^N x(j)\sin \left(\frac{\pi }{N+1}jk\right),
と比べて、最終的に以下の通り DST-I と DFT の関係が求まります。
y(k)=-{{\rm Im}}(Y(k+1))\sqrt{\frac{2}{N+1}}.
拡張した長さ $2N+2$ のベクトル $xx$ に対しての離散フーリエ変換 $Y(k)$ の $k=2$ から $k=N+1$ を下の式で変換したものが、求めたい DST の係数です。
DST-II の計算
ここでは $x(1)$から $x(N)$ を以下のように拡張します。
変換対象の元データを x （サイズ N )とします。
xx=[x(1),x(2),x(3),\cdots ,x(N),0,0,0,\cdots ,0,0]
全長 $2N$ です。長さ $2N$ の離散フーリエ変換を考えて
w(k)=\exp \left(\frac{-\pi i}{2N}(k-1)\right),
という重みをかけた式を展開していきます。半サンプル分データをシフトすることに対応する処理と理解しています。上の定義式から
Y(k)w(k)=\sum_{j=1}^{2N} xx(j)\exp \left(\frac{-\pi i}{2N}(2j-1)(k-1)\right)=\sum_{j=1}^N x(j)\exp \left(\frac{-\pi i}{2N}(2j-1)(k-1)\right)
と展開されます。j = N + 1 以降は xx = 0 であることを反映させました。また虚部だけを取ります。
{{{\rm Im}(Y(k)w(k))}}=-\sum_{j=1}^N x(j)\sin \left(\frac{\pi }{2N}(2j-1)(k-1)\right)
$k$ について微調整すると DST-II の定義式に近づきます。
{{{\rm Im}(Y(k+1)w(k+1))}}=-\sum_{j=1}^N x(j)\sin \left(\frac{\pi }{2N}(2j-1)k\right)
DST-IIの定義式
{{\rm DST-II}}:y(k)=\sqrt{\frac{2}{N}}\sum_{j=1}^N x(j)\frac{1}{\sqrt{1+\delta_{kN} }}\sin \left(\frac{\pi }{2N}(2j-1)k\right)
と比較して最終的に以下の通り DST-II と DFT の関係が求まりました。
y(k)=-{{{\rm Im}(Y(k+1)w(k+1))}}\sqrt{\frac{2}{N\left(1+\delta_{kN} \right)}}
DST-III の計算
最後は DST-III です。こちらは DST-II の逆変換。この場合は、 $x(1)$ から $x(N)$ を以下のように拡張すると考えます。
0,x(1),x(2),x(3),\cdots ,x(N)/\sqrt{2},0,0,\cdots ,0
全長 $2N$です。$x(N)$だけ $\sqrt{2}$ で割っておきます。DST-I と同じ展開ですが、今度は $xx(j)$ に重みをかけた $xx(j)w(j)$ の DFT から始めます。また $w(j)$は以下のように定義します;
w(j)=\exp \left(\frac{-\pi i}{2N}(j-1)\right)
重みを付けているので $Y^{\prime } (k)$とおきます。
\begin{array}{cl}
Y^{\prime } (k) & =\sum_{j=1}^{2N} xx(j)w(j)\exp \left(\frac{-2\pi i}{2N}(j-1)(k-1)\right),\\
& =\sum_{j=1}^{N+1} xx(j)\exp \left(\frac{-\pi i}{2N}(j-1)(2k-1)\right)
\end{array}
虚部を取ると
\begin{array}{cl}
{{\rm Im}}(Y^{\prime } (k)) & =-\sum_{j=1}^{N+1} xx(j)\sin \left(\frac{\pi }{2N}(j-1)(2k-1)\right)
\end{array}
$j$ について微調整すると
\begin{array}{cl}
{{\rm Im}}(Y^{\prime } (k)) & =-\sum_{j=0}^N xx(j+1)\sin \left(\frac{\pi }{2N}j(2k-1)\right)\\
& =-\sum_{j=1}^N xx(j+1)\sin \left(\frac{\pi }{2N}j(2k-1)\right)\\
& =-\sum_{j=1}^N x(j)\frac{1}{\sqrt{1+\delta_{jN} }}\sin \left(\frac{\pi }{2N}j(2k-1)\right)
\end{array}
j = 0 の時 $\sin$ の項はすべて 0 になるので除外しました。DST-III の定義式
{{\rm DST-III}}:y(k)=\sqrt{\frac{2}{N}}\sum_{j=1}^N x(j)\frac{1}{\sqrt{1+\delta_{jN} }}\sin \left(\frac{\pi }{2N}j(2k-1)\right)
と比較すると最終的に以下の通り DST-III と DFT の関係が求まります。
y(k)=-{{\rm Im}}(Y(k))\sqrt{\frac{2}{N}}
まとめ
これで安心して離散サイン変換 DST-I, DST-II, DST-III が計算できるようになりました。
ポワソン方程式の数値解法を例に「いろんな状況でポワソン方程式を高速に解く 」で活用しています。