ifftc.cpp

在TMS6713上实现IFFT的C语言程序

#include<stdio.h>
#include<math.h>
//#include<IFFTC.h>
/**************************************/
/*              定义                  */
/**************************************/
#define N 8                         //fft 点数
#define pi 3.141592653589793

typedef struct{
    float real;
    float image;
}Complex;

Complex ComMul(Complex p1,Complex p2);
Complex ComAdd(Complex p1,Complex p2);
Complex ComSub(Complex p1,Complex p2);

void FFTcompute(Complex data_in[], Complex data_out[], int len);
void ResetSrc(Complex data_src[], int len, int size);
int BitReverse(int src, int size);
int	ComputeGradeNum(int len);
void CalcW(Complex *W,int len);
void Real2Complex(float data_real[], Complex data_plex[], int len);


//******************************思路说明*********************************/
//GradeNum--蝶形运算级数; 
//CellNum--每一级子fft个数；
//Subfftlen--子fft长度；
//在fft蝶形运算中，其运算级数为log2(N)；
//自左向右每一级里又有CellNum个子fft运算(CellNum=2^(k-h),基2自左向右减少)；
//每一级的子fft运算长度为Subfftlen(Subfftlen=2^(h)),基2自左向右增加,根据fft两个source num经不同的运算得到两个目标值)；
//CellNum*Subfftlen=N
//----------------------------------------
//1)利用for循环自左向右计算每一级每个fft
//2)计算前需对输入序列重排
//3)蝶形运算为原位运算，每次算出来的结果直接放入原寄存器
//4)每一级的比例因子为W(j),       =W(j*CellNum);先把所用到的W值(N/2个)都计算出来，用到的时候直接取 
//                    W(Subfftlen) W(N)
/***********************************************************************/
main()
{

	float data_in[N];
	Complex data_out[N];
	Complex data_in_plex[N];
	int i;
	FILE* fp;


	 if((fp=fopen("C:\MATLAB6p5\work\input","r"))==NULL)
    {
        printf("cannot open file");
  //      exit(0);
    }

    for(i=0;i<N;i++)
        scanf("%f",&data_in[i]);

	Real2Complex(data_in, data_in_plex, N);

	FFTcompute(data_out, data_in_plex, N);
	
//-----------output data---------------------------//
	 for(i=0;i<N;i++)
        printf("%f, %f, \n",data_out[i].real, data_out[i].image);
}

/*****************************************************************/
/*                           子函数                              */
/*****************************************************************/
/*--------FFT core 算法--------*/
/*In:data_in[]--数据列
             N--数据列(FFT)长度              
/*Out:data_out[]--数据列after fft
/*-----------------------------*/
void FFTcompute(Complex data_out[], Complex data_in[], int len)
{
	
	Complex W[N/2];                   //比例因子W
	Complex Wmul;
	Complex temp;

	int tmp;
	int	index1, index2;
	int i, j, h;
	int CellNum, Subfftlen, GradeNum;

//-----------Preparation-------------------------//
	CalcW(W, len);								    //计算所用到的W值(N/2个)
	GradeNum = ComputeGradeNum(len);              //Compute 级数
//-----------输入序列重排(index比特反转)--------//
	ResetSrc(data_in, len, GradeNum);
	
//--------- fft computation----------------------//
	for(h=1; h<=GradeNum; h++)
	{
		CellNum = pow(2, (GradeNum-h));                                  
		Subfftlen = pow(2, h);

		for(i=0; i<CellNum; i++)
			for(j=0; j<Subfftlen/2; j++)									//做一次运算占用两个源数 
			{

				index1 = j+i*Subfftlen;										// fft数1index；参见蝶形运算图理解index
				index2 = j+(Subfftlen/2)+i*Subfftlen;						// fft数2index
				
				Wmul = ComMul(W[CellNum*j], data_in[index2]);                // W()
//				data_in[index1] = ComAdd(data_in[index1], Wmul);
				temp = ComAdd(data_in[index1], Wmul);
				data_in[index2] = ComSub(data_in[index1], Wmul);
				data_in[index1] = temp;
			}
	}


	for(tmp=0; tmp<N; tmp++)
	{
		data_out[tmp] = data_in[tmp];
	}


}
/*--------FFT级数--------*/
/*In:len--FFT长度              
/*Out:M--级数
/*-----------------------*/
int	ComputeGradeNum(int len)
{
	int M=0;
	M = int((log10(len))/(log10(2)));
	return M;
}
/*--------W比例因子计算--------*/
/*In:len--FFT长度              
/*Out:W--比例因子矩阵
/*-----------------------------*/
void CalcW(Complex *W,int len)
{
	int size = len/2;
	for(int tmp=0; tmp<size; tmp++)
	{
		W[tmp].real = cos(pi*tmp/size);
		W[tmp].image = -1*sin(pi*tmp/size);
	}
	
}
/*--------W比例因子计算--------*/
/*In:data_src--数据列;
          len--数据列长度
		 size--index比特宽度	
/*Out:data_src--重排后数据列
/*-----------------------------*/
void ResetSrc(Complex data_src[], int len, int size)
{
	for(int tmp=0; tmp<N; tmp++)
	{
		Complex data_temp;
		int index_reverse = 0;

		index_reverse = BitReverse(tmp, size);
		if(index_reverse>tmp)
		{
			data_temp = data_src[index_reverse];
			data_src[index_reverse] = data_src[tmp];
			data_src[tmp] = data_temp;
		}

	}

}
/*-------------(Num比特反转)-----------------*/
/*In:src--source num
/*   size--数的bit宽度               
/*Out:des--destination num
/*-------------------------------------------*/
int BitReverse(int src, int size)
{
	int temp = src;
	int des = 0;
	int tmp = size-1;
	for(tmp;tmp>=0;tmp--)                  //利用比特移位
	{
		des = ((temp&0x1)<<tmp)|des;
		temp = temp>>1;
	}
	return des;
}
/*--------复数计算(乘，加，减)---------------*/
/*In:Complex p1,p2--source num
/*Out:Res--Computation Result(Complex)
/*-------------------------------------------*/
Complex ComMul(Complex p1,Complex p2)
{
	Complex	Res;
	Res.real = p1.real*p2.real-p1.image*p2.image;
	Res.image = p1.real*p2.image+p1.image*p2.real;
	return Res;
}

Complex ComAdd(Complex p1,Complex p2)
{
    Complex Res;
    Res.real = p1.real + p2.real;
    Res.image = p1.image + p2.image;
    
    return Res;
}

Complex ComSub(Complex p1,Complex p2)
{
    Complex Res;
    Res.real = p1.real - p2.real;
    Res.image = p1.image - p2.image;
    
    return Res;
}
/*--------------实数转复数---------------*/
/*In:data_real[]--source real num
             len--length of data stream
/*Out:data_plex--destination complex num
/*-------------------------------------------*/
void Real2Complex(float data_real[], Complex data_plex[], int len)
{
	for(int tmp=0; tmp<len; tmp++)
	{
		data_plex[tmp].image = 0;
		data_plex[tmp].real = data_real[tmp];
	
	}
	
}