二维DCT变换

时间:2021-01-31 02:58:30

DCT(Discrete Consine Transform),又叫离散余弦变换,它的第二种类型,经常用于信号和图像数据的压缩。经过DCT变换后的数据能量非常集中,一般只有左上角的数值是非零的,也就是能量都集中在离散余弦变换后的直流和低频部分。

1. 一维DCT变换

一维DCT变换共有8中,其中最实用的是第二种形式,公式如下:

\[F(u)=c(u)\sum_{i=0}^{N-1}f(i)\cos{[\frac{(i+0.5)\pi}{N}u]}
\]

\[c(u)=\begin{cases}\sqrt{\frac{1}{N}}, & u=0 \\ \sqrt{\frac{2}{N}}, & u\neq0\end{cases}
\]

其中\(c(u)\)就是加上去一个系数,为了能使DCT变换变成正交矩阵。\(N\)是\(f(x)\)的总数。

2. 二维DCT变换

二维DCT变换是在一维的基础上再进行一次DCT变换,公式如下:

\[F(u,v)=c(u)c(v)\sum_{i=0}^{N-1}\sum_{i=0}^{N-1}f(i,j)\cos{[\frac{(i+0.5)\pi}{N}u]}\cos{[\frac{(i+0.5)\pi}{N}v]}
\]

这里只讨论了两个\(N\)相等的情况,也就是数据是方阵的形式,在实际应用中对不是方阵的数据都是先补齐再进行变换的。

写成矩阵形式:

\[\mathrm{F}=\mathrm{A}f\mathrm{A}^\mathrm{T}
\]

\[\mathrm{A}(i,j)=c(i)\cos{[\frac{(j+0.5)\pi}{N}i]}
\]

用MATLAB进行验证:

clear;
clc; X = round(rand(4) * 100); % 生成随机数据
A = zeros(4); % 变换矩阵 for i = 0 : 3
if i == 0
c = sqrt(1/4);
else
c = sqrt(2/4);
end
for j = 0 : 3
A(i + 1, j + 1) = c * cos(pi * (j + 0.5) * i / 4);
end
end Y = A * X * A'; % DCT变换
YY = dct2(X); % 使用MATALAB函数进行DCT变换 disp('使用公式进行DCT变换:')
disp(Y)
disp('使用MATLAB函数DCT变换:')
disp(YY)

输入结果:

使用公式进行DCT变换:
204.7500 -2.5322 27.2500 24.5909
32.1461 3.7448 -20.9667 24.5450
54.2500 -1.9287 -2.2500 -24.9079
12.9327 -40.4550 -25.1401 9.7552 使用MATLAB函数DCT变换:
204.7500 -2.5322 27.2500 24.5909
32.1461 3.7448 -20.9667 24.5450
54.2500 -1.9287 -2.2500 -24.9079
12.9327 -40.4550 -25.1401 9.7552

3. 二维DCT逆变换

DCT逆变换的公式如下:

\[f(i,j)=\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}c(u)c(v)F(u,v)\cos[\frac{(i+0.5)\pi}{N}u]\cos[\frac{(j+0.5)\pi}{N}v]
\]

\[c(u)=\begin{cases}\sqrt{\frac{1}{N}},&u=0\\\sqrt{\frac{2}{N},}&u\neq0\end{cases}
\]

矩阵形式的变换公式推到如下:

\[\begin{align}&\quad\mathrm{F}=\mathrm{A}f\mathrm{A}^{\mathrm{T}}\\&\because{\mathrm{A}^{-1}=\mathrm{A}^{\mathrm{T}}}\\&\therefore{f=\mathrm{A}^{-1}\mathrm{F}(\mathrm{A}^{T})^{-1}}=\mathrm{A}^{T}\mathrm{F}\mathrm{A}\end{align}
\]

用MATALAB进行验证:

clear;
clc; X = round(rand(4) * 100); % 生成随机数据
A = zeros(4); % 变换矩阵 for i = 0 : 3
if i == 0
c = sqrt(1/4);
else
c = sqrt(2/4);
end
for j = 0 : 3
A(i + 1, j + 1) = c * cos(pi * (j + 0.5) * i / 4);
end
end Y = A * X * A'; % DCT变换
XX = A'* Y* A; % DCT逆变换 disp('原始矩阵:')
disp(X)
disp('使用公式进行DCT逆变换:')
disp(XX)
disp('使用MATLAB函数DCT逆变换:')
disp(idct2(Y))

输出结果:

原始矩阵:
28 69 44 19
5 32 38 49
10 95 77 45
82 3 80 65 使用公式进行DCT逆变换:
28.0000 69.0000 44.0000 19.0000
5.0000 32.0000 38.0000 49.0000
10.0000 95.0000 77.0000 45.0000
82.0000 3.0000 80.0000 65.0000 使用MATLAB函数DCT逆变换:
28.0000 69.0000 44.0000 19.0000
5.0000 32.0000 38.0000 49.0000
10.0000 95.0000 77.0000 45.0000
82.0000 3.0000 80.0000 65.0000

4. DCT变换的可分离性

DCT变换是可分离的变换。通常根据可分离性,二维DCT可用两次一维DCT变换来完成,即

\[\begin{align}f(x,y)&\to F_{行}[f(x,y)]=F(x,v)\\&\to F(x,v)^{\mathrm{T}}\to F_{列}[f(x,v)^{\mathrm{T}}]=F(u,v)^{\mathrm{T}}\\&\to F(u,v)\end{align}
\]

先进行行变换,再进行列变换和先进行列变换,再进行行变换的结果是一样的。

Python scipy模块中的fftpack.dct()函数提供了一维DCT变换功能(默认是沿着矩阵的最后一个axis进行变换),下面使用Python代码进行验证。

import numpy as np
from scipy import fftpack def dct(mat2x2):
return fftpack.dct(fftpack.dct(mat2x2, norm='ortho').T, norm='ortho').T def dct2(mat2x2):
return fftpack.dct(fftpack.dct(mat2x2.T, norm='ortho').T, norm='ortho') if __name__ == '__main__':
sample = np.random.rand(3, 3)
print('先进行行变换,再进行列变换:')
print(dct(sample))
print('先进行列变换,再进行行变换:')
print(dct2(sample))

输出结果:

先进行行变换,再进行列变换:
[[ 1.3763706 -0.42355794 0.03903157]
[-0.18270004 0.06454257 -0.05273778]
[ 0.16962548 0.22247218 -0.06953193]]
先进行列变换,再进行行变换:
[[ 1.3763706 -0.42355794 0.03903157]
[-0.18270004 0.06454257 -0.05273778]
[ 0.16962548 0.22247218 -0.06953193]]

5. DCT用于图像压缩

对于二维灰度图像进行DCT变换,就能得到图像的频谱图:低阶(变化幅度小)的部分反映在DCT的左上方,高阶(变化幅度大)的部分反映在DCT的右下方。由于人眼对高阶部分不敏感,依靠低阶部分就能基本识别出图像内容,所以JPEG进行压缩的时候,基本上只存储DCT变换后的左上部分,而右下部分则直接丢弃。

MATALAB代码验证:

clear;
clc; im = imread('https://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png'); % 读入图像
figure(),
subplot(221),
imshow(im);
title('原始彩色图像'); grayim = rgb2gray(im);
dctim = dct2(grayim);
subplot(222),
% imshow(I,[]) displays the grayscale image I scaling the display based on the range of pixel values in I.
imshow(log(abs(dctim)), []),
title('DCT变换图像'); idctim = idct2(dctim);
subplot(223),
imshow(idctim, [])
title('DCT逆变换图像'); subplot(224),
imshow(grayim)
title('原始灰度图像');

运行结果:

![运行结果图片](二维DCT变换

参考文献

[1] 二维DCT变换:https://wuyuans.com/2012/11/dct2

[2] 余弦离散变换原理及应用:http://blog.csdn.net/shenziheng1/article/details/52965104

[3] MATLAB分析图像的离散余弦变换(DCT):http://blog.csdn.net/u013354805/article/details/52259471

[4] 图像DCT变换:https://feichashao.com/image_dct/