我是在Matlab上找到灵感

先去网站

http://www.fftw.org/

下载

然后调用实现fftw以及ifftw

然后在Matlab里面找到dct 以及idct的文件

把那几行Matlab改成c语言就行了

DCT变换代码：

// AnyLengFDCT.cpp : 定义控制台应用程序的入口点。

//

#include "stdafx.h"

#include "fftw3.h"//引用库文件

#include

#define PI 3.1415926

static double *as = NULL;

static double *ax = NULL;

static fftw_complex *in, *out;

static int n = 4;

void InitData();

void dft(int dn);

void fft_dct(double *a, int dn);

int _tmain(int argc, _TCHAR* argv[])

{

double a[4] = {1, 2, 3, 4};

InitData();

dft(n);

fft_dct(a, n);

return 0;

}

void InitData(){

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*n); //FFT变换前数据

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*n);//FFT变换后数据

as = (double*)malloc(n*sizeof(double));

ax = (double*)malloc(n*sizeof(double));

}

//*******************************************************************

//方法名称:dft

//说明：对所变幻的数据的长度进行DFT，因为每次变换的数据长度一致，所以单独拿出来，防止重复计算

//日期：2012.4.19

//参数说明:

//as: DFT变换后数据的实数

//ax:DFT变换后数据的虚部

//dn:DFT变换数据的个数

//*******************************************************************

void dft(int dn){

int i;

double x = 0.0;

double y;

//Compute weights to multiply DFT coefficients

for(i = 0; i < dn; i++){

y = -i*PI/(2*dn);

as[i] = exp(x)*cos(y)/sqrt(2.0*dn);

ax[i] = exp(x)*sin(y)/sqrt(2.0*dn);

}

as[0] = as[0]/sqrt(2.0);

for(i = 0; i < dn; i++){

printf("%f, %f\n", as[i], ax[i]);

}

}

//*******************************************************************

//方法名称:fft_dct

//说明：任意长度的快速DCT变换

//日期：2012.4.19

//参数说明:

//a: 数据

//dn:变换数据的大小

//*******************************************************************

void fft_dct(double *a, int dn){

int l;

int j;

int k;

fftw_plan p;

if(dn%2 == 1){

//奇数的DCT变换

for(l = 0; l < dn; l++){

in[l][0] = a[l];

in[dn+l][0] = a[dn-1-l];

}

p=fftw_plan_dft_1d(2*dn,in,out, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(p); /* repeat as needed*/

}else{

//偶数的DCT变换

for(l = 0; l < dn; l++){

as[l] = as[l]*2;

ax[l] = ax[l]*2;

//printf("%f, %f\n", as[i], ax[i]);

}

for(l = 0, j = 0,k = dn-1; l < dn; l++){

if(l%2){

in[k][0] = a[l];

--k;

}else{

in[j][0] = a[l];

++j;

}

}

p=fftw_plan_dft_1d(dn,in,out, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(p); /* repeat as needed*/

}

for(l = 0; l < dn; l++){

a[l] = as[l]*out[l][0] - ax[l]*out[l][1];

//printf("%f\n", a[l]);

}

fftw_destroy_plan(p);

fftw_free(in);

fftw_free(out);

free(as);

free(ax);

}

反DCT变换：

// AnyLengFIDCT.cpp : 定义控制台应用程序的入口点。

//

#include "stdafx.h"

#include "fftw3.h"//引用库文件

#include

#define PI 3.1415926

//static double *ias = NULL;

//static double *iax = NULL;

//static fftw_complex *in, *out;

//static int n = 6;

//void InitData();

//void idft(int dn);

//void ifft_idct(double *a, int dn);

void idft(int dn, double *ias, double *iax);

void ifft_idct(int dn, fftw_plan p, double *a, double *data, double *ias, double *iax, fftw_complex *in, fftw_complex *out);

int _tmain(int argc, _TCHAR* argv[])

{

int n = 6;

//double a[5] = {6.708204, -3.1495, 0, -0.283990, 0};

double a[6] = {8.573214, -4.162562, 0, -0.408248, 0, -0.080079};

double *ias = (double*)malloc(n*sizeof(double));

double *iax = (double*)malloc(n*sizeof(double));

double *data = (double*)malloc(n*sizeof(double));

fftw_plan p = NULL;

fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*n); //IFFT变换前数据

fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2*n);//IFFT变换后数据

idft(n, ias, iax);

ifft_idct(n, p, a, data, ias, iax, in, out);

for(int i = 0; i < n; i++)

printf("%d %f\n", i+1, data[i]);

return 0;

}

//*******************************************************************

//方法名称:idft

//说明：对所变幻的数据的长度进行IDFT，因为每次变换的数据长度一致，所以单独拿出来，防止重复计算

//日期：2012.4.19

//参数说明:

//as: DFT变换后数据的实数

//ax:DFT变换后数据的虚部

//dn:DFT变换数据的个数

//*******************************************************************

void idft(int dn, double *ias, double *iax){

int i;

double x = 0.0;

double y;

//Compute weights to multiply IDFT coefficients

for(i = 0; i < dn; i++){

y = i*PI/(2*dn);

ias[i] = exp(x)*cos(y)*sqrt(2.0*dn);

iax[i] = exp(x)*sin(y)*sqrt(2.0*dn);

}

// for(i =0; i < dn; i++){

// printf("%f, %f\n", ias[i], iax[i]);

// }

}

//*******************************************************************

//方法名称:fft_dct

//说明：任意长度的快速IDCT变换

//日期：2012.4.19

//参数说明:

//a: 数据

//dn:变换数据的大小

//*******************************************************************

void ifft_idct(int dn, fftw_plan p, double *a, double *data, double *ias, double *iax, fftw_complex *in, fftw_complex *out){

int l;

int j;

int k;

double ias_0;

if(dn%2 == 1){

//奇数的DCT变换

ias_0 = ias[0] * sqrt(2.0);

in[0][0] = ias_0*a[0];

in[0][1] = iax[0]*a[0];

for(l = 1; l < dn; l++){

in[l][0] = ias[l]*a[l];

in[l][1] = iax[l]*a[l];

//printf("1==%f, 2==%f\n", in[l][0], in[l][1]);

in[dn+l][0] = iax[l]*a[dn-l];

in[dn+l][1] = -ias[l]*a[dn-l];

//printf("a==%f, b==%f\n", in[dn+l][0], in[dn+l][1]);

}

p=fftw_plan_dft_1d(2*dn,in,out, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(p); /* repeat as needed*/

for(l = 0; l < dn; l++){

data[l] = out[l][0] / dn/2;

//out[l][1] = out[l][1] / dn;

//data[l] = out[l][0];

//printf("data == %f\n", data[l]);

}

}else{

//偶数的DCT变换

ias_0 = ias[0] / sqrt(2.0);

in[0][0] = ias_0*a[0];

in[0][1] = iax[0]*a[0];

for(l = 1; l < dn; l++){

in[l][0] = ias[l]*a[l];

in[l][1] = iax[l]*a[l];

//printf("%f, %f\n", in[l][0], in[l][1]);

}

p=fftw_plan_dft_1d(dn,in,out, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(p); /* repeat as needed*/

for(l = 0; l < dn; l++){

out[l][0] = out[l][0] / dn;

out[l][1] = out[l][1] / dn;

//printf("out == %f\n", out[l][0]);

}

for(l = 0, j = 0, k = dn-1; l < dn; l++){

if(l%2){

data[l] = out[k][0];

k--;

}else{

data[l] = out[j][0];

++j;

}

//printf("data == %f\n", data[l]);

}

}

}

注：记得下载FFTW以及IFFTW