我是在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