FFT算法实现

关于FFT算法的原理这里就不多说了,具体参考有关书籍。

DFT与FFT运算量的比较

N点DFT的运算量

 

复数乘法

复数加法

一个X(k)

N

N-1

N个X(k)(N点DFT)

N*N

N(N-1)

 

N点FFT的运算量

 

复数乘法

复数加法

N个X(k)

(N/2)*log2N

N*log2N

 

如何用计算机程序实现FFT呢

参考《数字信号处理》(西电版)P115的蝶形运算公式,可以推导出用计算机程序实现的公式。由于word中输入公式太麻烦了,这里先将推导过程写在纸上,然后拍照上传到这里。



 根据上面的推导结果,就可以用程序来实现FFT了。看下面的例程。

#include "math.h"
#define PI	3.14159
#define N   128

//real:实数
//imaginary:虚数
#pragma DATA_SECTION(Input,"ExRamMem");//这是DSP处理器(CCSv4开发环境)中的一种操作内存的方法
int Input[N] = {0.0};//用于存放待处理数据
#pragma DATA_SECTION(SinTab,"ExRamMem");
float SinTab[N] = {0.0};//用于FFT运算的正弦表
#pragma DATA_SECTION(CosTab,"ExRamMem");
float CosTab[N] = {0.0};//用于FFT运算的余弦表
#pragma DATA_SECTION(Real,"ExRamMem");
float Real[N] = {0.0};//待处理数据的实部
#pragma DATA_SECTION(Imag,"ExRamMem");
float Imag[N] = {0.0};//待处理数据的虚部
#pragma DATA_SECTION(Mag,"ExRamMem");
float Mag[N] = {0.0};//FFT的幅度谱

//函数功能:生成一个正弦波信号
void GenSin(void)
{
	int i;
	
	for ( i=0;i<N;i++ )
	{
		INPUT[i]=sin(PI*2*i/N*3)*1024;
	}
}

//函数功能:计算用于FFT运算的正、余弦表。
//(me)注意:旋转因子和正、余弦表不是同一个概念。
void TwiddleCmpt(void)
{
	int i;
	
	for ( i=0; i<N; i++ )
	{
		SinTab[i] = sin(PI*2*i/N);
		CosTab[i] = cos(PI*2*i/N);
	}
}

//输入参数:real:待处理的数据的实部
//		  imag:待处理的数据的虚部
//返回参数:暂时就返回OK
//函数功能:计算实部real和虚部imag的N点FFT。
//注意:因为FFT的运算为了节约内存,而将下一级的运算结果直接覆盖上一级,所以下面的程序中要采取temp_R、temp_I、temp_R_kb暂时存放一下。
int FFT(float real[N],float imag[N])
{
	//调试发现:定义long型变量与定义int型变量时的计算速度居然不一样,定义long型变量的计算速度要慢一点,即使是参与运算的相关变量都是long型的。但是没办法当FFT点数较大时,int型变量的范围是不够用的。
	long i,j=0,k,b,p,Nv2,Nm1,bm1,bx2;	//Nv2表示N/2,Nm2表示N-1,bm1表示b-1,bx2表示b乘2,PIx2表示PI乘2,定义这五个数是为了减少运算时间。
	float temp,temp_R,temp_I,temp_R_kb;
	long L,M;
	
//	f=N;
//	for(m=1;(f=f/2)!=1;m++)		//判断N是否能够被2整除,判断N是2的几次幂。
//	{}
	Nv2=N/2;
	Nm1=N-1;
	M=log(N)/log(2);
	
	//这里的倒序采用雷德算法。
	//=================== following code invert sequence ===================
	for(i=0; i<Nm1; i++)
	{
		if(i < j)
		{
			temp = real[j];
			real[j] = real[i];
			real[i] = temp;
		}
		k = Nv2;
		while(k <= j)
		{
			j = j - k;
			k = k/2;
		}
		j = j + k;
	}
	
	//===================following code FFT ===================
	for (L=1; L<=M; L++ )//for(1):控制运算级数
	{
		b = 1;
		i = L - 1;
		while (i > 0)
		{
			b=b*2; i--;
		}	// b= 2^(L-1),蝶尖距
		
		bm1 = b-1;	//在进行第二层循环前,先将bm1和bx2算出来。
		bx2 = b*2;
		for (j=0; j<=bm1; j++ )//for(2):同一级中的蝶群数。
		{
			p=1; i=M-L;
			while (i > 0)
			{
				p = p*2; i--;
			}
			p = p*j;	// p=pow(2,M-L)*j;
			
			for (k=j; k<N; k=k+bx2)//for(3):一个蝶群中有多少只蝴蝶
			{
				temp_R = real[k]; temp_I = imag[k]; temp_R_kb = real[k+b];
				real[k]  = real[k] + real[k+b] * CosTab[p] + imag[k+b] * SinTab[p];
				
				imag[k]  = imag[k] - real[k+b] * SinTab[p] + real[k+b] * CosTab[p];
				
				real[k+b]= temp_R  - real[k+b] * CosTab[p] - real[k+b] * SinTab[p];	//因为此时的real[k]已被上一条语句覆盖,但这里还是想用之前的real[k],所以才事先将real[k]存储在temp_R中,temp_I和temp_R_kb与此类似。
				
				imag[k+b]= temp_I  + temp_R_kb * SinTab[p] - real[k+b] * CosTab[p];
			} // END for (3) 
		} // END for (2) 
	} // END for (1) 
	for ( i=0; i<Nv2; i++ )//计算幅度谱
	{
		Mag[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);//FFT的输出是由实部和虚部组成,要想看频谱,还求实部和虚部平方和再开根号,这样得到就才是幅度谱。
	}
	asm("	NOP");
	return(OK);
}

int fft(void)
{
	Uint16 	i;
   	GenSin();
   	TwiddleCmpt();
   	for(i=0; i<N; i++)
   	{
   		Real[i] = Input[i];//把采集到的N个数据付给FFT的实部。
   		Imag[i] = 0.0;//虚部置0
   		Mag[i] = 0.0;//幅度置0
   	}
   	FFT(Real, Imag);
	return(OK);
}

正弦波生成函数GenSin()生成波形如下

经过fft处理后得到的幅度谱如下

注意:

1.  这里得到的幅度谱,其横坐标只是下标值,若想将横坐标转换为频率,还需要进行计算,公式是:f=n*Fs/N,其中n是对应处理结果Mag的下标值,Fs是采样率,N是fft点数,f即是n点对应的频率。

2.  如下图这个处理结果,给人的感觉是只有Mag[3]的值较大,其他都是0,其实有很多点的值不是等于0的,只是他们的值相对Mag[3]太小,所以显得他们是0。

3.  这里的信号源是交流信号,即不含直流分量,所Mag的第0点的值是0。而一般情况下AD的采样结果都是正值,也就是采样信号中会含有直流分量,这时候在第0点也会有值,只不过该点对应的频率是0。

CCSv4软件自带的频谱分析工具分析的fft处理结果如下


关于倒序

以N=128(2^7)为例,给出倒序前后的序号。

倒序前序号

倒序后序号

十进制

二进制

二进制

十进制

0

000 0 000

000 0 000

0

1

000 0 001

100 0 000

64

2

000 0 010

010 0 000

32

3

000 0 011

110 0 000

96

4

000 0 100

001 0 000

16

……

……

……

……

126

111 1 110

011 1 111

63

127

111 1 111

111 1 111

127

 

倒序除了上面的雷德算法也外,可以采用下面这种更为直接的方法进行,不过这种方法的效率显然没有雷德算法高。

{
	int x0,x1,x2,x3,x4,x5,x6,xx;
		
	//倒序
        //=================== following code invert sequence ===================//
        for(i=0; i<128; i++)
        {
	    x0=x1=x2=x3=x4=x5=x6=0;
            x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;
            xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
            imag[xx] = real[i];
        }
        for(i=0; i<128; i++ )
        {
            real[i] = imag[i];
            imag[i] = 0;
        }
}

 


  • 3
    点赞
  • 13
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值