前言
暂略。
前置知识
需要一定编程基础，理解连续/离散傅里叶变换原理与计算方法。
正文
首先，让我们来看一下DCT的公式：
$$ F(u)= \begin{cases} \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}f(x), & n\ne m \\ \frac{2}{\sqrt{N}}\sum_{x=0}^{N-1}f(x)\cos\frac{(2x+1)u\pi}{2N}, &n=m \end{cases} $$
具体推导过程可见前文。
上式即为DCT的公式表达。而在具体到计算时，例如要处理一个长度为8的时域信号时，可以展开为：
通过观察我们可以发现在变换的具体表达上很像矩阵相乘的样子。因此，我们可以进一步改写为：
其中，最左边的$[F(u)]$为变换后的矩阵，即为频域的矩阵。最右边的$f[x]$为变换前的时域矩阵。而中间的矩阵在N确定时为一确定矩阵，我们把它称之为DCT的变换矩阵A。
而A可以写为：
$$ A[x,y]= \begin{cases} \frac{1}{\sqrt{N}}, & y= 0 \\ \frac{2}{\sqrt{N}}\cos\frac{(2x+1)y\pi}{2N}, &y\ne 0 \end{cases} $$
所以，我们可以把一维的DCT变换表现为矩阵相乘的形式，即：$F=A[f(x)]$。
这是一维的情况，那么二维的呢？
我们首先来看下二维离散余弦变换的一般公式（为了写的方便稍微换了下写法）：
$$ F(u,v)=\alpha(u)\alpha(v)\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)\cos\frac{(2x+1)u\pi}{2M}\cos\frac{(2y+1)v\pi}{2N} $$
其中：
$$ \alpha(u)= \begin{cases} \frac{1}{\sqrt{N}}, & u= 0 \\ \frac{2}{\sqrt{N}}, &u\ne 0 \end{cases} \alpha(v)= \begin{cases} \frac{1}{\sqrt{N}}, & v= 0 \\ \frac{2}{\sqrt{N}}, &v\ne 0 \end{cases} $$
可以看到，二维的DCT表达起来和一维的DCT变换差不多，基本上就是在原始数据的x方向做一次一维DCT，然后再在y方向再做一次一维 DCT。而事实上也就是那样。
所以对于一个指定的$[f(x,y)]$，我们要做的就是对它的每一列都做一次DCT变换，再转置一次后再做一次DCT变换，最后再转置回来，就完成了二维的DCT。
具体可以表示为：
$$ [F(u,v)]=(A((A[f(x,y)])^T)) ^T $$
再因为转置矩阵的性质：
$$ (AB)^T = B^T A^T $$
所以我们可以进一步化简：
$$ [F(u,v)]=(A((A[f(x,y)])^T)) ^T\\ =(A([f(x,y)]^T A^T)) ^T\\ =([f(x,y)]^T A^T) ^T A^T\\ =A[f(x,y)]A^T $$
(写的格式有点问题见谅。。）
所以我们就得到了二维DCT变换公式：
$$ F(u,v)=A[f(x,y)]A^T $$
编码实现
为了进一步巩固学习内容，这里使用C语言写一个小DEMO，实现二维DCT的效果。代码除了math库和stdio库外不包含任何第三方库。编写完成后会使用matlab对所得结果进行验证。
首先，进行DCT需要用到矩阵运算，这里先给出简易的矩阵运算代码。
typedef struct {
int rows;
int cols;
double** p;
} Matrix;
Matrix initMatrix(int rows, int cols) {
Matrix mat;
mat.rows = rows;
mat.cols = cols;
mat.p = (double**)malloc(sizeof(double*)*rows);
// new double*[rows];//分配rows_num个指针
for (int i = 0; i < rows; ++i) {
mat.p[i] = (double*)malloc(sizeof(double)*cols);
// new double[cols];//为p[i]进行动态内存分配，大小为cols
}
return mat;
}
Matrix initMatrixWithValue(int rows, int cols, double* p) {
Matrix mat = initMatrix(rows, cols);
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
mat.p[i][j] = *p++;
}
}
return mat;
}
Matrix matrixTrans(Matrix mat) {
Matrix res = initMatrix(mat.cols, mat.rows);
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
res.p[j][i] = mat.p[i][j];
}
}
return res;
}
Matrix matrixMul(Matrix m1, Matrix m2){
Matrix res = initMatrix(m1.rows, m2.cols);
for (int i = 0; i < m1.rows; i++) {
for (int j = 0; j < m2.cols; j++) {
res.p[i][j] = 0;
for (int k = 0; k < m1.cols; k++) {
res.p[i][j] += (m1.p[i][k] * m2.p[k][j]);
}
}
}
return res;
}
//矩阵显示
void showMatrix(Matrix mat) {
//cout << rows_num <<" "<<cols_num<< endl;//显示矩阵的行数和列数
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
printf("%.4f ", mat.p[i][j]);
}
printf("\n");
}
printf("\n");
}上述代码实现了矩阵的初始化，加法，乘法和转置运算。
Matrix get4x4DCTMatrixL(){
Matrix mat = initMatrix(4, 4);
for(int i = 0;i < 4;i++){
for(int j = 0;j < 4;j++){
double c = (i==0) ? sqrt(1.0/4) : sqrt(2.0/4);
mat.p[i][j] = c * cos( (j + 0.5) * PI * i / 4);
}
}
return mat;
}这段代码即为生成DCT变换矩阵的代码。这里是一个处理4x4的数据的矩阵。之后若是要进行DCT变换，只要按照上文去做计算即可。
为了进行测试，这里我先用matlab生成一段数据用于测试。代码如下：
clc;clear
f = (rand(4,4)*100);
dct2(f);f =
40.1808 18.3908 90.2716 33.7719
7.5967 23.9953 94.4787 90.0054
23.9916 41.7267 49.0864 36.9247
12.3319 4.9654 48.9253 11.1203
ans =
156.9409 -54.8586 -28.9792 51.3964
43.0922 -19.6215 -0.1741 18.9064
-26.9619 28.4906 -3.5949 26.3422
-6.7750 39.6837 -3.5258 -9.3417接下来我们就在c语言中对上述数据进行计算，并对答案进行比较。代码如下：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926f
typedef struct {
int rows;
int cols;
double** p;
} Matrix;
Matrix initMatrix(int rows, int cols) {
Matrix mat;
mat.rows = rows;
mat.cols = cols;
mat.p = (double**)malloc(sizeof(double*)*rows);
// new double*[rows];//分配rows_num个指针
for (int i = 0; i < rows; ++i) {
mat.p[i] = (double*)malloc(sizeof(double)*cols);
// new double[cols];//为p[i]进行动态内存分配，大小为cols
}
return mat;
}
Matrix initMatrixWithValue(int rows, int cols, double* p) {
Matrix mat = initMatrix(rows, cols);
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
mat.p[i][j] = *p++;
}
}
return mat;
}
Matrix matrixTrans(Matrix mat) {
Matrix res = initMatrix(mat.cols, mat.rows);
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
res.p[j][i] = mat.p[i][j];
}
}
return res;
}
Matrix matrixMul(Matrix m1, Matrix m2){
Matrix res = initMatrix(m1.rows, m2.cols);
for (int i = 0; i < m1.rows; i++) {
for (int j = 0; j < m2.cols; j++) {
res.p[i][j] = 0;
for (int k = 0; k < m1.cols; k++) {
res.p[i][j] += (m1.p[i][k] * m2.p[k][j]);
}
}
}
return res;
}
//矩阵显示
void showMatrix(Matrix mat) {
//cout << rows_num <<" "<<cols_num<< endl;//显示矩阵的行数和列数
for (int i = 0; i < mat.rows; i++) {
for (int j = 0; j < mat.cols; j++) {
printf("%.4f ", mat.p[i][j]);
}
printf("\n");
}
printf("\n");
}
Matrix get4x4DCTMatrixL(){
Matrix mat = initMatrix(4, 4);
for(int i = 0;i < 4;i++){
for(int j = 0;j < 4;j++){
double c = (i==0) ? sqrt(1.0/4) : sqrt(2.0/4);
mat.p[i][j] = c * cos( (j + 0.5) * PI * i / 4);
}
}
return mat;
}
int main(){
Matrix dctm = get4x4DCTMatrixL();
double testF_a[4*4] = {
40.1808, 18.3908, 90.2716, 33.7719,
7.5967, 23.9953, 94.4787, 90.0054,
23.9916, 41.7267, 49.0864, 36.9247,
12.3319, 4.9654, 48.9253, 11.1203
};
Matrix testF = initMatrixWithValue(4, 4, testF_a);
Matrix tmp = matrixMul(dctm, testF);
Matrix ans = matrixMul(tmp, matrixTrans(dctm));
showMatrix(ans);
return 0;
}我们把它编译运行试试：
[ame@RaMizy ws]$ gcc test_dct.c -o test_dct -lm && ./test_dct
156.9409 -54.8586 -28.9792 51.3964
43.0922 -19.6215 -0.1741 18.9063
-26.9619 28.4906 -3.5949 26.3423
-6.7750 39.6837 -3.5258 -9.3417 可以看到，我们的代码计算结果与matlab的结果在合理的误差范围内基本一致。结果存在的微小误差可认为来自于不同的计算方式及精度带来的截断误差。
结尾
至此，已经基本讲完了DCT的矩阵相关。但由于时间所限，我也不可能讲的很完整，也难免会有些问题。就之后再补充吧。
原本也想讲下像是蝶式算法之类的快速算法，就之后再细写吧。reference先放这里了。
Reference
二维离散余弦变换（2D-DCT）DCT变换展开时用到的示例图片
DCT变换自学笔记
离散余弦变换（DCT）的来龙去脉
图像的DCT算法 matlab代码参考
快速离散余弦变换代码实现(FDCT)
H.264中整数DCT变换，量化，反量化，反DCT究竟是如何实现的？
H264___DCT蝶形算法____理解
H264-整数DCT变换和蝶形变换代码实现
评论 (0)