Discrete sine transform¶
Trigonometric relation¶
The relation
\[\stwiddle{
2 \pi
}{
\left( m + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}
=
\left( - 1 \right)^{l + 1}
\stwiddle{
2 \pi
}{
\left( n + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}\]
where \(m \equiv N - 1 - n\) is repeatedly used in this section, where the symbols \(n, k, l\) are all integers.
Derivation
\[ \begin{align}\begin{aligned}\stwiddle{
2 \pi
}{
\left( m + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}
&
=
\stwiddle{
2 \pi
}{
\left\{ \left( N - 1 - n \right) + \frac{1}{2} \right\}
\left( 2 k + l \right)
}{
2 N
}\\&
=
\stwiddle{
2 \pi
}{
\left\{ N - \left( n + \frac{1}{2} \right) \right\}
\left( 2 k + l \right)
}{
2 N
}\\&
=
\sin
\left(
\left( 2 k + l \right) \pi
-
2
\pi
\frac{
\left( n + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}
\right).\end{aligned}\end{align} \]
By utilizing
\[ \begin{align}\begin{aligned}\sin \left( \alpha - \beta \right)
&
=
\sin \alpha \cos \beta
-
\cos \alpha \sin \beta,\\\sin \left( \left( 2 k + l \right) \pi \right)
&
=
0,\\\cos \left( \left( 2 k + l \right) \pi \right)
&
=
\left( - 1 \right)^l,\end{aligned}\end{align} \]
we end up with
\[\stwiddle{
2 \pi
}{
\left( m + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}
=
\left( - 1 \right)^{l + 1}
\stwiddle{
2 \pi
}{
\left( n + \frac{1}{2} \right)
\left( 2 k + l \right)
}{
2 N
}.\]
Forward transform¶
We consider the forward transform (discrete sine transform of type II):
\[X_k
\equiv
2
\sum_{n = 0}^{N - 1}
x_n
\stwiddle{
2 \pi
}{
\left( n + \frac{1}{2} \right)
\left( k + 1 \right)
}{
2 N
}
\equiv
\dstii{k}{N}{x_0}{x_1}{x_{N - 1}}\]
for \(\seq{k}{0}{1}{N - 1}\).
Assuming that \(N\) is a multiple of \(2\), we consider the even and odd wavenumbers separately:
\[ \begin{align}\begin{aligned}X_{2 k}
=
&
2
\sum_{n = 0}^{N - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N},\\X_{2 k + 1}
=
&
2
\sum_{n = 0}^{N - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N},\end{aligned}\end{align} \]
for \(\seq{k}{0}{1}{N / 2 - 1}\).
The summation with respect to \(n\) is split into two parts as well following
\[ \begin{align}\begin{aligned}\sum_{n = 0}^{N - 1}
q_n
=
&
\sum_{n = 0}^{N / 2 - 1}
q_n
+
\sum_{n = N / 2}^{N - 1}
q_n,\\=
&
\sum_{n = 0}^{N / 2 - 1}
q_n
+
\sum_{n = 0}^{N / 2 - 1}
q_{N - 1 - n},\end{aligned}\end{align} \]
giving
\[ \begin{align}\begin{aligned}X_{2 k}
=
&
2
\sum_{n = 0}^{N / 2 - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
+
2
\sum_{n = 0}^{N / 2 - 1}
x_{N - 1 - n}
\stwiddle{2 \pi}{\left( N - 1 - n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\=
&
2
\sum_{n = 0}^{N / 2 - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
+
2
\sum_{n = 0}^{N / 2 - 1}
x_{N - 1 - n}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\\=
&
2
\sum_{n = 0}^{N / 2 - 1}
\left(
x_n
+
x_{N - 1 - n}
\right)
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N},\end{aligned}\end{align} \]
and
\[ \begin{align}\begin{aligned}X_{2 k + 1}
=
&
2
\sum_{n = 0}^{N / 2 - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}
+
2
\sum_{n = 0}^{N / 2 - 1}
x_{N - 1 - n}
\stwiddle{2 \pi}{\left( N - 1 - n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}\\=
&
2
\sum_{n = 0}^{N / 2 - 1}
x_n
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}
-
2
\sum_{n = 0}^{N / 2 - 1}
x_{N - 1 - n}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}\\=
&
2
\sum_{n = 0}^{N / 2 - 1}
\left(
x_n
-
x_{N - 1 - n}
\right)
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 \left( N / 2 \right)}\\=
&
\dstii{k}{N / 2}{x_0 - x_{N - 1}}{x_1 - x_{N - 2}}{x_{N / 2 - 1} - x_{N / 2}}\end{aligned}\end{align} \]
for \(\seq{k}{0}{1}{N / 2 - 1}\). Note that we adopt the trigonometric relations derived before.
By utilizing the product-to-sum identity:
\[\sin \alpha
\equiv
\frac{
\sin \left( \alpha + \beta \right)
+
\sin \left( \alpha - \beta \right)
}{
2 \cos \beta
}\]
with
\[\beta
=
2 \pi
\frac{
n + \frac{1}{2}
}{
2 N
},\]
we obtain
\[ \begin{align}\begin{aligned}X_{2 k}
=
&
2
\sum_{n = 0}^{N / 2 - 1}
\frac{
x_n
+
x_{N - 1 - n}
}{
2 \cos \beta
}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) k}{2 \left( N / 2 \right)}
+
2
\sum_{n = 0}^{N / 2 - 1}
\frac{
x_n
+
x_{N - 1 - n}
}{
2 \cos \beta
}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 \left( N / 2 \right)}\\=
&
\dstii{
k - 1
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
}
+
\dstii{
k
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
},\end{aligned}\end{align} \]
which can be applied directly for \(\seq{k}{1}{2}{N / 2 - 1}\).
For \(k = 0\) (i.e. a corner case), assigning \(k = 0\) to the first term reveals
\[\dstii{
- 1
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
}
=
0.\]
In summary, we obtain the following recurrence relation:
\[ \begin{align}\begin{aligned}X_{2 k}
&
=
\dstii{
k - 1
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
}\\&
+
\dstii{
k
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
},\\X_{2 k + 1}
&
=
\dstii{
k
}{
N / 2
}{
x_{0} - x_{N - 1}
}{
x_{1} - x_{N - 2}
}{
x_{N / 2 - 1} - x_{N / 2}
},\end{aligned}\end{align} \]
which is valid for \(\seq{k}{0}{1}{N / 2 - 1}\). Recall that
\[\dstii{
- 1
}{
N / 2
}{
\frac{
x_0
+
x_{N - 1}
}{
2 \cos \beta
}
}{
\frac{
x_1
+
x_{N - 2}
}{
2 \cos \beta
}
}{
\frac{
x_{N / 2 - 1}
+
x_{N / 2}
}{
2 \cos \beta
}
}
=
0.\]
Note that, \(X_0 = 2 x_0\) is satisfied for \(N = 1\).
src/dst.c¶
37static int dst2 (
38 const size_t nitems,
39 const size_t stride,
40 const double * const restrict table,
41 double * const restrict xs,
42 double * const restrict ys
43) {
44 if (1 == nitems) {
45 xs[0] *= 2.;
46 return 0;
47 }
48 if (0 == nitems % 2) {
49 // divide and conquer
50 const size_t nh = nitems / 2;
51 double * const buf0 = ys + 0;
52 double * const buf1 = ys + nh;
53 // create input buffers of DST-II
54 for (size_t n = 0; n < nh; n++) {
55 // c: 1 / [ 2 cos(beta) ]
56 const double c = table[(2 * n + 1) * stride];
57 const double value0 = xs[ n];
58 const double value1 = xs[nitems - 1 - n];
59 buf0[n] = c * (value0 + value1);
60 buf1[n] = 1. * (value0 - value1);
61 }
62 // solve two sub problems
63 dst2(nh, stride * 2, table, buf0, xs);
64 dst2(nh, stride * 2, table, buf1, xs);
65 // distribute results
66 for (size_t k = 0; k < nh; k++) {
67 // to even frequencies
68 // for the first element, "-1"-th element is zero
69 const double value0 = 0 == k ? 0. : buf0[k - 1];
70 const double value1 = buf0[k ];
71 xs[k * 2 + 0] = value0 + value1;
72 // to even frequencies
73 xs[k * 2 + 1] = buf1[k];
74 }
75 } else {
76 // fallback to N^2 DST2
77 for (size_t k = 0; k < nitems; k++) {
78 double * const y = ys + k;
79 *y = 0.;
80 for (size_t n = 0; n < nitems; n++) {
81 const double phase = 2. * pi * (n + 0.5) * (k + 1) / (2 * nitems);
82 *y += 2. * xs[n] * sin(phase);
83 }
84 }
85 for (size_t k = 0; k < nitems; k++) {
86 xs[k] = ys[k];
87 }
88 }
89 return 0;
90}
Backward transform¶
We consider the inverse transform (discrete sine transform of type III):
\[x_n
=
\frac{1}{N}
\left\{
\sum_{k = 0}^{N - 2}
\hat{X}_k
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 N}
+
\left( - 1 \right)^n
\frac{\hat{X}_{N - 1}}{2}
\right\}\]
for \(\seq{n}{0}{1}{N - 1}\). For later convenience, we define
\[\begin{split}X_k
=
\begin{cases}
\frac{\hat{X}_{N - 1}}{2} & k = N - 1, \\
\hat{X}_k & \text{otherwise},
\end{cases}\end{split}\]
to write the inverse transform as
\[x_n
=
\frac{1}{N}
\sum_{k = 0}^{N - 1}
X_k
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 N}
\equiv
\dstiii{n}{N}{X_0}{X_1}{X_{N - 1}}\]
for \(\seq{n}{0}{1}{N - 1}\).
Note that we perform this conversion a priori:
src/dst.c¶
215xs[nitems - 1] *= 0.5;
216dst3(nitems, 1, table, xs, ys);
Assuming that \(N\) is a multiple of \(2\), we decompose the right-hand-side terms into even and odd components:
\[ \begin{align}\begin{aligned}x_n
&
=
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
+
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k + 1}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}.\\&
=
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
+
\frac{1}{2}
\dstiii{n}{N / 2}{X_1}{X_3}{X_{N - 1}}.\end{aligned}\end{align} \]
By introducing \(m = N - 1 - n\) to utilize the symmetric nature of \(x_n\), we see
\[ \begin{align}\begin{aligned}x_{N - 1 - n}
&
=
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( m + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
+
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k + 1}
\stwiddle{2 \pi}{\left( m + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}\\&
=
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}
-
\frac{1}{2}
\dstiii{n}{N / 2}{X_1}{X_3}{X_{N - 1}}.\end{aligned}\end{align} \]
By adopting the product-to-sum identity:
\[\sin \alpha
\equiv
\frac{
\sin \left( \alpha + \beta \right)
+
\sin \left( \alpha - \beta \right)
}{
2 \cos \beta
}\]
with
\[\beta
=
2 \pi
\frac{
n + \frac{1}{2}
}{
2 N
},\]
the first term
\[\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 1 \right)}{2 N}\]
leads to
\[ \begin{align}\begin{aligned}&
\frac{1}{2 \cos \beta}
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k \right)}{2 N}
+
\frac{1}{2 \cos \beta}
\frac{1}{N}
\sum_{k = 0}^{N / 2 - 1}
X_{2 k}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( 2 k + 2 \right)}{2 N}\\=
&
\frac{1}{4 \cos \beta}
\frac{1}{N / 2}
\sum_{k = - 1}^{N / 2 - 2}
X_{2 k + 2}
\stwiddle{2 \pi}{\left( n + \frac{1}{2} \right) \left( k + 1 \right)}{2 \left( N / 2 \right)}
+
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}}\\=
&
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_2}{X_4}{X_N}
+
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_0}{X_2}{X_{N - 2}}\\=
&
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_0 + X_2}{X_2 + X_4}{X_{N - 2} + X_N},\end{aligned}\end{align} \]
where we let \(X_N = 0\).
In summary, we obtain the following recurrence relation:
\[ \begin{align}\begin{aligned}x_n
&
=
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_0 + X_2}{X_2 + X_4}{X_{N - 2} + X_N}
+
\frac{1}{2}
\dstiii{n}{N / 2}{X_1}{X_3}{X_{N - 1}},\\x_{N - 1 - n}
&
=
\frac{1}{4 \cos \beta}
\dstiii{n}{N / 2}{X_0 + X_2}{X_2 + X_4}{X_{N - 2} + X_N}
-
\frac{1}{2}
\dstiii{n}{N / 2}{X_1}{X_3}{X_{N - 1}},\end{aligned}\end{align} \]
for \(\seq{n}{0}{1}{N / 2 - 1}\). Note that \(x_0 = X_0\) for \(N = 1\).
src/dst.c¶
93static int dst3 (
94 const size_t nitems,
95 const size_t stride,
96 const double * const restrict table,
97 double * const restrict xs,
98 double * const restrict ys
99) {
100 if (1 == nitems) {
101 return 0;
102 }
103 if (0 == nitems % 2) {
104 // divide and conquer
105 const size_t nh = nitems / 2;
106 double * const buf0 = ys + 0;
107 double * const buf1 = ys + nh;
108 // create input buffers of DST-III
109 for (size_t k = 0; k < nh; k++) {
110 const double value0 = xs[k * 2 ];
111 const double value1 = nh - 1 == k ? 0. : xs[k * 2 + 2];
112 buf0[k] = value0 + value1;
113 buf1[k] = xs[k * 2 + 1];
114 }
115 // solve two sub problems
116 dst3(nh, stride * 2, table, buf0, xs);
117 dst3(nh, stride * 2, table, buf1, xs);
118 // combine results of sub problems
119 for (size_t n = 0; n < nh; n++) {
120 // c: 1 / [ 2 cos(beta) ]
121 const double c = table[(2 * n + 1) * stride];
122 const double value0 = 0.5 * c * buf0[n];
123 const double value1 = 0.5 * buf1[n];
124 xs[ n] = value0 + value1;
125 xs[nitems - 1 - n] = value0 - value1;
126 }
127 } else {
128 // fallback to N^2 DST3
129 for (size_t n = 0; n < nitems; n++) {
130 double * const y = ys + n;
131 *y = 0.;
132 for (size_t k = 0; k < nitems; k++) {
133 const double phase = 2. * pi * (n + 0.5) * (k + 1) / (2 * nitems);
134 *y += xs[k] * sin(phase);
135 }
136 *y /= 1. * nitems;
137 }
138 for (size_t n = 0; n < nitems; n++) {
139 xs[n] = ys[n];
140 }
141 }
142 return 0;
143}
Reference¶
Lee, “A New Algorithm to Compute the Discrete Cosine Transform”, IEEE T. Acoust. Speech, 1984