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

来源：互联网发布：石家庄软件培训学校编辑：程序博客网时间：2024/06/05 01:53

二维快速傅立叶变换

算法原理

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

X (n) = \sum n = 0 N - 1 x (n) W n k N ， 其 中 W n k N = e - j 2 Π N n k

其逆变换IDFT为

x (n) = 1 N \sum n = 0 N - 1 X (n) W - n k N

由公式可知，正逆变换的运算量是相同的，计算序列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
得出：

W n (N - k) N = W (N - n) k N = W - n k N

W N / 2 N = - 1

W (k + N / 2) N = - W k N

利用上述系数特性可以合并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) = \sum n = 0 N - 1 \sum i = 0 N / 2 - 1 x 1 (i) W i k N / 2 + W k N \sum j = 0 N / 2 - 1 x 2 (j) W j k N / 2

= X 1 (k) + W k N X 2 (k)

再利用系数W的周期性：

WikN/2=Wi(k+N/2)N/2
再由欧拉公式：

eix=cosx+isinx可推导出：

X 1 (k + N / 2) = X 1 (k) X 2 (k + N / 2) = X 2 (k) W N / 2 + k N = W N / 2 N W k N = - W k N

因此序列X(n)以N/2为界，前半部分为

X(k)=X1(k)+WkNX2(k)
后半部分为

X(k+N/2)=X1(k)−WkNX2(k)

这意味着我们只需要求出前半部分的X1和X2，就可以算出X(n)后半部分的所有X(k)
根据这两个式子正负交叉的特点，我们可以画出运算流程图如下：

Markdown

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

对长度为8的序列，流程图为：
Markdown
按照以下算法步骤完成计算：
- 迭代次数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
利用公式计算每一级迭代的节点
$x i (k + N 2 i) = x i - 1 (k) - W p N x i - 1 (k + N 2 i)$
重新排序
将节点xl(k)的k变为二进制数
$x l (k) = x l (k r - 1 k r - 2 \dots k 0)$
将二进制数按比特位倒转
$x l (k) = x l (k 0 k 1 \dots k r - 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:
Markdown
二维FFT:

图像FFT:

个人讨论

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

离散余弦变换

算法原理

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

由式可知，可以先求延拓后2N序列的FFT，然后再求得C(k)
同理，IDCT表达式为：
Markdown
可以先对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:
Markdown

二维DCT:
Markdown

图片DCT:
Markdown

个人讨论

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

完整代码

/*FFT DCTauthor: 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 4using 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;}

阅读全文

0 0