图像处理与模式识别作业二：快速傅立叶变换FFT与离散余弦变换DCT

二维快速傅立叶变换

算法原理

长度为N的有限序列x(n)的DFT为

X(n)=∑n=0N−1x(n)WnkN，其中WnkN=e−j2ΠNnk

其逆变换IDFT为

x(n)=1N∑n=0N−1X(n)W−nkN

由公式可知，正逆变换的运算量是相同的，计算序列X(n)的时间复杂度为O(N^2)，这样的计算量太庞大，计算速度太慢了，因此我们要利用DFT的一些系数特性来简化计算。

特性如下：

1. 共轭对称性：(WnkN)∗=(W−nkN)

2. 周期性：WnkN=W(n+N)kN=Wn(k+N)N

3. 可约性：WnkN=WnmkN=Wnk/mN/m

得出：

Wn(N−k)N=W(N−n)kN=W−nkNWN/2N=−1W(k+N/2)N=−WkN

利用上述系数特性可以合并DFT计算过程中的一些计算项，从而减少计算量。

根据项的抽选和合并的选择不同，FFT的实现方法也不止一种，下面给出按时间抽选DIT的基-2 FFT算法，也就是最经典的库利-图基算法的思路：

首先解释基-2的含义，即序列长度N=2^r
之所以这样设置是因为算法中使用了二分递归，长度为2的幂次方时运算最方便。
当序列长度为非基-2时，我们会将序列延拓到基-2长度，多出的位置上补0。

算法部分：

先将x(n)序列按n的奇偶性分成0，2，4…和1，3，5…两组子序列，每组的长度为N/2-1

即X(k)=∑N−1n=0x(n)WnkN=∑N/2−1i=0x1(i)(W2N)ik+WkN∑N/2−1j=0x2(j)(W2N)jk

利用系数W的可约性：W2N=W1N/2

将上式化简为X(k)=∑n=0N−1∑i=0N/2−1x1(i)WikN/2+WkN∑j=0N/2−1x2(j)WjkN/2=X1(k)+WkNX2(k)

再利用系数W的周期性：WikN/2=Wi(k+N/2)N/2

再由欧拉公式：eix=cosx+isinx可推导出：

X1(k+N/2)=X1(k)X2(k+N/2)=X2(k)WN/2+kN=WN/2NWkN=−WkN

因此序列X(n)以N/2为界，前半部分为X(k)=X1(k)+WkNX2(k)

后半部分为X(k+N/2)=X1(k)−WkNX2(k)

这意味着我们只需要求出前半部分的X1和X2，就可以算出X(n)后半部分的所有X(k)

根据这两个式子正负交叉的特点，我们可以画出运算流程图如下：



该流程图称为“蝶形计算图”

对长度为8的序列，流程图为：



按照以下算法步骤完成计算：

- 迭代次数r确定

- 对偶节点的计算

- 加权系数W_N^P的计算

- 重新排序

r= log2 N ，r取整数

对于节点x_I (k)，I为迭代次数，k为流程图的行数（也是未排序的数组下标），则xI(k)的对偶节点为xI(k+N/2I) 例：对于长度为8的序列，节点x(0)的对偶节点为x(4)

加权系数W_N^P的计算：先要确定p

(1)把k值写成r位二进制数

(2)把这个二进制数右移r-l位

(3)二进制数按比特位倒转

(4)倒转后的二进制数变为十进制数

例：

求k=2,l=2,N=8的加权系数WP8的计算：

(1)把k值写成r位二进制数：010

(2)把这个二进制数右移r-l=1位：001

(3)二进制数按比特位倒转：100

(4)倒转后的二进制数变为十进制数：4

利用公式计算每一级迭代的节点xi(k+N2i)=xi−1(k)−WpNxi−1(k+N2i)

重新排序

将节点xl(k)的k变为二进制数

xl(k)=xl(kr−1kr−2…k0)

将二进制数按比特位倒转

xl(k)=xl(k0k1…kr−1)

将二进制数变为十进制数，即为重新排序位置

例：

x3(0)→ x3(000) → x3(000) → x3(0)

x3(1)→ x3(001) → x3(100) → x3(4)

x3(2)→ x3(010) → x3(010) → x3(2)

x3(3)→ x3(011) → x3(110) → x3(6)

x3(5)→ x3(101) → x3(101) → x3(5)

x3(6)→ x3(110) → x3(011) → x3(3)

x3(7)→ x3(111) → x3(111) → x3(7)

源码解读

我自己定义了复数类型，并实现了相关的加、减、乘、交换、复制运算

class myComplex {
public :
double real;
double imagin;
};

对于算法中所需要的二进制反转、计算P、计算WPN也编写了专门的函数

/* 二进制反转 */
int reverseNum(int k, int r){
unsigned int ret = 0;
for (int i = 0; i < r; i++){
ret <<= 1;
ret |= (k >> i) & 1;
}
return ret;
}

int getP(int k, int I, int r){
int t = k >> (r - I);
return reverseNum(t, r);
}

void getWn(int n, int N, myComplex * dst){
double x = 2.0 * M_PI * n / N;
dst->imagin = -sin(x);
dst->real = cos(x);
}

按照算法流程图，一维递归实现FFT：

void FFT(myComplex* src, int N, int I, int r, int idx, int MAX_N){
if (N == 1)return;
myComplex *tt1 = new myComplex();
myComplex *tt2 = new myComplex();
myComplex* t1 = new myComplex();
//cout << "idx=" << idx << endl;
for (int i = 0; i < (N >> 1); i ++){
int k = i + idx;
int dst = k + (MAX_N >> I);
int P = getP(k, I, r);
int dstP = getP(dst, I, r);

//伪代码temp[k] = src[k] + getWn(P) * src[dst];
getWn(P, MAX_N, t1);
Mul_Complex(t1, &src[dst], t1);
Add_Complex(&src[k], t1, tt1);
//伪代码temp[dst] = src[k] - getWn(dstP) * src[dst];
getWn(dstP, MAX_N, t1);
Mul_Complex(t1, &src[dst], t1);
//Sub_Complex(&src[k], t1, &temp[dst]); //  不需要用到减法，因为W4 = -W0
Add_Complex(&src[k], t1, tt2);

Copy_Complex(tt1, &src[k]);
Copy_Complex(tt2, &src[dst]);
}
//free(temp);
delete(tt1);
delete(tt2);
delete(t1);
FFT(src, N>>1, I+1, r, idx, MAX_N);
FFT(src, N>>1, I+1, r, idx + (N>>1), MAX_N);
}

经过查阅资料了解到，二维FFT是对原矩阵先每行做一维FFT，再对每一列做一次FFT。由于数组指针的特点，对每列做操作需要先将矩阵转置，然后做每一行的FFT。

结果分析

一维FFT:



二维FFT:



图像FFT:

个人讨论

由于傅立叶变换的特点，可以在转换得到的频谱中添加信息从而达到图片防伪的效果，对于个别有规律的图片噪声也可以通过对频谱的操作完成去噪。

离散余弦变换

算法原理

(摘自《信号与系统》郑君里第二版)

由式可知，可以先求延拓后2N序列的FFT，然后再求得C(k)

同理，IDCT表达式为：



可以先对2N长度的序列C(k)的每一项进行预处理，乘上ejkπ2N，然后对序列进行IFFT，二维的计算方法和二维FFT相同

源码解读

void DCT(double *src, double *dst, int N){
double dTemp;
myComplex *temp = (myComplex*)malloc(sizeof(myComplex)*N*2);
myComplex *temp2 = (myComplex*)malloc(sizeof(myComplex)*N * 2);
memset(temp, 0, sizeof(myComplex)*N *2);
memset(temp2, 0, sizeof(myComplex)*N*2);
for (int i = 0; i < N; i++){
temp[i].real = src[i];
temp[N + i].real = 0;
}

dTemp = 1.0 / sqrt(N);

for (int i = 0; i < N; i++)
dst[0] += temp[i].real;
dst[0] *= dTemp;

realFFT(temp, temp2, N<<1);
dTemp *= sqrt(2);

myComplex *t = new myComplex();
for (int i = 1; i < N; i++){
getWn(i, N<<2, t);
Mul_Complex(t, &temp2[i], t);
dst[i] = dTemp * (t->real);
}

free(temp);
free(temp2);
}

void IDCT(double *src, double *dst, int N){
double N1 = 1 / sqrt(N);
double N2 = sqrt(2) / sqrt(N);
myComplex * t = new myComplex();
myComplex * tC = new myComplex();
myComplex *temp = (myComplex*)malloc(sizeof(myComplex)*N * 2);
myComplex *temp2 = (myComplex*)malloc(sizeof(myComplex)*N * 2);
double *tsrc = (double*)malloc(sizeof(double)*N * 2);
memset(tsrc, 0, sizeof(double)*N * 2);
for (int i = 0; i < N; i++)
tsrc[i] = src[i];
for (int i = 0; i < (N<<1); i++){
getWn(-i, N << 2, t); // 这里要用-i

tC->real = tsrc[i];
tC->imagin = 0;

Mul_Complex(t, tC, &temp[i]);
}
IFFT(temp, temp2, N << 1);

for (int i = 0; i < N; i++){
dst[i] = (N1 - N2) * tsrc[0] + N2 * temp2[i].real * (N<<1);
}

free(temp);
free(temp2);
free(tsrc);
}

结果分析

一维DCT:



二维DCT:



图片DCT:

个人讨论

DCT的优势在于，相对于FFT转换后减少了大量冗余信息，且能量集中于左上角，右下角信息损失对于还原后的图片影响很小，因而多用于有损压缩。一个很有名的应用实例为JPEG图片的压缩算法。

完整代码

/*
FFT DCT
author: Billow_Tau
编译环境：Visual Studio 2013
第三方库：opencv
*/
#define _USE_MATH_DEFINES
#include <iostream>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
#include <math.h>
#include<iomanip>

#define NUM 4

using namespace std;
using namespace cv;

class myComplex {
public :
double real;
double imagin;
};

void Show_Complex(myComplex * src, int size_n){
if (size_n == 1){
if (src->imagin >= 0.0)
printf("%.4lf+%.4lfj  ", src->real, src->imagin);
else
printf("%.4lf%.4lfj  ", src->real, src->imagin);

}
else if (size_n>1){
for (int i = 0; i<size_n; i++)
if (src[i].imagin >= 0.0){
printf("%.4lf+%.4lfj  ", src[i].real, src[i].imagin);
}
else
printf("%.4lf%.4lfj  ", src[i].real, src[i].imagin);
}
}
void MyProc(myComplex *src, int nr, int nc){
const int A = 0;
int half_nr = (nr >> 1) - A;
int half_nc = (nc >> 1) - A;
for (int i = 0; i < (A << 1); i++){
for (int j = 0; j < (A << 1); j++){
src[(half_nr + i)*nc + (half_nc + j)].real = 0;
src[(half_nr + i)*nc + (half_nc + j)].imagin = 0;
}
}
}
/* 复数加法 */
void Add_Complex(myComplex * src1, myComplex *src2, myComplex *dst){
dst->imagin = src1->imagin + src2->imagin;
dst->real = src1->real + src2->real;
//cout << "add:" << endl;
//Show_Complex(dst, 1);
//cout << endl;
}
/* 复数减法 */
void Sub_Complex(myComplex * src1, myComplex *src2, myComplex *dst){
dst->imagin = src1->imagin - src2->imagin;
dst->real = src1->real - src2->real;
//cout << "sub:" << endl;
//Show_Complex(dst, 1);
//cout << endl;
}
/* 复数乘法  */
void Mul_Complex(myComplex * src1, myComplex *src2, myComplex *dst){
double r1 = 0.0, r2 = 0.0;
double i1 = 0.0, i2 = 0.0;
r1 = src1->real;
r2 = src2->real;
i1 = src1->imagin;
i2 = src2->imagin;
dst->imagin = r1*i2 + r2*i1;
dst->real = r1*r2 - i1*i2;
//cout << "mul:" << endl;
//Show_Complex(dst, 1);
//cout << endl;
}
/* 复数赋值 */
void Copy_Complex(myComplex * src, myComplex *dst){
dst->real = src->real;
dst->imagin = src->imagin;
}
void Swap_Complex(myComplex *a, myComplex*b){
myComplex *t = new myComplex();
Copy_Complex(a, t);
Copy_Complex(b, a);
Copy_Complex(t, b);
}
/* 二进制反转 */
int reverseNum(int k, int r){
unsigned int ret = 0;
for (int i = 0; i < r; i++){
ret <<= 1;
ret |= (k >> i) & 1;
}
return ret;
}
int getP(int k, int I, int r){
int t = k >> (r - I);
return reverseNum(t, r);
}
void getWn(int n, int N, myComplex * dst){
double x = 2.0 * M_PI * n / N;
dst->imagin = -sin(x);
dst->real = cos(x);
}
/*取共轭*/
void getConjugate(myComplex* src){
src->imagin = -(src->imagin);
}
void FFT(myComplex* src, int N, int I, int r, int idx, int MAX_N){
if (N == 1)return;
myComplex *tt1 = new myComplex();
myComplex *tt2 = new myComplex();
myComplex* t1 = new myComplex();
//cout << "idx=" << idx << endl;
for (int i = 0; i < (N >> 1); i ++){
int k = i + idx;
int dst = k + (MAX_N >> I);
int P = getP(k, I, r);
int dstP = getP(dst, I, r);

//伪代码temp[k] = src[k] + getWn(P) * src[dst];
getWn(P, MAX_N, t1);
Mul_Complex(t1, &src[dst], t1);
Add_Complex(&src[k], t1, tt1);
//伪代码temp[dst] = src[k] - getWn(dstP) * src[dst];
getWn(dstP, MAX_N, t1);
Mul_Complex(t1, &src[dst], t1);
//Sub_Complex(&src[k], t1, &temp[dst]); //  不需要用到减法，因为W4 = -W0
Add_Complex(&src[k], t1, tt2);

Copy_Complex(tt1, &src[k]);
Copy_Complex(tt2, &src[dst]);
}
//free(temp);
delete(tt1);
delete(tt2);
delete(t1);
FFT(src, N>>1, I+1, r, idx, MAX_N);
FFT(src, N>>1, I+1, r, idx + (N>>1), MAX_N);
}
void realFFT(myComplex *src, myComplex *dst, int N){
int r = log10(N) / log10(2); // log2(N) = lg(N) / lg(2)

//Show_Complex(temp, N);

FFT(src,N,1,r,0, N);
for (int i = 0; i < N; i++){
//cout << reverseNum(i, r) << ' ';
Copy_Complex(&src[i], &dst[reverseNum(i, r)]);
}
}
void IFFT(myComplex* src, myComplex *dst, int N){
double t = 1.0 / N;
for (int i = 0; i < N; i++){
getConjugate(&src[i]);
}
realFFT(src, dst, N);
for (int i = 0; i < N; i++){
getConjugate(&dst[i]);
dst[i].real *= t;
dst[i].imagin *= t;
}
}
void FFT2D(myComplex *src, myComplex *dst, int nr, int nc){
myComplex* temp = (myComplex*)malloc(sizeof(myComplex)* nr*nc);
for (int i = 0; i < nr*nc; i++){
temp[i].real = src[i].real;
temp[i].imagin = src[i].imagin;
}
// 先对每行进行一维FFT
for (int i = 0; i < nr; i++)
realFFT(&src[i*nc], &temp[i*nc], nc);

//cout << "------------------------mid:--------------------------" << endl;
//for (int i = 0; i < NUM; i++, cout << endl)
//  Show_Complex(&temp[i*NUM], NUM);
//cout << "------------------------mid:--------------------------" << endl;

// 然后转置
for (int i = 0; i < nr; i++){
for (int j = i + 1; j < nc; j++){
Swap_Complex(&temp[j*nr + i], &temp[i*nc + j]);
}
}
// 再次对每行进行一维FFT，等同于对转置前的矩阵每列进行FFT
for (int i = 0; i < nr; i++)
realFFT(&temp[i*nc], &dst[i*nc], nc);

free(temp);

// 转置回来
for (int i = 0; i < nr; i++){
for (int j = i + 1; j < nc; j++){
Swap_Complex(&dst[j*nr + i], &dst[i*nc + j]);
}
}
}

void IFFT2D(myComplex *src, myComplex *dst, int N){
myComplex* temp = (myComplex*)malloc(sizeof(myComplex)* N*N);
// 先转置
for (int i = 0; i < N; i++){
for (int j = i + 1; j < N; j++){
Swap_Complex(&src[j*N + i], &src[i*N + j]);
}
}
// 再次对每行进行一维IFFT，等同于对转置前的矩阵每列进行IFFT
for (int i = 0; i < N; i++)
IFFT(&src[i*N], &temp[i*N], N);
// 取矩阵的共轭
//for (int i = 0; i < N*N; i++){
//  getConjugate(&temp[i]);
//}
// 再次转置
for (int i = 0; i < N; i++){
for (int j = i + 1; j < N; j++){
Swap_Complex(&temp[j*N + i], &temp[i*N + j]);
}
}
// 再次对每行进行一维IFFT，等同于对原矩阵每行进行IFFT
for (int i = 0; i < N; i++){
IFFT(&temp[i*N], &dst[i*N], N);
}
}
void FFT_Shift(myComplex *src, int nr, int nc){
for (int i = 0; i < nr; i++){
for (int j = 0; j < nc; j++){
if ((i + j) % 2){
src[i*nc + j].real = -src[i*nc +j].real;
}
}
}
}
void ImageFFT(Mat &src, myComplex *dst){
int nr = src.rows;
int nc = src.cols * src.channels();
myComplex *image_data = (myComplex*)malloc(sizeof(myComplex)*nr*nc);
for (int i = 0; i < nr; i++){
uchar* data = src.ptr<uchar>(i);
for (int j = 0; j < nc; j++){
image_data[i*nc + j].real = data[j];
image_data[i*nc + j].imagin = 0.0;
}
}
FFT_Shift(image_data, nr, nc);
FFT2D(image_data, dst, nr, nc);
free(image_data);
}
void ImageIFFT(myComplex *src, Mat &dst, int nr, int nc){
myComplex *image_data = (myComplex*)malloc(sizeof(myComplex)*nr*nc);
IFFT2D(src, image_data, nr);
FFT_Shift(image_data, nr, nc);
for (int i = 0; i < nr; i++){
uchar *data = dst.ptr<uchar>(i);
for (int j = 0; j < nc; j++)
data[j] = image_data[i*nc + j].real;
}
free(image_data);
}
/*幅度拉伸*/
void NormalSize(double* src, double *dst, int nr, int nc){
double maxm = 0.0, minm = 256;
for (int i = 0; i < nr*nc; i++){
maxm = src[i] > maxm ? src[i] : maxm;
minm = src[i] < minm ? src[i] : minm;
}
double step = 255.0 / (maxm - minm);
for (int i = 0; i < nr*nc; i++){
dst[i] = (src[i] - minm) * step;
dst[i] *= 45.9*log((double)(1 + dst[i]));
//dst[i] *= log((double)(1 + dst[i]));
if (dst[i] > 255) dst[i] = 255;
else if (dst[i] < 0)    dst[i] = 0;
}
}
/*频谱获得*/
void getAmplitudespectrum(myComplex *src, int nr, int nc, Mat &dst){
//myComplex *despe = (myComplex*)malloc(sizeof(myComplex)*nc * nr);
double *despe = (double*)malloc(sizeof(double)*nc * nr);
if (despe == NULL)exit(0);

double real = 0.0, imagin = 0.0;
for (int i = 0; i < nr; i++){
for (int j = 0; j < nc; j++){
real = src[i*nc + j].real;
imagin = src[i*nc + j].imagin;
despe[i*nc + j] = sqrt(real * real + imagin * imagin );
}
}

NormalSize(despe, despe, nr, nc);
for (int i = 0; i < nr; i++){
uchar *data = dst.ptr<uchar>(i);
for (int j = 0; j < nc; j++){
data[j] = despe[i*nc + j];
}
}
free(despe);
}
void Image(Mat &img){

// image
Mat result(img.cols, img.rows, CV_8U);
Mat result2(img.cols, img.rows, CV_8U);

myComplex *dst = (myComplex*)malloc(sizeof(myComplex)* img.rows * img.cols);
ImageFFT(img, dst);
MyProc(dst, img.rows, img.cols);
getAmplitudespectrum(dst, img.rows, img.cols, result);

ImageIFFT(dst, result2, img.rows, img.cols);

namedWindow("原图", CV_WINDOW_AUTOSIZE);
imshow("原图", img);
namedWindow("傅里叶谱", CV_WINDOW_AUTOSIZE);
imshow("傅里叶谱", result);
namedWindow("逆傅里叶图", CV_WINDOW_AUTOSIZE);
imshow("逆傅里叶图", result2);
waitKey(0);
}
void fft2d(double *src){
myComplex *dst = (myComplex*)malloc(sizeof(myComplex)* NUM*NUM);
myComplex * temp = (myComplex*)malloc(sizeof(myComplex)* NUM*NUM);

for (int i = 0; i < NUM*NUM; i++){
temp[i].real = src[i];
temp[i].imagin = 0.0;
}
cout << "data:" << endl;
for (int i = 0; i < NUM; i++, cout << endl)
Show_Complex(&temp[i*NUM], NUM);
cout << "\nFFT2D:" << endl;

FFT2D(temp, dst, NUM, NUM);

for (int i = 0; i < NUM; i++, cout << endl)
Show_Complex(&dst[i*NUM], NUM);
cout << "\nIFFT2D:" << endl;
IFFT2D(dst, temp, NUM);
for (int i = 0; i < NUM; i++, cout << endl)
Show_Complex(&temp[i*NUM], NUM);

}
void fft1d(double *src, int N){
myComplex *dst = (myComplex*)malloc(sizeof(myComplex)* N);

myComplex* temp = (myComplex*)malloc(sizeof(myComplex)* N);
for (int i = 0; i < N; i++){
temp[i].real = src[i];
temp[i].imagin = 0.0;
}
realFFT(temp, dst, N);
Show_Complex(dst, N);

cout << "\n------------------------" << endl;

myComplex* temp2 = (myComplex*)malloc(sizeof(myComplex)* NUM);
for (int i = 0; i < NUM; i++){
temp2[i].real = src[i];
temp2[i].imagin = 0.0;
}
IFFT(dst, temp2, NUM);
Show_Complex(temp2, NUM);
free(temp);
free(temp2);
free(dst);
}

void DCT(double *src, double *dst, int N){
double dTemp;
myComplex *temp = (myComplex*)malloc(sizeof(myComplex)*N*2);
myComplex *temp2 = (myComplex*)malloc(sizeof(myComplex)*N * 2);
memset(temp, 0, sizeof(myComplex)*N *2);
memset(temp2, 0, sizeof(myComplex)*N*2);
for (int i = 0; i < N; i++){
temp[i].real = src[i];
temp[N + i].real = 0;
}

dTemp = 1.0 / sqrt(N);

for (int i = 0; i < N; i++)
dst[0] += temp[i].real;
dst[0] *= dTemp;

realFFT(temp, temp2, N<<1);
dTemp *= sqrt(2);

myComplex *t = new myComplex();
for (int i = 1; i < N; i++){
getWn(i, N<<2, t);
Mul_Complex(t, &temp2[i], t);
dst[i] = dTemp * (t->real);
}

free(temp);
free(temp2);
}
void IDCT(double *src, double *dst, int N){
double N1 = 1 / sqrt(N);
double N2 = sqrt(2) / sqrt(N);
myComplex * t = new myComplex();
myComplex * tC = new myComplex();
myComplex *temp = (myComplex*)malloc(sizeof(myComplex)*N * 2);
myComplex *temp2 = (myComplex*)malloc(sizeof(myComplex)*N * 2);
double *tsrc = (double*)malloc(sizeof(double)*N * 2);
memset(tsrc, 0, sizeof(double)*N * 2);
for (int i = 0; i < N; i++)
tsrc[i] = src[i];
for (int i = 0; i < (N<<1); i++){
getWn(-i, N << 2, t); // 这里要用-i

tC->real = tsrc[i];
tC->imagin = 0;

Mul_Complex(t, tC, &temp[i]);
}
IFFT(temp, temp2, N << 1);

for (int i = 0; i < N; i++){
dst[i] = (N1 - N2) * tsrc[0] + N2 * temp2[i].real * (N<<1);
}

free(temp);
free(temp2);
free(tsrc);
}
void dct1D(double *src, int N){
double dst[8], ans[8];
cout << "data:" << endl;
for (int i = 0; i < N; i++)
cout << src[i] << ' ';
DCT(src, dst, N);
cout << "\n------------------" << endl;
cout << "DCT:" << endl;
for (int i = 0; i < N; i++)
cout << setprecision(3) << dst[i] << " ";
cout << "\n------------------" << endl;

cout << "IDCT:" << endl;
IDCT(dst, ans, N);
for (int i = 0; i < N; i++)
cout << setprecision(3) << ans[i] << " ";
}
void DCT2D(double *src , double *dst, int nr, int nc){
double *temp = (double *)malloc(sizeof(double)*nr*nc);
for (int i = 0; i < nr; i++){
DCT(&src[i*nc], &temp[i*nc], nc);
}
// 然后转置
for (int i = 0; i < nr; i++){
for (int j = i + 1; j < nc; j++){
double t = temp[j*nr + i];
temp[j*nr + i] = temp[i *nc + j];
temp[i*nc + j] = t;
//Swap_Complex(&temp[j*nr + i], &temp[i*nc + j]);
}
}
for (int i = 0; i < nc; i++){
DCT(&temp[i*nr], &dst[i*nr], nr);
}
for (int i = 0; i < nc; i++){
for (int j = i + 1; j < nr; j++){
double t = dst[i*nr + j];
dst[i*nr + j] = dst[j *nr + i];
dst[j*nr + i] = t;
}
}
free(temp);
}
void IDCT2D(double *src, double*dst, int nr, int nc){
double *temp = (double *)malloc(sizeof(double)*nr*nc);
for (int i = 0; i < nc; i++){
for (int j = i + 1; j < nr; j++){
double t = src[i*nr + j];
src[i*nr + j] = src[j *nr + i];
src[j*nr + i] = t;
}
}
for (int i = 0; i < nr; i++){
IDCT(&src[i*nc], &temp[i*nc], nc);
}
for (int i = 0; i < nr; i++){
for (int j = i + 1; j < nc; j++){
double t = temp[j*nr + i];
temp[j*nr + i] = temp[i *nc + j];
temp[i*nc + j] = t;
//Swap_Complex(&temp[j*nr + i], &temp[i*nc + j]);
}
}
for (int i = 0; i < nc; i++)
IDCT(&temp[i*nr], &dst[i*nr], nr);
}
void dct2D(double *src, int N){
double dst[16], ans[16];
memset(dst, 0, sizeof(double)* 16);
memset(ans, 0, sizeof(double)* 16);
cout << "data:" << endl;
for (int i = 0; i < NUM; i++, cout << endl){
for (int j = 0; j < NUM; j++)
cout << src[i*NUM + j] << ' ';
}
cout << "\n----------------" << endl;
cout << "DCT2D:" << endl;
DCT2D(src, dst, NUM, NUM);
for (int i = 0; i < NUM; i++, cout << endl){
for (int j = 0; j < NUM; j++)
cout << dst[i*NUM + j] << ' ';
}
cout << "\n----------------" << endl;
cout << "IDCT2D:" << endl;
IDCT2D(dst, ans, NUM, NUM);
for (int i = 0; i < NUM; i++, cout << endl){
for (int j = 0; j < NUM; j++)
cout << ans[i*NUM + j] << ' ';
}
}
void ImageDCT(Mat &src, double *dst){
int nr = src.rows;
int nc = src.cols;
double *temp = (double*)malloc(sizeof(double)*nr*nc);
for (int i = 0; i < nr; i++){
uchar *data = src.ptr<uchar>(i);
for (int j = 0; j < nc; j++){
temp[i*nc + j] = data[j];
}
}
DCT2D(temp, dst, nr, nc);
}
void ShowAmplitudespectrum(double *src, Mat &dst, int nr, int nc){
uchar maxn = 0, minn = 255;
double* temp = (double*)malloc(sizeof(double)*nr*nc);
for (int i = 0; i < nr*nc; i++){
maxn = (src[i]> maxn) ? src[i] : maxn;
minn = (src[i] < minn) ? src[i] : minn;
}
for (int i = 0; i < nr*nc; i++){
temp[i] = 255.0 / (maxn - minn) * src[i];
}
for (int i = 0; i < nr; i++){
uchar *data = dst.ptr<uchar>(i);
for (int j = 0; j < nc; j++){
data[j] = temp[i*nc + j];
}
}
free(temp);
}
void ImageIDCT(double *src, Mat &dst){
int nr = dst.rows;
int nc = dst.cols;
double* temp = (double*)malloc(sizeof(double)*nr*nc);
IDCT2D(src, temp, nr, nc);
for (int i = 0; i < nr; i++){
uchar *data = dst.ptr<uchar>(i);
for (int j = 0; j < nc; j++){
data[j] = temp[i*nc + j];
}
}
free(temp);
}
void Image_dct(Mat &img){
double *dst = (double*)malloc(sizeof(double)*img.rows*img.cols);
memset(dst, 0, sizeof(double)*img.cols*img.rows);
ImageDCT(img, dst);
for (int i = 0; i < 512; i++)
cout << dst[i] << ' ';

Mat result(img.cols, img.rows, CV_8U);
Mat result2(img.cols, img.rows, CV_8U);

ShowAmplitudespectrum(dst, result, img.rows, img.cols);

ImageIDCT(dst, result2);

imshow("离散余弦变换", result);
namedWindow("离散余弦变换", CV_WINDOW_AUTOSIZE);
imshow("还原", result2);
namedWindow("还原", CV_WINDOW_AUTOSIZE);

waitKey(0);
}
int main(){

Mat img = imread("D:\\opencv\\project\\lena.jpg", 1);
cvtColor(img, img, CV_BGR2GRAY);

//double src[8] = { 0.80902013, 0.42188063, 0.84532871, 0.24776496, 0.82053718, 0.09025071,
//  0.55021263, 0.55021845 };
//double src[16] = { 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8 };

//1D
//fft1d(src,NUM);

// 2D
//double src[16] = { 47, 80, 44, 34, 23, 21, 3, 57, 5, 96, 54, 41, 93, 24, 36, 20 };
//double src[16] = {1,2,3,0,4,5,6,0,7,8,9,0,0,0,0,0};

//fft2d(src);

//image
//Image(img);

// dct1D
//dct1D(src, NUM);

//dct2D
//dct2D(src, NUM);

//image dct
//Image_dct(img);

system("pause");
return 0;
}