FFT算法实现

最新推荐文章于 2024-05-07 20:12:00 发布

取取经

最新推荐文章于 2024-05-07 20:12:00 发布

阅读量3.2k

点赞数 3

分类专栏：嵌入式文章标签： fft 算法

本文链接：https://blog.csdn.net/yahstudio/article/details/38711049

版权

嵌入式专栏收录该内容

5 篇文章 0 订阅

订阅专栏

关于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;
        }
}