DTFT,DFT,DCT的关系
DTFT是离散时间傅里叶变换,用来表达连续的信号的频谱。
然后理解DFT:
DFT是离散傅里叶变换,针对的是离散的信号和频谱。DFT是DTFT变化而来,其实就是将连续时间t变成了nT. 为什么要这样做呢,因为计算机是在数字环境下工作的,它不可能看见或者处理现实中连续的信号,只能够进行离散计算,在真实性上尽可能地逼近连续信号。所以DFT是为了我们能够去用工具分析信号而创造出来的,通常我们直接用DTFT的机会很少。
然后再理解FFT:
首先,DCT是DFT的一种形式。所谓“余弦变换”,是在DTFT傅立叶级数展开式中,如果被展开的函数是实偶函数,那么其傅立叶级数中只包含余弦项,再将其离散化(DFT)可导出余弦变换,因此称之为离散余弦变换(DCT)。其实DCT属于DFT的一个子集。DCT用于语音和图像处理比较多。
当u,v = 0 时,离散余弦正变换(DCT)后的系数若为F(0,0)=1,则离散余弦反变换(IDCT)后的重现函数 f(x,y)=1/8,是个常 数值,所以将 F(0,0) 称为直流(DC)系数;当 u,v≠0 时,正变换后的系数为 F(u,v)=0,则反变换后的重现函数 f(x,y) 不是常数,此时 正变换后的系数 F(u,v) 为交流(AC)系数。
矩阵表示法
更为简洁的定义方法是采用矩阵式定义。根据以上公式定义可知,离散余弦变换的系数矩阵可以写成如下:
如果令N=4,那么由一维解析式定义可得如下展开式。
写成矩阵式
二维变换的函数如下:
离散余弦变换(DCT)公式:
离散余弦逆变换(IDCT)公式:
二维离散余弦变换的正变换公式为:
在图像的压缩编码中,N一般取8.
二维DCT的反变换公式为:
以上各式中的系数:
也就是:
代码实现:
#include <iostream>
#include <math.h>
using namespace std;
#define NUM 8
#define PI 3.1415926
short round(double a)
{
if (a >= 0)
{
return (short)(a + 0.5);
}
else
{
return (short)(a - 0.5);
}
}
// DCT - Discrete Cosine Transform
void DCT(short data[NUM][NUM])
{
short output[NUM][NUM];
double ALPHA, BETA;
short u = 0;
short v = 0;
short i = 0;
short j = 0;
for(u = 0; u < NUM; u++)
{
for(v = 0; v < NUM; v++)
{
if(u == 0)
{
ALPHA = sqrt(1.0 / NUM);
}
else
{
ALPHA = sqrt(2.0 / NUM);
}
if(v == 0)
{
BETA = sqrt(1.0 / NUM);
}
else
{
BETA = sqrt(2.0 / NUM);
}
double tmp = 0.0;
for(i = 0; i < NUM; i++)
{
for(j = 0; j < NUM; j++)
{
tmp += data[i][j] * cos((2*i+1)*u*PI/(2.0 * NUM)) * cos((2*j+1)*v*PI/(2.0 * NUM));
}
}
output[u][v] = round(ALPHA * BETA * tmp);
}
}
memset(data, 0, NUM * NUM * sizeof(short));
memcpy(data, output, NUM * NUM * sizeof(short));
}
// Inverse DCT
void IDCT(short data[NUM][NUM])
{
short output[NUM][NUM];
double ALPHA, BETA;
short u = 0;
short v = 0;
short i = 0;
short j = 0;
for(i = 0; i < NUM; i++)
{
for(short j = 0; j < NUM; j++)
{
double tmp = 0.0;
for(short u = 0; u < NUM; u++)
{
for(v = 0; v < NUM; v++)
{
if(u == 0)
{
ALPHA = sqrt(1.0 / NUM);
}
else
{
ALPHA = sqrt(2.0 / NUM);
}
if(v == 0)
{
BETA = sqrt(1.0 / NUM);
}
else
{
BETA = sqrt(2.0 / NUM);
}
tmp += ALPHA * BETA * data[u][v] * cos((2*i+1)*u*PI/(2.0 * NUM)) * cos((2*j+1)*v*PI/(2.0 * NUM));
}
}
output[i][j] = round(tmp);
}
}
memset(data, 0, NUM * NUM * sizeof(short));
memcpy(data, output, NUM * NUM * sizeof(short));
}
int main(void)
{
short i = 0;
short j = 0;
short u = 0;
short v = 0;
// 8 x 8 的图像数据
short input[NUM][NUM] =
{
{89, 101, 114, 125, 126, 115, 105, 96},
{97, 115, 131, 147, 149, 135, 123, 113},
{114, 134, 159, 178, 175, 164, 149, 137},
{121, 143, 177, 196, 201, 189, 165, 150},
{119, 141, 175, 201, 207, 186, 162, 144},
{107, 130, 165, 189, 192, 171, 144, 125},
{97, 119, 149, 171, 172, 145, 117, 96},
{88, 107, 136, 156, 155, 129, 97, 75}
};
DCT(input);
cout << "The result of DCT:" << endl;
for(u = 0; u < NUM; u++)
{
for(v = 0; v < NUM; v++)
{
cout << input[u][v] << '\t';
}
cout << endl;
}
IDCT(input);
cout << "The result of IDCT:" << endl;
for(i = 0; i < NUM; i++)
{
for(short j = 0; j < NUM; j++)
{
cout << input[i][j] << '\t';
}
cout << endl;
}
return 0;
}
运算结果如下:
The result of DCT:
1125 -32 -185 -7 2 -1 -2 2
-22 -16 45 -3 -2 0 -2 -2
-165 32 17 2 1 -1 -3 0
-7 -4 0 2 2 -1 -1 2
-2 0 0 3 0 0 2 1
3 1 1 -1 -2 0 2 0
0 0 2 -1 -1 2 1 -1
0 3 1 -1 2 1 -2 0
The result of IDCT:
89 101 114 125 126 115 105 96
97 115 131 147 149 135 123 113
114 134 159 178 175 164 149 137
121 143 177 196 201 189 165 150
119 141 175 201 207 186 162 144
107 130 165 189 192 171 144 125
96 119 150 171 172 145 116 96
88 107 136 156 155 129 97 75
另一个例子:
#include "highgui.h"
#include <math.h>
#include <cv.h>
#include "cxcore.h"
#define cvCvtPlaneToPix cvMerge
double PSNR_B = 0;
double PSNR_G = 0;
double PSNR_R = 0;
double PSNR;
int main(int argc, char* argv[])
{
const char* imagename = "F:/c/59.jpg";
IplImage *src;
CvScalar SrcPixel;
CvScalar DstPixel;
double sumB = 0;
double sumG = 0;
double sumR = 0;
double mseB;
double mseG;
double mseR;
src = cvLoadImage(imagename, 1);
if (!src)
{
printf("can't open the image...\n");
return -1;
}
// YUV颜色空间
IplImage* YUVImage = cvCreateImage(cvSize(src->width, src->height), src->depth, 3);
IplImage* dst = cvCreateImage(cvSize(src->width, src->height), src->depth, 3);
// YUV颜色空间各通道
IplImage* Y = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);
IplImage* U = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);
IplImage* V = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);
//cvNamedWindow( "Origin Image", CV_WINDOW_AUTOSIZE );
cvCvtColor(src, YUVImage, CV_BGR2YUV); //BGR→YUV
cvSplit(YUVImage, Y, U, V, NULL);//分割通道
CvMat* MatY = cvCreateMat(Y->height, Y->width, CV_64FC1);
CvMat* MatU = cvCreateMat(V->height, U->width, CV_64FC1);
CvMat* MatV = cvCreateMat(V->height, V->width, CV_64FC1);
CvMat* DCTY = cvCreateMat(Y->height, Y->width, CV_64FC1);
CvMat* DCTU = cvCreateMat(U->height, U->width, CV_64FC1);
CvMat* DCTV = cvCreateMat(V->height, V->width, CV_64FC1);
cvScale(Y, MatY);
cvScale(U, MatU);
cvScale(V, MatV);
cvDCT(MatY, DCTY, CV_DXT_FORWARD); //余弦变换
cvDCT(MatU, DCTU, CV_DXT_FORWARD); //余弦变换
cvDCT(MatV, DCTV, CV_DXT_FORWARD); //余弦变换
//Y 通道压缩
for (int i = 0; i < Y->height; i++)
{
for (int j = 0; j < Y->width; j++)
{
double element = CV_MAT_ELEM(*DCTY, double, i, j);
if (abs(element) < 10)
CV_MAT_ELEM(*DCTY, double, i, j) = 0;
}
}
// U 通道压缩
for (int i = 0; i < U->height; i++)
{
for (int j = 0; j < U->width; j++)
{
double element = CV_MAT_ELEM(*DCTU, double, i, j);
if (abs(element) < 20)
CV_MAT_ELEM(*DCTU, double, i, j) = 0;
}
}
// V 通道压缩
for (int i = 0; i < V->height; i++)
{
for (int j = 0; j < V->width; j++)
{
double element = CV_MAT_ELEM(*DCTV, double, i, j);
if (abs(element) < 20)
CV_MAT_ELEM(*DCTV, double, i, j) = 0;
}
}
cvDCT(DCTY, MatY, CV_DXT_INVERSE); //余弦反变换
cvDCT(DCTU, MatU, CV_DXT_INVERSE);
cvDCT(DCTV, MatV, CV_DXT_INVERSE);
cvScale(MatY, Y);
cvScale(MatU, U);
cvScale(MatV, V);
cvMerge(Y, U, V, NULL, YUVImage);
cvCvtColor(YUVImage, dst, CV_YUV2BGR); //YUV→BGR
// 计算前后两幅图像的PSNR值
for (int i = 0; i < src->height; i++)
{
for (int j = 0; j < src->width; j++)
{
SrcPixel = cvGet2D(src, i, j);
DstPixel = cvGet2D(dst, i, j);
sumB += (SrcPixel.val[0] - DstPixel.val[0]) * (SrcPixel.val[0] - DstPixel.val[0]);
sumG += (SrcPixel.val[1] - DstPixel.val[1]) * (SrcPixel.val[1] - DstPixel.val[1]);
sumR += (SrcPixel.val[2] - DstPixel.val[2]) * (SrcPixel.val[2] - DstPixel.val[2]);
}
}
mseB = sumB / ((src->width) * (src->height)); //计算均方差
mseG = sumG / ((src->width) * (src->height));
mseR = sumR / ((src->width) * (src->height));
PSNR_B = 10.0 * (log10(255.0 * 255.0 / mseB));
PSNR_G = 10.0 * (log10(255.0 * 255.0 / mseG));
PSNR_R = 10.0 * (log10(255.0 * 255.0 / mseR));
PSNR = (PSNR_B + PSNR_G + PSNR_R) / 3;
printf("PSNR:%d ", PSNR_B);
cvShowImage("YImage", Y);
cvShowImage("UImage", U);
cvShowImage("VImage", V);
cvShowImage("DstImage", dst);
cvSaveImage("F:/dstdemo.jpg", dst);
while (1)
{
if (cvWaitKey(0) == 27) break;
}
cvDestroyWindow("YImage");
cvDestroyWindow("UImage");
cvDestroyWindow("VImage");
cvDestroyWindow("DstImage");
cvReleaseImage(&Y);
cvReleaseImage(&U);
cvReleaseImage(&V);
cvReleaseImage(&src);
cvReleaseImage(&dst);
cvReleaseImage(&YUVImage);
system("pause");
return 0;
}
主要应用
离散余弦变换,尤其是它的第二种类型,经常被信号处理和图像处理使用,用于对信号和图像(包括静止图像和运动图像)进行有损数据压缩。这是由于离散余弦变换具有很强的"能量集中"特性:大多数的自然信号(包括声音和图像)的能量都集中在离散余弦变换后的低频部分,而且当信号具有接近马尔科夫过程(Markov processes)的统计特性时,离散余弦变换的去相关性接近于K-L变换(Karhunen-Loève 变换--它具有最优的去相关性)的性能。
例如,在静止图像编码标准JPEG中,在运动图像编码标准MJPEG和MPEG的各个标准中都使用了离散余弦变换。在这些标准制中都使用了二维的第二种类型离散余弦变换,并将结果进行量化之后进行熵编码。这时对应第二种类型离散余弦变换中的n通常是8,并用该公式对每个8x8块的每行进行变换,然后每列进行变换。得到的是一个8x8的变换系数矩阵。其中(0,0)位置的元素就是直流分量,矩阵中的其他元素根据其位置表示不同频率的交流分类。
一个类似的变换, 改进的离散余弦变换被用在高级音频编码(AAC for Advanced Audio Coding),Vorbis 和 MP3 音频压缩当中。
离散余弦变换也经常被用来使用谱方法来解偏微分方程,这时候离散余弦变换的不同的变量对应着数组两端不同的奇/偶边界条件。