離散フーリエ変換と各種離散サイン変換（DST-I ～ DST-III）


# はじめに


以前、離散サイン変換（DST-I）と `fft` との関係は確認しましたが、DST-II と DST-III のタイプの離散サイン変換も必要になったので、改めて `fft` 関数を使って計算する方法を確認してみました。




何通りもやり方がありますし各所にコードも転がっている（例えば [File Exchange: Wavelet Scale Spectra](https://jp.mathworks.com/matlabcentral/fileexchange/35092-wavelet-scale-spectra?s_eid=PSM_29435) ）んですが、やっぱり確認しておかないと・・ということで備忘録代わりに。



   -  参照１：[いろんな状況でポワソン方程式を高速に解く](https://qiita.com/eigs/items/fee977f4acbf1a4b3472) 
   -  参照２：[離散フーリエ変換から離散サイン変換までの式展開](https://qiita.com/eigs/items/940c61f07134bb62dd3b) 
   -  参照３：[離散フーリエ変換から離散コサイン変換までの式展開](https://qiita.com/eigs/items/ff8e77d1e5d8844b4236) 



ちなみに DST-I については MATLAB の Partial Differential Equation Toolbox に `dst` という関数があります。


# 定義式


まずは定義式から。長さ $N$の信号 $x$ に対して離散フーリエ変換を $Y(k)$、離散コサイン変換を $y(k)$ とおきます。$\delta_{kN}$ はクロネッカーデルタです。



```math
\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 に掲載されている定義式は以下の通り。



```math
\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 の計算


これは既に [離散フーリエ変換から離散サイン変換までの式展開](https://qiita.com/eigs/items/940c61f07134bb62dd3b) で一例を記載していますが、式展開がもう少し楽になる方法を改めて記載します。




変換対象の元データを x （サイズ N )とします。



```math
xx=[0,x(1),x(2),x(3),\cdots ,x(N),0,0,0,\cdots ,0,0]
```


と拡張した長さ $2N+2$ のベクトル $xx$ に対しての離散フーリエ変換から式を展開していきます。




上の FFT の定義式から



```math
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 の定義式に近づきます。



```math
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)
```


虚部を取ると



```math
\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 の定義



```math
y(k)=\sqrt{\frac{2}{N+1}}\sum_{j=1}^N x(j)\sin \left(\frac{\pi }{N+1}jk\right),
```


と比べて、最終的に以下の通り DST-I と DFT の関係が求まります。



```math
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 )とします。



```math
xx=[x(1),x(2),x(3),\cdots ,x(N),0,0,0,\cdots ,0,0]
```


全長 $2N$ です。長さ $2N$ の離散フーリエ変換を考えて



```math
w(k)=\exp \left(\frac{-\pi i}{2N}(k-1)\right),
```


という重みをかけた式を展開していきます。[半サンプル分データをシフト](https://dsp.stackexchange.com/questions/10250/how-to-circularly-shift-a-signal-by-a-fraction-of-a-sample)することに対応する処理と理解しています。上の定義式から



```math
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 であることを反映させました。また虚部だけを取ります。



```math
{{{\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 の定義式に近づきます。



```math
{{{\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の定義式



```math
{{\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 の関係が求まりました。



```math
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)$ を以下のように拡張すると考えます。



```math
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)$は以下のように定義します;



```math
w(j)=\exp \left(\frac{-\pi i}{2N}(j-1)\right)
```


重みを付けているので $Y^{\prime } (k)$とおきます。



```math
\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}
```


虚部を取ると



```math
\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$ について微調整すると



```math
\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 の定義式



```math
{{\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 の関係が求まります。



```math
y(k)=-{{\rm Im}}(Y(k))\sqrt{\frac{2}{N}}
```
  
# まとめ


これで安心して離散サイン変換 DST-I, DST-II, DST-III が計算できるようになりました。
ポワソン方程式の数値解法を例に「[いろんな状況でポワソン方程式を高速に解く](https://qiita.com/eigs/items/fee977f4acbf1a4b3472) 」で活用しています。