任意长度DCT变换C语言实现方法

我是在Matlab上找到灵感
先去网站 http://www.fftw.org/  下载
然后调用实现fftw以及ifftw
然后在Matlab里面找到dct 以及idct的文件

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

DCT变换代码：

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

#include "stdafx.h"
#include "fftw3.h"//引用库文件
#include <string>
#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 <string>
#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