timefft.h

32位嵌入式系统实现

#include "stdio.h"
#include "math.h"
#ifndef FFT_Def
#define FFT_Def
#define N 512						//FFT点数

#if (N==2048)
#define FindIndex FindIndex11bit
#elif (N==1024)
#define FindIndex FindIndex10bit
#elif  (N==512)
#define FindIndex FindIndex9bit
#elif  (N==256)
#define FindIndex FindIndex8bit
#elif  (N==128)
#define FindIndex FindIndex7bit
#elif (N==64)
#define FindIndex FindIndex6bit
#endif
struct Complex{							//构造复数结构
float real;
float imag;
};
typedef struct Complex COMPLEX;
#define uchar unsigned char
//寻找对应当前下标的反转下标

int FindIndex11bit(int k){
	int j;
	k=k&(0x7ff);
	j=((k&0x001)<<10)+((k&0x002)<<8)+((k&0x004)<<6)+((k&0x008)<<4)+((k&0x010)<<2)+((k&0x020))+((k&0x040)>>2)+((k&0x080)>>4)+((k&0x100)>>6)+((k&0x200)>>8)+((k&0x400)>>10);
	return j;
}

int FindIndex10bit(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;
}

int FindIndex9bit(int k){
	int j;
	k=k&(0x1ff);
	j=((k&0x001)<<8)+((k&0x002)<<6)+((k&0x004)<<4)+((k&0x008)<<2)+((k&0x010))+((k&0x020)>>2)+((k&0x040)>>4)+((k&0x080)>>6)+((k&0x100)>>8);
	return j;
}

uchar FindIndex8bit(uchar k){
	uchar j;
	k=k&(0xff);
	j=((k&0x001)<<7)+((k&0x002)<<5)+((k&0x004)<<3)+((k&0x008)<<1)+((k&0x010)>>1)+((k&0x020)>>3)+((k&0x040)>>5)+((k&0x080)>>7);
	return j;
}

uchar FindIndex7bit(uchar k){
	uchar j;
	k=k&(0x7f);
	j=((k&0x001)<<6)+((k&0x002)<<4)+((k&0x004)<<2)+((k&0x008))+((k&0x010)>>2)+((k&0x020)>>4)+((k&0x040)>>6);
	return j;
}

uchar FindIndex6bit(uchar k){
	uchar j;
	k=k&(0x7f);
	j=((k&0x001)<<5)+((k&0x002)<<3)+((k&0x004)<<1)+((k&0x008)>>1)+((k&0x010)>>3)+((k&0x020)>>5);
	return j;
}

//反转下标
void ReverseIndex(COMPLEX *x){
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){					//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[N/2],XkWn;				//wn为旋转因子数组，XkWn为临时存储当前蝶形的旋转因子的临时变量			
	for (k=0;k<N/2;k++){				//初始化wn数组
			wn[k].real=cos(2*3.14159*k/N);
			wn[k].imag=-sin(2*3.14159*k/N);
		}
	ReverseIndex(x);					//输入反序
	for(k=N;k>1;k=(k>>1)){				//第一个循环，代表log2(k)阶的变换
		SubBlockNum=k>>1;				//子块个数为所做点数的一半
		SubBlockStep=SubBlockStep<<1;	//子块间同等地位的元素间隔以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值
		}
	}

}

#endif