傅里叶变换与离散余弦变换技术报告

概述

图像处理的方法主要包括两大类，分别是在空间域的处理方法和频率域的处理方法。第一次大作业实现基于点运算的灰度处理是在空间域的处理，而这次大作业主要实现了图像从空间域到频率域的变换。频率域的图像信息对于图像的特征提取、空间频域滤波、图像恢复有着重要的价值，所以这种变换有着重要的意义。

自从傅里叶变换提出之后，其在物理学、电子类学科、信号处理、概率论等领域中得到了广泛的应用。对应于计算机离散的特性，傅里叶变换表现为离散傅里叶变换DFT(Discreate Fourier Transform)。到了1965年由J.W.库利和T.W.图基提出基2的快速傅里叶变换FFT(Fast Fourier Transform)，这种算法使得计算机计算离散傅里叶变换时复数乘法的次数显著减小，一维傅里叶变换的时间复杂度由 $O(n^2)$ 降低到 $O(n log(n))$ ，使得图像的实时处理成为可能。

在傅里叶级数展开时如果被展开的函数只有实部，即其傅里叶级数中只包含余弦项，可导出离散余弦变换DCT(Discreate Cosine Transform)，近年来由于专用集成电路设计的发展，DCT在图像编码中得到了广泛的应用，同时也成为了JPEG、MPEG等国际公认编码中的重要算法，此外在视频压缩中最常用的变换方法亦为DCT。

基于以上考虑，这次大作业主要实现二维空间域到频率域的快速傅里叶变换FFT，以及离散余弦变换DCT的算法，同上次作业相同，此次开发语言仍使用C#。

傅里叶变换与离散余弦变换技术报告概述目录一、快速傅里叶变换（FFT）1. 数学原理1. 一维傅里叶变换1. 一维连续傅里叶变换2. 一维离散傅里叶变换2.二维离散傅里叶变换（DFT）3. 一维快速傅里叶变换（FFT）1. 一维快速傅里叶变换算法的实现2. 一维快速傅里叶逆变换算法的实现4. 二维FFT及其的逆变换的思想及实现2. C#实现0. 图像的扩展1. fft2. ifft3. fft24. ifft25. fftshift3. 效果展示二、离散余弦变换（DCT）1. 数学原理1. 离散余弦变换VS.傅里叶变换2. 一维离散余弦变换3. 二维离散余弦变换4. 离散余弦变换的快速实现1. 一维DCT及其逆变换算法的快速实现2. 二维DCT及其逆变换算法的快速实现2. C#实现0. 图像的扩展1. dct2. idct3. dct24. idct23. 效果展示参考资料附录确定蝶式流程图对偶顶点的算法

一、快速傅里叶变换（FFT）

1. 数学原理

1. 一维傅里叶变换

1. 一维连续傅里叶变换

设 $f(t)$ 是一个连续的信号，或者仅有有限个间断点、极值点，满足Dirichlet条件，且满足

$\int_{-\infty}^{\infty}\mid f(t)\mid^2dt<\infty\tag{1.1}$

则傅里叶变换存在，定义为：

$F(w)=\int_{-\infty}^{\infty}f(t)e^{-jwt}dt\tag{1.2}$

且傅里叶反变换为：

$f(t)= {1 \over 2\pi}\int_{-\infty}^{\infty}F(w)e^{jwt}dw\tag{1.3}$

由欧拉公式有：

$e^{jwt}=cos(wt)+jsin(wt)\tag{1.4}$

其中 $f(t)$ 表示空间域为实函数，而 $F(w)$ 表示频率域为复函数，且由实部和虚部组成。满足：

$F(w)=R(w)+jI(w)=\mid F(w)\mid e^{j\phi(w)}\tag{1.5}$

$F(w)$ 称为信号 $f(w)$ 的概率密度函数，在信号处理中常常习惯把 $\mid F(w)\mid$ 称为幅度谱或者频谱，而把 $\phi(w)$ 称为相位谱或简称相谱，其中两种谱的导出公式为：

$\mid F(w)\mid = \sqrt{R^2(w)+I^2(w)} \;,\phi(w)=\arctan{I(w) \over R(w)}\tag{1.6}$

2. 一维离散傅里叶变换

对于计算机上处理的离散信号来说，采用离散傅里叶变换，其中 $N$ 次离散采样序列可以表示为 $\left\lbrace f(0),f(1),f(2),\cdots ,f(N-1)\right\rbrace$ ，令 $x$ 为离散实域变量，令 $u$ 为离散频域变量，则一维离散傅里叶变换（DFT）可以表示为：

$F(u)=\sum_{x=0}^{N-1}f(x)e^{-j2\pi ux \over N},\; u=0,1,\cdots,N-1\tag{1.7}$

且离散傅里叶反变换可以表示为：

$f(x)={1 \over N}\sum_{u=0}^{N-1}F(u)e^{j2\pi ux \over N},\; x=0,1,\cdots,N-1\tag{1.8}$

2.二维离散傅里叶变换（DFT）

对于二维的图像信号，其离散傅里叶变换可以表示为：

$F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi ({ux \over M}+{vy \over N})},\; u=0,1,\cdots,M-1,\; v=0,1,\cdots,N-1\tag{1.9}$

且离散傅里叶反变换可以表示为：

$f(x,y)={1 \over MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi ({ux \over M}+{vy \over N})},\; x=0,\cdots,M-1,\; y=0,\cdots,N-1\tag{1.10}$

与一维情况类似，频谱和相谱的导出公式为：

$\mid F(u,v)\mid = \sqrt{R^2(u,v)+I^2(u,v)} \;,\phi(u,v)=\arctan{I(u,v) \over R(u,v)}\tag{1.11}$

二维离散傅里叶变换的基本性质有：

可分离性
平移性
周期性
共轭对称性
旋转不变性

3. 一维快速傅里叶变换（FFT）

根据傅里叶变换的性质，进一步推导快速傅里叶变换（FFT），令 $W_N =e^{-j{2\pi \over N}}$ ,根据式(1.4)和式(1.5)，对于 $N$ 点序列的一维傅里叶变换对为：

$\left\{ \begin{array}{ll} F(u)= \sum_{x=0}^{N-1}f(x)W_N^{ux} &,\; u=0,1,\cdots,N-1 \\ f(x)={1 \over N}\sum_{u=0}^{N-1}F(u)W_N^{-ux} &,\; x=0,1,\cdots,N-1 \\ \end{array} \right.\tag{1.12}$

离散傅里叶变换有大量的重复运算，令矩阵：

$W_N=[W^{nm}]=\begin{bmatrix} W^0 & W^0 & W^0 & \cdots & W^0 \\ W^0 & W^1 & W^2 & \cdots & W^{N-1} \\ W^0 & W^2 & W^4 & \cdots & W^{2(N-1)}\\ \vdots & \vdots& \vdots & \ddots & \vdots \\ W^0 & W^{N-1} & W^{2(N-1)} & \cdots & W^{(N-1)(N-1)}\\ \end{bmatrix}\tag{1.13}$

令 $F_N=[F(0),F(1),\cdots,F(N-1)]^T$ ， $f_N=[f(0),f(1),\cdots,f(N-1)]^T$ ，则一维离散傅里叶正变换可写为矩阵形式，即：

$F_N=W_Nf_N\tag{1.14}$

由于 $W_N$ 具有周期性， $W_N$ 的因子取值有如下特点：

$\left\{ \begin{array}{ll} W^0 &=1\\ W^{N \over 2} &=-1\\ W_N^{N+r} &=W_N^r\\ W_N^{{N \over 2} + r} &=-W_N^r\\ W_{2N}^{2r} &=W_N^r\\ \end{array} \right.\tag{1.15}$

对于 $N+8$ 的DFT可以写成如下矩阵形式：

$\begin{bmatrix} F(0) \\F(1) \\ F(2) \\F(3) \\ F(4) \\F(5) \\ F(6) \\F(7) \\ \end{bmatrix} =\begin{bmatrix} 1& 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1& W^1 & W^2 & W^3&-1&-W^1&-W^2 & -W^3 \\ 1& W^2 & -1 & -W^2&1&W^2&-1 & -W^2 \\ 1& W^3 & -W^2 & W^1&-1&-W^3&W^2 & -W^1 \\ 1& -1 & 1 & -1 & 1 & -1 & 1 & -1 \\ 1& -W^1 & W^2 & -W^3&-1&W^1&-W^2 & W^3 \\ 1& -W^2 & -1 & W^2&1&-W^2&-1 & W^2 \\ 1& -W^3 & -W^2 & -W^1&-1&W^3&W^2 & W^1 \\ \end{bmatrix} \begin{bmatrix} f(0) \\f(1) \\ f(2) \\f(3) \\ f(4) \\f(5) \\ f(6) \\f(7) \\ \end{bmatrix} \tag{1.16}$

经过一系列的初等变换可得：

$\begin{bmatrix} F(0) \\F(1) \\ F(2) \\F(3) \\ F(4) \\F(5) \\ F(6) \\F(7) \\ \end{bmatrix} =\begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & -1 & W^2 & -W^2&W^1&-W^1&W^3 & -W^3 \\ 1 & 1 & -1 & -1&W^2&W^2&-W^2& -W^2 \\ 1 & -1 & -W^2 & W^2&W^3&-W^3&W^1 & -W^1 \\ 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 \\ 1 & -1 & W^2 & -W^2&-W^1&W^1&-W^3 & W^3 \\ 1 & 1 & -1 & -1&-W^2&-W^2&W^2 & W^2 \\ 1 & -1 & -W^2 & -W^2&-W^3&W^3&-W^1 & W^1 \\ \end{bmatrix} \begin{bmatrix} f(0) \\f(4) \\ f(2) \\f(6) \\ f(1) \\f(5) \\ f(3) \\f(7) \\ \end{bmatrix} \tag{1.17}$

可以画出蝶式流程图：

fft蝶式流程图

其中 $m$ 是级，表示当前算法是第几趟执行，对于具有 $N$ 个采样点的蝶式流程图具有 $M=Log_2(N)$ 级，其中 $m=0,1,\cdots,M-1$ 。例如：上图有 $8$ 个采样点，共 $M=3$ 级， $m=0,1,2$ 从左到右依次为 $m=0$ 级、 $m=1$ 级、 $m=2$ 级。

此外因子 $W$ 具有如下分布规律：

$\left. \begin{array}{} m=0 ,& W_2^r=W_N^{{N \over 2}r},\qquad &r=0\\ m=1,& W_4^r=W_N^{{N \over 4}r},\qquad &r=0,1\\ m=2,& W_8^r=W_N^{{N \over 8}r},\qquad &r=0,1,2,3\\ \cdots & \cdots\qquad \cdots\\ m=M-1 ,& W_N^r=W_N^{{N \over N}r},\qquad &r=0,1,2,\dots,{N\over 2}-1\\ \end{array} \right. \tag{1.18}$

由于 $N=2^M$ 故由上面可以看出规律为：

$\left. \begin{array}{} m=k ,& W_{2^{k+1}}^r=W_N^{{N \over 2^{k+1}}r},\qquad &r=0,1,2,\dots,2^k-1\\ \end{array} \right. \tag{1.19}$

对于上图 $N=8$ 级因子 $W$ 分布规律为：

$\left. \begin{array}{} m=0 ,& W_2^r=W_N^{4r},\qquad &r=0\\ m=1,& W_4^r=W_N^{2r},\qquad &r=0,1\\ m=2,& W_8^r=W_N^{r},\qquad &r=0,1,2,3\\ \end{array} \right. \tag{1.20}$

1. 一维快速傅里叶变换算法的实现

从面的分析可以得出以下几个结论，用于指导我们实现FFT算法。

蝶式流程图的最大级数 $M$ 与采样点个数 $N$ 的关系为 $M=Log_2(N)$
对于每一级计算的 $W$ 下标都不同，由于傅里叶级数具有周期性，故可以统一变换成为下表均为采样点个数 $N$ 的形式，所以仅需要计算 $N\over 2$ 个 $W_N^r$ ，即仅需计算 $W_N^r,r=0,1,2,\cdots,{N\over 2}-1$ 。
频域和实域的对应点映射关系发生了变化，对应点的计算由PPT上的二进制算法给出，在这个报告的附录中会做一个简要的介绍。
通过蝶形运算计算由实域变换到频域相应的值

2. 一维快速傅里叶逆变换算法的实现

由公式(1.7)和(1.8)我们发现，求FFT逆变换的算法步骤如下：

先求出 $F(u)$ 的共轭 $F^*(u)$ 代入FFT正变换的公式。得到 $Nf^*(x)$ 。
再将结果求共轭并除以系数 $N$ 即可得到所求的 $f(x)$ 。
IFFT实现了信号从频域到实域的转换。

4. 二维FFT及其的逆变换的思想及实现

FFT算法是基2的快速算法，而图像的长宽并不一定是基2的，先对图像进行扩展。
根据傅里叶变换的性质，二维的傅里叶变换可以通过两步的一维傅里叶变换实现。
可以先在 $x$ 方向上实现一维的傅里叶变换，再将结果在 $y$ 方向上实现另一维的傅里叶变换。
IFFT2的实现与上述描述类似。

2. C#实现

二维的快速傅里叶变换实现了如下的 $5$ 个函数(为了与MATLAB保持一致，采用相应的名称)，分别为一维的快速傅里叶变换及其逆变换，二维的快速傅里叶变换和其逆变换，最后一个函数实现了把低频的频谱从四个角移到中心：

public Complex[] fft(Complex[] rawArray,int N)
public Complex[] ifft(Complex[] rawArray, int N)
public Complex[] fft2(Bitmap currentBitmap)
public byte[] ifft2(Complex[] fftBuffer)
public byte[] fftshift(byte[] imgBuf)

同时实现了 $3$ 个图像绘制的相关函数，其中第一个函数把相继执行FFT2和IFFT2的图显示出来，进而与原图对比，后两个函数分别画出FFT2变换后的频率谱和相位谱：

public Bitmap GetIfft2Bitmap(byte[] imgBuffer)
public Bitmap GetAmplitudeBitmap(byte[] imgBuffer)
public Bitmap GetPhaseBitmap(byte[] imgBuffer)

0. 图像的扩展

根据前面的描述可知，我们实现的FFT算法是基 $2$ 的快速算法，而图像的长宽并不一定是基 $2$ 的所以要对图像进行扩展使得算法能够得以执行。


 
if (currentBitmap != null)
{
  int i = 1;
  //宽度扩展
  while (i < currentBitmap.Width)
    i *= 2;
  imageWidth = i;
  //高度扩展
  i = 1;
  while (i < currentBitmap.Height)
    i *= 2;
  imageHeight = i;
}

1. fft

在第一步描述的fft算法流程具体实现步骤如下：

计算最大级数


 
int stageNum=Convert.ToInt32(Math.Log(N, 2));//求fft的级数

计算 $N\over 2$ 个 $W_N^r$


 
Complex[] W = new Complex[N / 2];
//计算W(i,N) 算一半
for (int i = 0; i < N / 2; i++)
{
  double angle = -i * 2.0 * Math.PI / N;
  W[i] = new Complex((float)Math.Cos(angle), (float)Math.Sin(angle));
}

计算频域和实域的对应点（算法讨论见附录）


 
Complex[] tmpBuffer1 = new Complex[N];
Complex[] tmpBuffer2 = new Complex[N];
//复制待处理数组(浅拷贝)
rawArray.CopyTo(tmpBuffer1, 0);
//确定数组元素的对偶节点并重新排序
for (int i = 0; i < N; i++)
{
  int originalPos = i;
  int postPos = 0;
  for (int j = 0; j < stageNum; j++)
  {
    //逆向添加每一位，postPos初始为0
    postPos <<= 1;
    //bit test 取得末位，并添加到postPos的末位
    postPos |= (originalPos - ((originalPos >> 1) << 1));
    //正向遍历每一位
    originalPos >>= 1;
  }
  tmpBuffer2[postPos] = tmpBuffer1[i];
}
tmpBuffer2.CopyTo(tmpBuffer1, 0);

蝶形运算计算由实域变换到频域相应的值


 
//蝶形运算
for (int stage = 0; stage < stageNum; stage++)
{
  int interval = 1 << stage;                      //组间间隔
  int groupNum = 1 << (stageNum - stage - 1);     //总的组数
  int groupSize = N / groupNum;                   //组的大小
  for (int group = 0; group < groupNum; group++)
  {
    int groupStart = groupSize * group;         //本组前半部分的第一个元素
    int groupMid = groupStart + groupSize / 2;  //本组后半部分的第一个元素
    int groupTail = groupStart + groupSize;     //本组的上界
    for (int r = 0; r < interval; r++)          //本组后半部分每一个元素是上次输出元素的W(groupNum*r,N)倍（复数乘积）
    {
      tmpBuffer1[groupMid + r] = tmpBuffer1[groupMid + r] * W[groupNum * r];
    }
    for (int k = groupStart; k < groupMid; k++) //本组前半部分为复数加
    {
      tmpBuffer2[k] = tmpBuffer1[k] + tmpBuffer1[k + groupSize / 2];
    }
    for (int k = groupMid; k < groupTail; k++)  //本组后半部分为复数减
    {
      tmpBuffer2[k] = tmpBuffer1[k - groupSize / 2] - tmpBuffer1[k];
    }
  }
  tmpBuffer2.CopyTo(tmpBuffer1, 0);
}

2. ifft

先求出 $F(u)$ 的共轭 $F^*(u)$ 并代入FFT正变换的公式。得到 $Nf^*(x)$ 。


 
Complex[] tmpBuffer = new Complex[N];
rawArray.CopyTo(tmpBuffer, 0);
//求F(u)的共轭
for (int i = 0; i < N; i++) 
{
  tmpBuffer[i].Imaginary = -tmpBuffer[i].Imaginary;
}
//调用快速傅里叶变换
tmpBuffer = fft(tmpBuffer, N);

再将结果( $Nf^*(x)$ )求共轭并除以系数 $N$ 即可得到所求的 $f(x)$ 。


 
//求结果的共轭并除以系数N
for (int i = 0; i < N; i++)
{
  tmpBuffer[i].Real = tmpBuffer[i].Real / N;
  tmpBuffer[i].Imaginary = -tmpBuffer[i].Imaginary / N;
}

3. fft2

先在 $y$ 方向上实现一维的快速傅里叶变换


 
Complex[] tmpBuffer1;
Complex[] tmpBuffer2;
tmpBuffer1 = new Complex[imageHeight];
tmpBuffer2 = new Complex[imageHeight];
for (int u = 0; u < imageWidth; u++)
{
  for (int v = 0; v < imageHeight; v++)
  {
    tmpBuffer1[v] = fftBuffer[v * imageWidth + u];
  }
  tmpBuffer2 = fft(tmpBuffer1, imageHeight);
  for (int v = 0; v < imageHeight; v++)
  {
    fftBuffer[v * imageWidth + u] = tmpBuffer2[v];
  }
}

再在 $x$ 方向上实现一维的快速傅里叶变换


 
tmpBuffer1 = new Complex[imageWidth];
tmpBuffer2 = new Complex[imageWidth];
for (int u = 0; u < imageHeight; u++)
{
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    tmpBuffer1[v] = fftBuffer[u * imageWidth + v];
  }
  tmpBuffer2 = fft(tmpBuffer1, imageWidth);
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    fftBuffer[u * imageWidth + v] = tmpBuffer2[v];
  }
}

计算幅度谱


 
//将幅度谱图以图像形式存入imgAmplitudeBuffer
double t;
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    t = fftBuffer[i * imageWidth + j].Abs();
    t = t / 500;
    if (t > 255)
      imgAmplitudeBuffer[i * imageWidth + j] = 255;
    else
      imgAmplitudeBuffer[i * imageWidth + j] = (byte)(t + 0.5);
  }
}

计算相位谱


 
//将相位谱图以图像形式存入imgBuffer
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    t = fftBuffer[i * imageWidth + j].Angle() + 2 * Math.PI;  //0≤t≤4π
    t = t * 255 / (4 * Math.PI);
    if (t > 255)
      imgPhaseBuffer[i * imageWidth + j] = 255;
    else
      imgPhaseBuffer[i * imageWidth + j] = (byte)(t + 0.5);
  }
}

4. ifft2

先在 $x$ 方向上实现一维的快速傅里叶逆变换


 
Complex[] tmpBuffer1;
Complex[] tmpBuffer2;
tmpBuffer1 = new Complex[imageHeight];
tmpBuffer2 = new Complex[imageHeight];
for (int u = 0; u < imageWidth; u++)
{
  for (int v = 0; v < imageHeight; v++)
  {
    tmpBuffer1[v] = fftBuffer[v * imageWidth + u];
  }
  tmpBuffer2 = ifft(tmpBuffer1, imageHeight);
  for (int v = 0; v < imageHeight; v++)
  {
    fftBuffer[v * imageWidth + u] = tmpBuffer2[v];
  }
}

再在 $y$ 方向上实现一维的快速傅里叶逆变换


 
tmpBuffer1 = new Complex[imageWidth];
tmpBuffer2 = new Complex[imageWidth];
for (int u = 0; u < imageHeight; u++)
{
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    tmpBuffer1[v] = fftBuffer[u * imageWidth + v];
  }
  tmpBuffer2 = ifft(tmpBuffer1, imageWidth);
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    fftBuffer[u * imageWidth + v] = tmpBuffer2[v];
  }
}

提取恢复后的图像


 
byte[] imgBuffer = new byte[imageWidth * imageHeight];
//将变换后的图像存入imgBuffer
double t;
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    t = fftBuffer[i * imageWidth + j].Abs();
    if (t > 255)
      imgBuffer[i * imageWidth + j] = 255;
    else
      imgBuffer[i * imageWidth + j] = (byte)(t + 0.5);
  }
}

5. fftshift

把低频的频谱从四个角移到中心


 
int i0, j0;
byte[] tmpBuffer = new byte[imgBuf.Length];
imgBuf.CopyTo(tmpBuffer, 0);
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    if (i < imageHeight / 2)
      i0 = i + imageHeight / 2;
    else
      i0 = i - imageHeight / 2;
    if (j < imageWidth / 2)
      j0 = j + imageWidth / 2;
    else
      j0 = j - imageWidth / 2;
    imgBuf[i * imageWidth + j] = tmpBuffer[i0 * imageWidth + j0];
  }
}

3. 效果展示

主要的绘制代码如下：


 
if (firstPainted)
{
  fftBuffer = fft2(currentBitmap);
  //对fftBuffer进行深拷贝
  for (int i = 0; i < fftBuffer.Length; i++)
  {
    fftTmpBuffer[i] = new Complex(fftBuffer[i]);
  }
  currentAmplitudeBitmap = GetAmplitudeBitmap(fftshift(imgAmplitudeBuffer));
  currentPhaseBitmap = GetPhaseBitmap(fftshift(imgPhaseBuffer));
  currentIfft2Bitmap = GetIfft2Bitmap(ifft2(fftTmpBuffer));
  firstPainted = false;
}

说明

由于对图像进行了扩展，所以幅度谱和相位谱的大小跟原来的图略有差别，在执行完ifft2之后得到的图像虚部都非常的小，取原图像大小的实部进行图像绘制，结果如下： FFT

二、离散余弦变换（DCT）

1. 数学原理

1. 离散余弦变换VS.傅里叶变换

在图像压缩领域中，傅里叶变换的一个不足之处在于，它的参数都是复数，在数据的描述上相当于实数的两倍，而且不易计算，因此，人们希望有一种计算的数据量不太大而且压缩效果又很好的变换，在此思想的推动下，出现了离散余弦变换。离散余弦变换与傅里叶变换相比仅采用实部而没有虚部。

2. 一维离散余弦变换

实域函数 $f(x)$ 的一维离散余弦变换定义为：

$C(u)=a(u)\sum_{x=0}^{N-1}f(x)cos{(2x+1)u\pi \over 2N},u=0,1,2,\cdots,N-1\tag{2.1}$

一维离散余弦变换的反变换定义为：

$f(x)=\sum_{u=0}^{N-1}a(u)C(u)cos{(2x+1)u\pi \over 2N},x=0,1,2,\cdots,N-1\tag{2.2}$

其中， $a(u)$ 为归一化加权系数，由下式定义：

$a(u)=\left\lbrace \begin{array}{} \sqrt{1\over N},\qquad u=0\\ \sqrt{2\over N},\qquad u=1,2,3,\cdots,N-1\\ \end{array} \right.\tag{2.3}$

3. 二维离散余弦变换

将一维离散余弦变换扩展到二维离散余弦变换，其定义为：

$C(u,v)=a(u)a(v)\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}f(x,y)cos[{(2x+1)u\pi \over 2N}]cos[{(2y+1)v\pi \over 2N}],x,y=0,1,\dots,N-1\tag{2.4}$

二维离散余弦变换的反变换定义为：

$f(x,y)=\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}a(u)a(v)C(u,v)cos[{(2x+1)u\pi \over 2N}]cos[{(2y+1)v\pi \over 2N}],u,v=0,1,\dots,N-1\tag{2.5}$

二维离散余弦变换是可分离的，因此二维的离散余弦变换及其逆变换可以逐次应用一维DCT算法加以计算。

4. 离散余弦变换的快速实现

1. 一维DCT及其逆变换算法的快速实现

DCT算法可以通过FFT来实现，下面通过对一维离散余弦变换公式的推导得出其与傅里叶变换之间的关系。

由公式(2.1)可得：

$\begin{array}{} F(u) &=a(u)\sum_{x=0}^{N-1}f(x)cos{(2x+1)u\pi \over 2N}\\ &=a(u)Re\{\sum_{x=0}^{N-1}f(x)e^{-j{(2x+1)u\pi}\over 2N}\}\\ \end{array} \tag{2.6}$

对实域数据作如下延拓：

$fe(x)=\left \lbrace\begin{array}{}f(x),\quad &x=0,1,2,\cdots,N-1\\0, &x=N,N+1,\cdots,2N-1\end{array} \right.\tag{2.7}$

则 $fe(x)$ 的离散余弦变换可以写为：

$\begin{array}{} F(u)=Fe(u) &=a(u)\sum_{x=0}^{2N-1}fe(x)cos{(2x+1)u\pi \over 2N}\\ &=a(u)Re\{\sum_{x=0}^{2N-1}fe(x)e^{-j{(2x+1)u\pi}\over 2N}\}\\ &=a(u)Re\{e^{-j{u\pi\over 2N}}\sum_{x=0}^{2N-1}fe(x)e^{-j{2xu\pi \over 2N}}\}\\ \end{array} \tag{2.8}$

对比式(1.7)和(2.8)发现 $\sum_{x=0}^{2N-1}fe(x)e^{-j{2xu\pi \over 2N}}$ 为 $2N$ 点的离散傅里叶变换。

故一维离散余弦变换的算法可以梳理如下：

将实域数据按照公式(2.7)延拓至长度为 $2N$ 点的实域数据。
计算 $2N$ 点的离散傅里叶变换，采用FFT来实现，每个结果与 $W=e^{-j{\pi \over 2N}}$ 做复数乘法。
取上一步计算结果的实部并乘以归一化加权系数 $a(u)$ 。

同理对于离散余弦变换IDCT 先在频率空间将 $Fe(u)$ 做如下延拓：

$Fe(u)=\left \lbrace\begin{array}{}F(u),\quad &u=0,1,2,\cdots,N-1\\0, &u=N,N+1,\cdots,2N-1\end{array} \right.\tag{2.9}$

则离散余弦反变换可以写为：

$\begin{array}{} f(x)=fe(x) &=\sum_{u=0}^{2N-1}a(u)Fe(u)cos{(2x+1)u\pi \over 2N}\\\\ &=\sqrt{1\over N}Fe(0)+\sqrt{2\over N}\sum_{u=1}^{2N-1}Fe(u)cos{(2x+1)u\pi \over 2N}\\ &=\sqrt{1\over N}Fe(0)+\sqrt{2\over N}Re\{\sum_{u=1}^{2N-1}Fe(u)e^{j{{2xu\pi}\over 2N}}e^{j{u\pi \over 2N}}\}\\ &=(\sqrt{1\over N}-\sqrt{2\over N})Fe(0)+\sqrt{2\over N}Re\{\sum_{u=0}^{2N-1}[Fe(u)e^{j{u\pi \over 2N}}]e^{j{2xu\pi\over 2N}}\}\\ &=(\sqrt{1\over N}-\sqrt{2\over N})Fe(0)+2N\sqrt{2\over N}Re\{{1\over 2N}\sum_{u=0}^{2N-1}[Fe(u)e^{j{u\pi \over 2N}}]e^{j{2xu\pi\over 2N}}\}\\ \end{array} \tag{2.10}$

对比式(1.8)和(2.10)发现 ${1\over 2N}\sum_{u=0}^{2N-1}[Fe(u)e^{j{u\pi \over 2N}}]e^{j{2xu\pi\over 2N}}$ 为函数 $[Fe(u)e^{j{u\pi \over 2N}}]$ 在 $2N$ 点上的离散傅里叶逆变换。

故一维离散余弦逆变换的算法可以梳理如下：

将频域数据按照公式(2.9)延拓至长度为 $2N$ 点的频域数据。
将原始数据与 $W_u=e^{j{u\pi \over 2N}}$ 做复数乘法，得到中间函数 $Fe(u)W_u$ 。
计算 $2N$ 点的离散傅里叶逆变换，采用IFFT来实现。
取上一步计算结果的实部并乘以 $2N\sqrt{2\over N}$ ，之后与常数 $d_0=(\sqrt{1\over N}-\sqrt{2\over N})Fe(0)$ 相加得到最终结果。

2. 二维DCT及其逆变换算法的快速实现

DCT算法是基2的算法，而图像的长宽并不一定是基2的，先对图像进行扩展。
根据离散余弦变换的性质，二维的离散余弦变换可以通过两步的一维离散余弦变换实现。
可以现在 $x$ 方向上实现一维的离散余弦变换，再将结果在 $y$ 方向上实现另一维的离散余弦变换。
IDCT2的实现与上述描述类似。

2. C#实现

二维的离散余弦变换实现了如下的 $6$ 个函数(为了与MATLAB保持一致，采用相应的名称)，分别为一维的快速傅里叶变换及其逆变换(与上述fft略有区别，这一fft在函数中对于数据进行了延拓，ifft类似进行了延拓，具体实现之后会提到)，一维的离散余弦变换及其逆变换，以及二维的离散余弦变换及其逆变换：

public Complex[] fft(Complex[] rawArray,int N)
public Complex[] ifft(Complex[] rawArray, int N)
public double[] dct(Complex[] rawArray, int N)
public double[] idct(Complex[] rawArray, int N)
public double[] dct2(Bitmap currentBitmap)
public byte[] idct2(double[] dctBuffer)

同时实现了 $2$ 个图像绘制的相关函数，其中前一个函数把相继执行了DCT2和IDCT2的图显示出来，进而与原图对比，后一个函数画出DCT2变换后的频率谱：

public Bitmap GetIdct2Bitmap(byte[] imgBuffer)
public Bitmap GetAmplitudeBitmap(byte[] imgBuffer)

0. 图像的扩展

根据前面的描述可知，我们实现的DCT算法是基 $2$ 的算法，而图像的长宽并不一定是基 $2$ 的所以要对图像进行扩展使得算法能够得以执行。


 
//判断图像是否需要扩展
if (needExpand())
{
  isExpanded = true;
  int i = 1;
  //宽度扩展
  while (i < currentBitmap.Width)
    i *= 2;
  imageWidth = i;
  //高度扩展
  i = 1;
  while (i < currentBitmap.Height)
    i *= 2;
  imageHeight = i;
}
else
{
  isExpanded = false;
  imageWidth = currentBitmap.Width;
  imageHeight = currentBitmap.Height;
}

1. dct

将实域数据按照公式(2.7)延拓至长度为 $2N$ 点的实域数据(在fft内实现)。


 
rawArray.CopyTo(tmpBuffer1,0);
for (int i = N / 2; i < N; i++)
{
  tmpBuffer1[i] = new Complex(0, 0);
}

计算 $2N$ 点的离散傅里叶变换，采用FFT来实现，每个结果与 $W=e^{-j{\pi \over 2N}}$ 做复数乘法，将计算结果的实部乘以归一化加权系数 $a(u)$ 。


 
Complex[] fftBufer = null;
double[] dctBuffer = new double[N];
double au = Math.Sqrt(1.0 / N);
fftBufer = fft(rawArray, 2 * N);
dctBuffer[0] = au * fftBufer[0].Real;
au *= Math.Sqrt(2.0);
for (int i = 1; i < N; i++)
{
  double angle = -i * Math.PI / N;
  Complex W = new Complex((float)Math.Cos(angle), (float)Math.Sin(angle));
  dctBuffer[i] = au * (W * fftBufer[i]).Real;
}

2. idct

将频域数据按照公式(2.9)延拓至长度为 $2N$ 点的频域数据(在ifft中实现)。


 
rawArray.CopyTo(tmpBuffer1,0);
for (int i = N / 2; i < N; i++)
{
  tmpBuffer1[i] = new Complex(0, 0);
}

将原始数据与 $W_u=e^{j{u\pi \over 2N}}$ 做复数乘法，得到中间函数 $Fe(u)W_u$ 。


 
double tmpNum = Math.Sqrt(2.0 / N);
for (int i = 0; i < N; i++)
{
  double angle = i * Math.PI / (2.0 * N);
  Complex W = new Complex((float)Math.Cos(angle), (float)Math.Sin(angle));
  tmpBuffer[i] = rawArray[i] * W;
}

计算 $2N$ 点的离散傅里叶逆变换，采用IFFT来实现。


 
ifftBufer = ifft(tmpBuffer, 2 * N);

取上一步计算结果的实部并乘以 $2N\sqrt{2\over N}$ ，之后与常数 $d_0=(\sqrt{1\over N}-\sqrt{2\over N})Fe(0)$ 相加得到最终结果。


 
double d0 = (Math.Sqrt(1.0 / N) - Math.Sqrt(2.0 / N)) * ifftBufer[0].Real;
for (int i = 1; i < N; i++)
{
  dctBuffer[i] = d0 + tmpNum * ifftBufer[i].Real * 2.0 * N;
}

3. dct2

先在 $y$ 方向上实现一维的离散余弦变换


 
Complex[] tmpBuffer1;
double[] tmpBuffer2;
tmpBuffer1 = new Complex[imageHeight];
tmpBuffer2 = new double[imageHeight];
for (int u = 0; u < imageWidth; u++)
{
  for (int v = 0; v < imageHeight; v++)
  {
    tmpBuffer1[v] = new Complex((float)dctBuffer[v * imageWidth + u],0);
  }
  tmpBuffer2 = dct(tmpBuffer1, imageHeight);
  for (int v = 0; v < imageHeight; v++)
  {
    dctBuffer[v * imageWidth + u] = tmpBuffer2[v];
  }
}

再在 $x$ 方向上实现一维的离散余弦变换


 
tmpBuffer1 = new Complex[imageWidth];
tmpBuffer2 = new double[imageWidth];
for (int u = 0; u < imageHeight; u++)
{
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    tmpBuffer1[v] = new Complex((float)dctBuffer[u * imageWidth + v],0);
  }
  tmpBuffer2 = dct(tmpBuffer1, imageWidth);
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    dctBuffer[u * imageWidth + v] = tmpBuffer2[v];
  }
}

计算幅度谱


 
//将幅度谱图以图像形式存入imgAmplitudeBuffer
double t;
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    t = Math.Abs(dctBuffer[i * imageWidth + j]);
    if (t > 255)
      imgAmplitudeBuffer[i * imageWidth + j] = 255;
    else
      imgAmplitudeBuffer[i * imageWidth + j] = (byte)(t + 0.5);
  }
}

4. idct2

先在 $x$ 方向上实现一维的离散余弦逆变换


 
tmpBuffer1 = new Complex[imageHeight];
tmpBuffer2 = new double[imageHeight];
for (int u = 0; u < imageWidth; u++)
{
  for (int v = 0; v < imageHeight; v++)
  {
    tmpBuffer1[v] = new Complex((float)dctBuffer[v * imageWidth + u],0);
  }
  tmpBuffer2 = idct(tmpBuffer1, imageHeight);
  for (int v = 0; v < imageHeight; v++)
  {
    dctBuffer[v * imageWidth + u] = tmpBuffer2[v];
  }
}

再在 $y$ 方向上实现一维的离散余弦逆变换


 
tmpBuffer1 = new Complex[imageWidth];
tmpBuffer2 = new double[imageWidth];
for (int u = 0; u < imageHeight; u++)
{
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    tmpBuffer1[v] = new Complex((float)dctBuffer[u * imageWidth + v],0);
  }
  tmpBuffer2 = idct(tmpBuffer1, imageWidth);
  for (int v = 0; v < imageWidth; v++)//性能瓶颈
  {
    dctBuffer[u * imageWidth + v] = tmpBuffer2[v];
  }
}

提取恢复后的图像


 
byte[] imgBuffer = new byte[imageWidth * imageHeight];
//将幅度谱图以图像形式存入imgBuffer
double t;
for (int i = 0; i < imageHeight; i++)
{
  for (int j = 0; j < imageWidth; j++)
  {
    t = dctBuffer[i * imageWidth + j]*255.0/400.0;
    if (t > 255)
      imgBuffer[i * imageWidth + j] = 255;
    else
      imgBuffer[i * imageWidth + j] = (byte)(t + 0.5);
  }
}

3. 效果展示

主要的绘制代码如下：


 
if (firstPainted)//防止窗口resize时候重新计算
{
  dctBuffer = dct2(currentBitmap);
  //对dctBuffer进行深拷贝
  dctTmpBuffer = (double[])dctBuffer.Clone();
  currentAmplitudeBitmap = GetAmplitudeBitmap(imgAmplitudeBuffer);
  currentIdct2Bitmap = GetIdct2Bitmap(idct2(dctTmpBuffer));
  firstPainted = false;
}

说明

由于对图像进行了扩展，所以幅度谱的大小跟原来的图略有差别，结果如下：

参考资料

Visual C++数字图像处理谢凤英
C#数字图像处理算法典型实例
图像的DCT变换
图像变换2离散余弦变换

附录

确定蝶式流程图对偶顶点的算法

将节点 $X_I(k)$ 的 $k$ 变为二进制数
$X_I(k)=X_I(k_{r-1}k_{r-2}\cdots k_0)$
将二进制数按比特位倒转
$X_I(k')=X_I(k_0k_1\cdots k_{r-1})$
将二进制数变换为十进制数，即为重新排序的位置，例：
$\begin{array}{} X_3(0)&\to X_3(000)&\to X_3(000)&\to X_3(0)\\ X_3(1)&\to X_3(001)&\to X_3(100)&\to X_3(4)\\ X_3(2)&\to X_3(010)&\to X_3(010)&\to X_3(2)\\ X_3(3)&\to X_3(011)&\to X_3(110)&\to X_3(6)\\ X_3(4)&\to X_3(100)&\to X_3(001)&\to X_3(1)\\ X_3(5)&\to X_3(101)&\to X_3(101)&\to X_3(5)\\ X_3(6)&\to X_3(110)&\to X_3(011)&\to X_3(3)\\ X_3(7)&\to X_3(111)&\to X_3(111)&\to X_3(7)\\ \end{array}$