timefft.c

基于DSP 6711 浮点的FFT计算

/**************************
	按时间抽取的FFT变换 
	输入反序,输出正序   
	吴臻志作于2006/12/8
***************************/
#include "stdio.h"
#include "math.h"
struct complex{							//构造复数结构
float real;
float imag;
};
typedef struct complex COMPLEX;

//寻找对应当前下标的反转下标
int FindIndex(int k){
	int j;
	k=k&(0x3ff);
	j=((k&0x001)<<9)+((k&0x002)<<7)+((k&0x004)<<5)+((k&0x008)<<3)+((k&0x010)<<1)+((k&0x020)>>1)+((k&0x040)>>3)+((k&0x080)>>5)+((k&0x100)>>7)+((k&0x200)>>9);
	return j;
}
//反转下标
void ReverseIndex(COMPLEX *x,int N){
int i,j;
COMPLEX tmp;
for (i=0;i<N;i++){
	j=FindIndex(i);
	if (j>i){
		tmp.real=x[i].real;
		tmp.imag=x[i].imag;
		x[i].real =x[j].real;
		x[i].imag =x[j].imag;
		x[j].real=tmp.real;
		x[j].imag=tmp.imag;
		}
	}

}
//按时间抽取法的fft变换，输入反序，输出正序
void fft(COMPLEX *x,int N){				//x为欲变换数组，也当作输出数组，N为点数，支持1024（其它的点数请自行修改反转下标函数即可）
	int j,k,u=0,l=0,wi=0;				//j第二层循环（子块中的每个蝶形的循环计数）
										//k第一层循环（横向fft变换阶数，为log2（N）
										//u 蝶形上标x[upper],l 蝶形下标x[lower]，wi旋转因子下标wn[wi]
	int SubBlockNum,SubBlockStep=1;		//SubBlockNum当前k层子块数量，SubBlockStep当前k层不同子块的相同位置元素间间隔
	COMPLEX wn[512],XkWn;				//wn为旋转因子数组，XkWn为临时存储当前蝶形的旋转因子的临时变量			
	for (k=0;k<N/2;k++){				//初始化wn数组
			wn[k].real=(float)cos(2*3.14159*k/N);
			wn[k].imag=(float)sin(2*3.14159*k/N);
		}
	ReverseIndex(x,N);					//输入反序
	for(k=N;k>1;k=(k>>1)){				//第一个循环，代表log2(k)阶的变换
		SubBlockNum=k>>1;				//子块个数为所做点数的一半
		SubBlockStep=SubBlockStep*2;	//子块间同等地位的元素间隔以2为倍数递增
		wi=0;							//旋转因子初始化
		for(j=0;j<SubBlockStep>>1;j++){	//第二层循环，更新j值，j表示各个子块的第j个蝶形。因为每个子块的同地位蝶形具有相同的wn，所以用第二层循环控制wn
			for(u=j;u<N;u+=SubBlockStep){		//第三层循环，循环于各个子块间的第j个蝶形，计算所有蝶形。直到下标u越界。(u>N)
						l=u+(SubBlockStep>>1);	//下标l计算
						XkWn.real=(x[l].real*wn[wi].real)-(x[l].imag*wn[wi].imag);//蝶形x[u]=x[u]+x[l]*Wn,x[l]=x[u]-x[l]*Wn的复数计算
						XkWn.imag=(x[l].imag*wn[wi].real)+(x[l].real*wn[wi].imag);
						x[l].real=(x[u].real-XkWn.real);
						x[l].imag=(x[u].imag-XkWn.imag);
						x[u].real=(x[u].real+XkWn.real);
						x[u].imag=(x[u].imag+XkWn.imag);
			}
			wi+=SubBlockNum;			//第二层循环更新wi值
		}
	}

}

void test2(){
int k;
int N;	
COMPLEX x[4];
		N=4;
		for (k=0;k<N;k++){
			x[k].real=k;
			x[k].imag=0;
		}
		fft(x,N);
		printf("");
}
void test(){							//测试函数，用于测试sin抽样函数的1024点fft
int k;
float i;
int N;	
COMPLEX x[1024];
		N=1024;
		i=0;
		for (k=0;k<N;k++){
			x[k].real=(float)sin(i);
			x[k].imag=0;
			i+=0.1;
		}

		fft(x,N);						//请在此设置断点，可观察变换前的sin抽样图形。图像参数start address:x,Buffer Size &Display Data Size:1024,IncreaseIndex:2,数据类型float							

}										//请在此设置断点，可观察变换后的频谱，图像参数不变