上节分享了二维图像离散傅里叶变换，本节来继续讲频域空间的另一种变换–二维离散余弦变换（Discrete Cosine Transform，DCT）。从运算方式上来讲，离散傅里叶变换计算的对象为复数，但离散余弦变换的对象为实数。虽然离散余弦变换没有离散傅里叶变换的功能强大，但是离散余弦变换的计算速度要比对象为复数的离散傅里叶变换快得多，并且已经被广泛应用到图像压缩编码、语音信号处理等众多领域。
01
—
离散余弦变换原理
一般来说，二维离散余弦变换的定义由如下公式表示：
相应的，二维离散余弦逆变换的定义由如下公式表示：
02
—
离散余弦变换Python实现
参考上述离散余弦变换给出的公式，我们来写其具体实现：
def DCT2D(x):'''Discrete space cosine transformx: Input matrix'''N1, N2 = x.shapeX = np.zeros((N1, N2))n1, n2 = np.mgrid[0:N1, 0:N2]for w1 in range(N1):for w2 in range(N2):l1 = (2/N1)**0.5 if w1 else (1/N1)**0.5l2 = (2/N2)**0.5 if w2 else (1/N2)**0.5cos1 = np.cos(np.pi*w1*(2*n1 + 1)/(2*N1))cos2 = np.cos(np.pi*w2*(2*n2 + 1)/(2*N2))X[w1, w2] = l1*l2*np.sum(x*cos1*cos2)return X
运行效果如下：
03
—
离散余弦逆变换Python实现
参考上述离散余弦逆变换给出的公式，我们实现其python代码：
def iDCT2D(X, shift=True):'''Inverse discrete space cosine transformX: Input spectrum matrix'''N1, N2 = X.shapex = np.zeros((N1, N2))k1, k2 = np.mgrid[0:N1, 0:N2]l1 = np.ones((N1, N2))*(2/N1)**0.5l2 = np.ones((N1, N2))*(2/N2)**0.5l1[0] = (1/N1)**0.5; l2[:,0] = (1/N2)**0.5for n1 in range(N1):for n2 in range(N2):cos1 = np.cos(np.pi*k1*(2*n1 + 1)/(2*N1))cos2 = np.cos(np.pi*k2*(2*n2 + 1)/(2*N2))x[n1, n2] = np.sum(l1*l2*X*cos1*cos2)return x
运行效果如下：
04
—
扩展
尽管DCT最常用的领域为图像压缩，但是我们可以尝试往经DCT变换后的图像中添加一些滤波器，再进行逆变换，我们来做个实验看一下其具体效果。
这里选用的filter如下：
测试代码如下：
u, v = np.mgrid[0:N1, 0:N2]/max(N1, N2)r = (u**2 + v**2)**0.5theta = np.arctan2(v, u)H = np.exp(-3*r**2)*(np.cos(4*2*theta)/2 + 1/2)image__ = iDCT2D(H*IMAGE)Hx__ = np.array([H, abs(image__*0.5 + 0.5)])panel(Hx__, (2, 1), text_color='green',texts=['Filter', 'Filtered image'])
运行结果如下：
05
—
总结
与离散傅里叶变换DFT相比，离散余弦变换DCT有两个主要优点：
DCT是一种比 DFT 具有更好的计算效率的实变换，DFT的定义是一个复变换。
DCT不引入间断，同时在时域信号中施加周期性。如上图所示：在 DFT 中，当时间信号被截断并假定为周期性时，在时域中引入了不连续性，并引入了相应的频域操作。但是，当截断时域信号时，即使假设对称性，也不会在 DCT 中引入不连续性。
注：完整代码可公众号内回复 DCT 即可获取