3
$\begingroup$

I am trying to digest this paper.

[2002] Anna Usakova. Using Of Discrete Orthogonal Transforms For Convolution

So far, I've been successful with DFT and DHT. Could someone give a hand with DCT?

I am trying to reproduce the convolution theorem for DCT-I and DCT-II. But my results never match the same as the result I get with DFT and DHT.

enter image description here

I am reading reference 9 of the paper

RAO, K. R.—YIP, P.: Discrete Cosine Transform — Algo- rithms, Advantages, Applications, Academic Press, Inc., Lon- don, 1990.

But it puts the conv. theorem of DCT-II in other terms.

I have also read other works, such as

[1994] Martucci. Symmetric Convolution and the Discrete Sine and Cosine Transforms

[2011] Roma. A tutorial overview on the properties of the discrete cosine transform for encoded image and video processing

and others. But I am interested on the Conv. Theorem on a similar fashion as the first two references.

$\endgroup$
1
  • $\begingroup$ Did you ever solve this problem? If so, could you write your solution as an answer to your own question? I am very interested in the answer! $\endgroup$ – James Dec 2 '20 at 13:01
4
$\begingroup$

I did it for the DCT-I and DCT-II.

At first I thought it is about circular convolution, but it is not. After desperate attempts to do it by myself I found this article: https://dr.ntu.edu.sg/bitstream/10356/94137/1/Convolution%20Using%20Discrete%20Sine%20and%20Cosine%20Transforms.pdf.

In this paper there is a derivation of circular convolution using DTT's, but it is constructed from sum of $C_1^{-1}$ and $S_1^{-1}$, so, it's not what we are interested in.

What I suggest is to look at Fourier transform first, in order to get the idea:

Let we have

$\displaystyle F[t] = \sum_{k=0}^{N-1} f[k]e^{-\frac{i2π}{N}tk}$,

$\displaystyle G[t] = \sum_{k=0}^{N-1} g[k]e^{-\frac{i2π}{N}tk}$.

So then

$$\displaystyle ℱ\{f*g\}[t] = \sum_{x=0}^{N-1} f[x]e^{-\frac{i2π}{N}tx}\cdot\sum_{y=0}^{N-1} g[y]e^{-\frac{i2π}{N}ty} = \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f[x]g[y]e^{-\frac{i2π}{N}t(x+y)}$$

Now, let $x+y = m$, $m = 0\ldots2N-2$ and then group all $(x,y)$ which share the same $m$. So you will get:

$$\displaystyle \sum_{m=0}^{2N-2}\biggl(\sum_{0\leq x,y\leq N-1\\x+y=m}f[x]g[y] \biggl) e^{-\frac{i2π}{N}tm} \tag{1}$$

Because $x, y$ are both positive and less than $N$, then sum inside is just linear convolution at $m$. Let me denote it as $A[m]$, then everything is going to be rewritten in that form

$$\sum_{m=0}^{2N-2} A[m] e^{-\frac{i2π}{N}tm} \tag{2}$$

DFT kernel is N-periodic that is why we can group terms another time:

$$\displaystyle \sum_{m=0}^{N-1} (A[m] + A[m+N])e^{-\frac{i2π}{N}tm}$$

$A[m]$ captures $A[0]\ldots A[N-1]$, $A[m+N]$ captures $A[N]\ldots A[2N-1]$. But remember that $A[m]$ is linear convolution of two sequences of size $N$ and therefore size is $N + N - 1$ then $A[2N-1] = 0$. So these two sums are equivalent. Also it is a DFT of $A[m] + A[m+N]$. What is it? Right, circular convolution, as expected.

Now plug in DHT kernel instead of DFT into $(1)$ and you will get the same result ─ circular convolution. That is because DHT and DFT kernels are quite similar.

Well, if you have 2 signals $f$, $g$ with lengths $Q$ and $P$, and $A = f*g$, with $Q+P-1 = L$. Now, add zeros to $f$ up to $L$, do the same for $g$. Now $f$, $g$, and $A$ are $L$ periodic signals, so we can take the DFT of them. DFT of $A$ is going to be:

$$\sum_{m=0}^{L-1} A[m] e^{-\frac{i2π}{L}tm} = \sum_{m=0}^{L-1}\biggl(\sum_{0\leq x,y\leq L-1\\x+y=m}f[x]g[y] \biggl) e^{-\frac{i2π}{L}tm} = \sum_{x=0}^{L-1}\sum_{y=0}^{L-1}f[x]g[y]e^{-\frac{i2π}{L}t(x+y)} \tag{3}$$

So, via $(2)$ we can get the rule and expression for the convolution and $(3)$ shows that zero-padding up to convolution length is required to get linear convolution using the same rule. We're going to apply it to the DCT's and DST's.

Because there are many varieties of DCT's and DST's, I will mention what I used:

DCT-I (N+1) points:

$\displaystyle C_1\{f\}[n] = \sum_{k=0}^{N} c_k f[k]cos\Bigl[\frac{π}{N}nk\Bigl], \text{where}\ c_k = \left\{ \begin{array}{l} \frac{1}{2}, k≡0\ (mod\ N), \\0,\ \text{otherwise}\end{array} \right. \tag{4}$

DST-I (N+1) points: $$S_1\{f\}[n] = \sum_{k=0}^{N} f[k]sin\Bigl[\frac{π}{N}nk\Bigl]=\sum_{k=0}^{N} c_k f[k] sin\Bigl[\frac{π}{N}nk\Bigl] \tag{5}$$

(N-1) for data: $f[1] \ldots f[N-1]$ it doesn't matter what is placed in $f[0]$ and $f[N]$.

DCT-II (N+1) points:

$$C_2\{f\}[n] = \sum_{k=0}^{N} f[k]cos\Bigl[\frac{π}{N}\Bigl(k+\frac{1}{2}\Bigl)n\Bigl] \tag{6}$$

N for data: $f[0]\ldots f[N-1]$ and $f[N] = 0$: DST-II (N+1) points, N for data: $f[0]\ldots f[N-1]$ and $f[N] = 0$:

$$S_2\{f\}[n] = \sum_{k=0}^{N} f[k]sin\Bigl[\frac{π}{N}\Bigl(k+\frac{1}{2}\Bigl)n\Bigl] \tag{7}$$

Inverse DCT-I: $$C_1^{-1} = \frac{2}{N}C_1 \tag{8}$$ Inverse DST-I: $$S_1^{-1} = \frac{2}{N}S_1 \tag{9}$$ Inverse DCT-II (N+1) points: $$C_2^{-1}\{F_{C_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{C_2}[n]cos\Bigl[\frac{π}{N}n\Bigl(k+\frac{1}{2}\Bigl)\Bigl] \tag{10}$$ $C_2^{-1}\{C_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $f[N-1]$.

Inverse DST-II (N+1) points: $$C_2^{-1}\{F_{S_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{S_2}[n]sin\Bigl[\frac{π}{N}n\Bigl(k+\frac{1}{2}\Bigl)\Bigl] \tag{11}$$ $S_2^{-1}\{S_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $-f[N-1]$. Also, they are not interchangeable.

Everything mentioned above is just Wikipedia ones, but modified to achieve more similar look for easier coding and derivation.

So, now again, applying the rule $(2)$ for the DCT-I we get:

$$\sum_{m=0}^{2N} A[m] c_m cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{m=0}^{N} p_m(A[m] + A[2N-m]) c_m cos\Bigl[\frac{π}{N}nm\Bigl]$$ $p_m = \frac{1}{2}$ if $m = N$ else $p_m = 1$

from the other hand we have $$\sum_{m=0}^{2N} \frac{A[m]}{c_m} c_m cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]cos\Bigl[\frac{π}{N}n(x+y)\Bigl]$$ $$\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]cos\Bigl[\frac{π}{N}n(x+y)\Bigl] = \\\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\biggl(cos\Bigl[\frac{π}{N}nx\Bigl]cos\Bigl[\frac{π}{N}ny\Bigl] - sin\Bigl[\frac{π}{N}nx\Bigl]sin\Bigl[\frac{π}{N}ny\Bigl]\biggl) = \\ \biggl(\sum_{x=0}^{N} f[x]cos\Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]cos\Bigl[\frac{π}{N}ny\Bigl]\biggl) - \biggl(\sum_{x=0}^{N} f[x]sin \Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]sin\Bigl[\frac{π}{N}ny\Bigl]\biggl)$$

It is now easy to see that it can be rewriten in a such way:

$$C_1\{\frac{f[m]}{c_m}\}[n]\cdot C_1\{\frac{g[m]}{c_m}\}[n] - S_1\{f[m]\}[n]\cdot S_1\{g[m]\}[n]$$

Let $\frac{f[m]}{c_m} = f^{/}$ then we can write everything in this form: $$\frac{p_m}{c_m}(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$

alternative version $$δ_{0m}(A[0] + A[2N]) + (A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$ or even $$(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1} - \frac{f[0]g[0] + f[N]g[N]}{2}$$ Now, remember the $(3)$? If you add zeroes, so the size of f and g is equal to the convolution's length, then you have: $$(f*g)^{/}⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$

It looks a little bit different to what is presented in the paper, but it is because I used different version of DCT.

For DCT-II we have the rule (remember that $f[N] = g[N] = 0$):

$$\sum_{m=0}^{2N} A[m] cos\Bigl[\frac{π}{N}\Bigl(m+\frac{1}{2}\Bigl)n\Bigl] = \sum_{m=0}^{2N-2} A[m] cos\Bigl[\frac{π}{N}\Bigl(m+\frac{1}{2}\Bigl)n\Bigl] = \sum_{m=0}^{N-1} (A[m] + A[2N-1-m]) cos\Bigl[\frac{π}{N}\Bigl(m+\frac{1}{2}\Bigl)n\Bigl] $$

And the expression is going to be: $$\sum_{m=0}^{2N} A[m] cos\Bigl[\frac{π}{N}\Bigl(m+\frac{1}{2}\Bigl)n\Bigl] = \sum_{m=0}^{2N} A[m] \Bigl(cos\Bigl[\frac{πn}{2N}\Bigl]cos\Bigl[\frac{π}{N}kn\Bigl] - sin\Bigl[\frac{πn}{2N}\Bigl]sin\Bigl[\frac{π}{N}kn\Bigl]\Bigl) = \\ cos\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}) - sin\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G_{S_1} + F_{S_1}G^{/}_{C_1}) = \\ =F^{/}_{C_1}(cos\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1} - sin\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}) - F_{S_1}(cos\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}+ sin\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1}) = \\ =F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2} = G^{/}_{C_1}F_{C_2} - G_{S_1}F_{S_2} $$

Now we finally can write it as: $$A[m] + A[2N-1-m] ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$ And again, if you add enough zeroes, then it will be: $$f*g ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$

You can play with it here: https://repl.it/@dogtyping/dct#main.py — to make sure that it does work.

Thanks for the question, it was a lot of fun.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.