二维DCT变换 | Python实现

引言

最近专业课在学信息隐藏与数字水印，上到了变换域隐藏技术，提到了其中的DCT变换，遂布置了一个巨烦人的作业，让手动给两个\(8\times8\)的矩阵做二维DCT变换，在苦逼的算了一小时后，我决定放弃，转而决定写脚本来解决，\((๑•̀ㅂ•́)و✧\)，正好看网上好像只有matlab的脚本，好像没人用Python来写这个，遂打算搞一个（你就是纯粹为了偷懒不做作业\((*￣rǒ￣)）\)

二维DCT变换原理

还是要普及一下的嘛，毕竟让我头疼了一下午的东西，当然也要好好给你们分享一下啦ԅ(¯﹃¯ԅ)
DCT(Discrete Cosine Transform)，又叫离散余弦变换，它的第二种类型，经常用于信号和图像数据的压缩。经过DCT变换后的数据能量非常集中，一般只有左上角的数值是非零的，也就是能量都集中在离散余弦变换后的直流和低频部分
1、一维DCT变换
要弄懂二维离散余弦变换，首先我们需要先了解它在一维下的情况，具体表达式如下:
\(F(0)=\dfrac{1}{\sqrt{N}}\sum_{x=0}^{N-1} f(x)\) \(...........1\)
\(F(u)=\sqrt{\dfrac{2}{N}}\sum_{x=0}^{N-1} f(x) \cos{\dfrac{(2x+1)u\pi}{2N}}\) \(...........2\)
式中\(F(u)\)是第\(u\)个余弦变换值，\(u\)是广义频率变量，\(u=1,2,….,N-1;f(x)\)是时域\(N\)点序列。\(x= 1,2,….,N-1;\)
2、二维DCT变换
二维离散余弦变换可由下列表达式表示：
\(F(0,0)=\dfrac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y)\) \(...........3\)
\(F(0,v)=\dfrac{\sqrt{2}}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2y+1)v\pi}{2N}}\) \(...........4\)
\(F(u,0)=\dfrac{\sqrt{2}}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}\) \(...........5\)
\(F(u,v)=\dfrac{2}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}\cdot \cos{\dfrac{(2y+1)v\pi}{2N}}\) \(...........6\)
\(6\)是二维离散余弦变换的正变换公式，其中\(f(x,y)\)是空间域一个\(N*N\)的二维向量元素，即一个\(N*N\)的矩阵，\(x,y = 0,1,2，…，N-1;F(U,V)\)是经计算后得到的变换域矩阵，\(u,v = 0,1,2，….，N-1\).求和可分性是二维离散余弦变换的一个重要特征，因此我们可以用下式表示\(6\)：
\(F(u,v)=\dfrac{2}{N}\sum_{x=0}^{N-1} \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}{\sum_{y=0}^{N-1} f(x,y) }\cdot \cos{\dfrac{(2y+1)v\pi}{2N}} ..........7\)
由一维和二维的离散余弦变换公式性质可以推导得到二维离散余弦变换也可以写成矩阵相乘形式:
\(F=A[f(x,y)A^{T}]..............8\)
\(A\)为一维离散余弦变换的变换系数矩阵，\(A^{T}\)是\(A\)的转置矩阵
对图像进行二维离散余弦变换\(（2D-DCT）\)的步骤
1.获得图像的二维数据矩阵\(f(x,y)\)；
2.求离散余弦变换的系数矩阵\(A\)；
3.求系数矩阵对应的转置矩阵\(A^{T}\)；
4.根据公式\(F=A\times[f(x,y)]\times A^{T}\)计算离散余弦变换；
注：公式的大致推导,由于\(A\)都是正交阵，所以\(A^{-1}=A^{T}\),故\(Y=A\times[f(x,y)]\times A^{T} ，f(x,y)=A^{T}\times Y \times A=A^{T}\times A\times[f(x,y)]\times A^{T}\times A\)

参考链接：`https://www.cnblogs.com/latencytime/p/10228938.html`

参考链接：`https://blog.csdn.net/allen_sdz/article/details/83279210`

Python编程实现

大家注意上述的第\(8\)个式子，将变换的两个\(\sum_{i=0}^n\)转变成了变换矩阵和转置矩阵以及代转换矩阵之间乘积的问题（注意这里的乘是指的叉乘）
遂可以发现\(DCT\)后的矩阵应等于\(A\times X \times A^{T}\)
此外，对于\(x,y\)同时为0的情况，其参数要单独考虑
大致的思路就是，由于直接计算变换比较繁琐，所以我们就先对于一个单位阵进行操作运算，将其变成一个DCT的变换矩阵，而后在与代算矩阵和转置矩阵叉乘即可。
至此我们就可以着手写脚本了，这里主要是用了两个库numpy和math

# -*- coding:utf-8 -*-
# Author:Konmu
# DCT二维变换

from numpy import array as matrix, arange,zeros,transpose,matmul,ones
from math import sqrt,cos,pi

'''
作业代转化矩阵1
a=matrix([[0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0]])
'''
a=ones((8,8))#生成单位阵
for i in range(8):
    a[i]=matrix([128]*8)
# 生成全是128的矩阵（作业代转化矩阵2）
'''
测试数据
a=matrix([[61,19,50,20],
          [82,26,61,45],
          [89,90,82,43],
          [93,59,53,97]])
'''
A=zeros((8,8))#生成0矩阵
shape=a.shape[1]#获取维数
for i in range(8):
    for j in range(8):
        if(i == 0):
            x=sqrt(1/shape)
        else:
            x=sqrt(2/shape)
        A[i][j]=x*cos(pi*(j+0.5)*i/shape)#与维数相关

A_T=A.transpose()#矩阵转置
Y1=matmul(A,a)#矩阵叉乘
Y=matmul(Y1,A_T)
print(Y)
'''
想要近似值可以尝试这样输出
for i in range(shape):
    for j in range(shape):
        print('{:^8.4f}'.format(Y[i][j]),end='\n')
print()
'''