DCT算法的原理和优化

未经博主容许,禁止转载!web

离散余弦变换(DCT for Discrete Cosine Transform)是与傅里叶变换相关的一种变换,它与离散傅里叶变换相似,可是只使用实数。算法

这种变化常常被信号处理和图像处理使用,用于对信号和图像(包括静止图像和运动图像)进行有损压缩。在压缩算法中,现将输入图像划分为8*8或16*16的图像块,对每一个图像块做DCT变换;而后舍弃高频的系数,并对余下的系数进行量化以进一步减小数据量;最后使用无失真编码来完成压缩任务。解压缩时首先对每一个图像块做DCT反变换,而后将图像拼接成一副完整的图像。ubuntu

大多数天然信号(包括声音和图像)的能量都集中在余弦变换后的低频部分。因为人眼对于细节信息不是很敏感,所以信息含量更少的高频部分能够直接去掉,从而在后续的压缩操做中得到较高的压缩比。并发

  • 1. 一维DCT变换
  • 2. 二维DCT变换
  • 3. DCT反变换
  • 4. 怎样提升二维DCT变换的速度

1. 一维DCT变换

一维DCT变换共有8种,最实用的一种公式以下:
这里写图片描述svg

其中N是一维数据的元素总数,c(u)系数使得DCT变换矩阵成为正交矩阵,正交特性在二维DCT变换中更能体现其优点。一维DCT变换的复杂度是O(n^2)工具

2. 二维DCT变换

二维DCT变换公式以下:性能

这里写图片描述

咱们将公式变换一下:
这里写图片描述ui

能够发现,二维DCT变换实际上是在一维DCT变换的基础上,再作一次一维DCT变换。二维DCT也能够写成矩阵相乘的形式:编码

这里写图片描述

矩阵A和f必须为长宽均为N的方阵,f不是方阵的状况能够补零再对其作DCT变换。生成矩阵A的matlab代码以下:spa

A=zeros(N);
 for i=0:N
    for j=0:N
        if i==0
            a=sqrt(1/N);
        else
            a=sqrt(2/N);
        end            
        A(i+1,j+1)=a*cos(pi*(2*j+1)*i/(2*N));
    end
end

二维DCT变换的复杂度达到O(n^4),因此进行DCT变换的矩阵不宜过大。在实际处理图片的过程当中,须要先把矩阵分块,通常分为8*8或16*16大小,这样DCT变换不至于耗费过多的时间。

3. DCT反变换

DCT变换是无损并可逆的,反变换的公式以下:

这里写图片描述

一样,DCT反变换也能够写成矩阵相乘的形式:

这里写图片描述

DCT反变换的时间复杂度与DCT变换相同。

4. 怎样提升DCT变换的速度

因为DCT变换的高效,这种变换方法已经普遍应用于信号处理和图像处理,成为许多图像编码国际标准的核心。但随着信息时代的来临和信息量的飞速增加,人们对图像处理的速度要求愈来愈高,虽然离散余弦变换已经比较高效了,可是可否在DCT变换的基础上,继续提高其速度呢?

咱们考虑用并行计算的方法提高速度。并行计算指的是同时使用多种计算资源解决计算问题的一种高效的计算手段。图形处理器(Graphics Processing Unit,缩写GPU)很是擅长大规模并发计算,因此我考虑在GPU上实现并行的二维DCT变换,看可否有效提高二维DCT变换的速度。

咱们观察二维DCT变换的公式,不难发现每个DCT系数的计算,不依赖于任何一个DCT系数,也就是说,任意多个DCT系数的计算,均可以独立进行。因此我考虑对整个矩阵分线程块计算,线程块内再以每个系数为单位分为线程,每个DCT系数由一个线程来完成。记矩阵的大小为length*length,块大小为block_len* block_len。

DCT算法的CUDA代码以下,环境为ubuntu 14.04,CUDA版本8.0:

int tidy = blockIdx.x * blockDim.x + threadIdx.x;
    int tidx = blockIdx.y * blockDim.y + threadIdx.y;
    int index = tidx * length + tidy;
    int i;
    float tmp;
    float beta, alfa;
    if (tidx == 0)
        beta = sqrt(1.0 / length);
    else
        beta = sqrt(2.0 / length);
    if (tidy == 0)
        alfa = sqrt(1.0 / length);
    else
        alfa = sqrt(2.0 / length);
    if (tidx < length && tidy < length) {
        for (i = 0; i < length * length; i++) {
            int x = i / length;
            int y = i % length;
            tmp += ((double) f[i])
                    * cos((2 * x + 1) * tidx * PI / (2.0 * length))
                    * cos((2 * y + 1) * tidy * PI / (2.0 * length));
        }
        F[index] = (double) alfa * beta * tmp;
    }
    __syncthreads();
}

因为GPU的控制单元能够把多个访问合并成少的访问,同时,计算每一个系数须要访问的数据都彻底相同,致使同一个块内的线程对于原矩阵的访问能够达到很是高的效率,因此我直接将块数设置为1,让全部计算系数的线程处于一个块内。

经过实验发现,基于CPU实现的DCT算法随着矩阵的增大,处理时间急剧的上升,而基于GPU实现的并行二维DCT变换算法,耗时很是短,计算速度惊人。

这里写图片描述

Matrix Length DCT Based on GPU(ms) DCT Based on CPU(ms) speedup
8 2.111 0.111 0.053
16 5.503 0.336 0.061
32 35.1296 2.844 0.081
64 1.105 16.338 14.785
128 2.312 218.512 94.291
256 2.151 1534.86 213.556
320 1.961 3034.05 1547.195

从上表中咱们发现,当矩阵比较小的时候,GPU并行计算没有带来优点,反而让DCT变换时间更长了。经过CUDA提供的程序性能分析工具,观察到GPU并行程序的大部分时间不是花在了计算上面,而是CUDA任务的建立上。当矩阵增大以后,CUDA程序的运行时间再也不主要取决于Event的建立。而基于CPU的DCT变换的时间基原本自于计算,因此随着矩阵的增大,时间基本呈指数增加。

这里写图片描述

经过以上内容,我完整的了解了一维和二维DCT变换的详细过程,而且经过GPU并行计算为二维DCT变换加速,实验结果显示,对于比较大的矩阵,基于GPU的二维DCT变换速度获得了巨大的提高,可是在矩阵比较小的应用场景下,CPU更有优点。