python离散余弦变换(DCT)

变换类型

离散余弦变换(DCT)相当于是傅里叶变换的实部，在scipy.fft中提供了dct函数及其逆变换idct。以dct为例，其函数定义如下，其中type表示余弦变换的类别，norm表示归一化模式。

dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, orthogonalize=None)

dct支持4种类型的余弦变换，下面列举出type分别为1,2,3,4时，且norm选用不同参数时，对应的DCT的具体形式。

• type=1

$$y[k]=x_0+(-1{)}^k{x}_{N-1}+2\sum _{n=1}^{N-2}x[n]\cos (\frac{\pi nk}{N-1}),0⩽k<N$$

• type=2

若norm=None，则

$$y[k]=2\sum _{n=0}^{N-1}x[n]\cos (\frac{\pi (2n+1)k}{2N}),0⩽k<N$$

当norm='ortho’时，表示采用正交DCT，此时$y[k]$需要乘以一个尺度因子$f$

$$f=\{\begin{array}{rl}& \sqrt{1/(4N)},ifk=0\\ & \sqrt{1/(2N)},otherwise\end{array}$$

此时，DCT的基函数${\varphi }_k[n]=2f\cos \frac{\pi (2n+1)k}{2N}$将彼此正交

$$\sum _{n=0}^{N-1}{\psi }_k[n]{\psi }_l[n]={\delta }_{lk}$$

• type=3

当norm=None时

$$y[k]=x_0+2\sum _{n=0}^{N-1}x[n]\cos (\frac{\pi n(2k+1)}{2N}),0⩽k<N$$

否则，若norm=‘ortho’，则

$$y[k]=\frac{x_0}{\sqrt{N}}+\frac{2}{\sqrt{N}}\sum _{n=0}^{N-1}x[n]\cos (\frac{\pi n(2k+1)}{2N}),0⩽k<N$$

• type=4

当norm=None时，

$$y[k]=2\sum _{n=0}^{N-1}x[n]\cos (\frac{\pi (2n+1)(2k+1)}{4N}),0⩽k<N$$

若norm=‘ortho’，则

$$y[k]=\sqrt{\frac{2}{N}}\sum _{n=0}^{N-1}x[n]\cos (\frac{\pi (2n+1)(2k+1)}{4N}),0⩽k<N$$

示例

示例如下

from scipy.fft import dct, idct
import matplotlib.pyplot as plt
N = 100
t = np.linspace(0,20,N, endpoint=False)
x = np.exp(-t/3)*np.cos(2*t)
y = dct(x, norm='ortho')

window = np.zeros(N)
window[:20] = 1
yr = idct(y*window, norm='ortho')
sum(abs(x-yr)**2) / sum(abs(x)**2)  # 0.0009872817275276098
plt.plot(t, x, '-bx')
plt.plot(t, yr, 'ro')
window = np.zeros(N)
window[:15] = 1
yr = idct(y*window, norm='ortho')
sum(abs(x-yr)**2) / sum(abs(x)**2)  #0.06196643004256714
plt.plot(t, yr, 'g+')
plt.legend(['x', '$x_{20}$', '$x_{15}$'])
plt.grid()
plt.show()

其中，$x(t)$是原始信号，$y(t)$是经过余弦变换之后的信号。window的作用是将$y$的部分数值置零，从而使用部分$y(t)$经由逆变换得到$x(t)$，最终结果如下。相比于保留20个数值的情况，只保留15个数值显然引入了更大的抖动。