DSP 基2 FFT 算法详细代码带注释

#include “DSP2833x_Device.h“
#include “DSP2833x_Examples.h“
#include  


/*********************************************************
傅里叶变换:           （ ∫为积分符号，从负无穷到正无穷）

		F（w）=         ∫f（t）exp（-i*w*t）dt
		f（t）= （1/2*pi）∫F（w）exp（i*w*t）dw

		F（s）= ∫f（t）exp（-i*2*pi*s*t）dt
		f（t）= ∫F（s）exp（i*2*pi*s*t）ds
 
离散傅里叶变换（DFT）   （ ∑为连加符  WN=exp（-i*2*pi/N））

		F（sj） =     ∑k Ak*（WN）^（jk）
		Ak = （1/N）* j∑ F（sj）*（WN）^（-jk）

		Ak=[f（t（k-N）） + f（tk）] * delta_t

		delta_t为采样时域间隔tk=k*delta_t                               
        delta_s为F（sj）频谱对应频率sj的间隔sj=j*delta_s

        上式中Ak并不是原时序序列，故频谱也不是频谱序列，再经变换，可得

		F（sj）=（1/N）*∑k f（tk）WN^（jk） 
		f（tk）=      ∑j F（sj）WN^（-jk） 
       
	    这里f（tk）即为原输入时间序列，F（sj）为对应频谱
	
		为了计算方便，假设如下

        F（j）= ∑k f（tk）WN^（jk）
        计算完F（j）后，再乘以（1/N）即得到F（sj）  
       

*********************************************************/

 
#define pi 3.141593       // float小数点后6位
#define N_SIGN 512        // 信号合成和采样点数直接影响delta_t
#define N_COMPOUND 10     // 合成谐波的数量
#define TIMES_T    10     // 完整信号周期个数与合成信号频率无关，只会影响采样速率和视窗
                          // 周期越数越多，采样速度越低，采样时域间隔越大。
						  // 当采样速率 f 不满足奈奎斯特定律（f >= 2 * max frequence of signal）时，无法还原

int N_FFT=N_SIGN;         // FFT点数（频域和时域点数一样）
float input[N_SIGN];      // 输入的实信号序列
float output[N_SIGN];     // 输出的实FFT幅值（复数的模）


//合成信号的正弦频率分布		
int s_compound[10]={135791113151719};  
//合成信号的正弦幅值分布     
float Fs_compound[10]={（float）10/（float）1 （float）10/（float）3 （float）10/（float）5
					   （float）10/（float）7 （float）10/（float）9 （float）10/（float）11
					   （float）10/（float）13（float）10/（float）15（float）10/（float）17
					   （float）10/（float）19};   


struct Complex		      // 定义复数结构体
{
   	float realimag;
};

struct Complex WNP;                   // 定义蝶形旋转因子
struct Complex WNP_XJB;               // 存放旋转因子与F（J+B）（即蝶形的下方入口）的乘积
struct Complex input_complex[N_SIGN]; // 采样输入的实数转化为复数






/************     复乘函数  ************/
struct Complex MUL_Complex（struct Complex astruct Complex b）
{
   	struct Complex c;
   	c.real=a.real*b.real-a.imag*b.imag;
   	c.imag=a.real*b.imag+a.imag*b.real;
   	return（c）;
}





/********************************
功能：计算复数的模
形参：*sample指向需要取模的复数结构体
      N为取模点数
	  *output存放取模后数值的数组
*********************************/
void ModelComplex（struct Complex *sampleint Nfloat *output）
{
   	int i;
   	for（i=0;i    {
		output[i]=（float）0;
		// 乘2因为正负谱的需要，乘（1/N）是因为要从F（j）求F（sj）F（sj）才为真实的频谱
    	output[i]=sqrt（sample[i].real*sample[i].real+sample[i].imag*sample[i].imag）*2/N;
	}
}








/****************   Cyrus Cheung的基2FFT算法解析   *************
基2FFT原理:
	对于N=2^M个点的FFT，共有L=M级运算，每级有2^（M-1）即N/2个蝶形运算组成；
	旋转因子为WNPWNP=exp（-i*2*pi*p/N）;
	P=J*2^（M-L）J=01......2^（L-1）-1，即每级有 2^（L-1）个不同的旋转因子；

	对于同一P，参与本级蝶形运算的次数为2^M-L，B=2^（L-1）为同一蝶形的两输入点间的距离，即距离B等于不同旋转因子的个数；
	每级第一次调用WNP的蝶形结第一个输入节点的位置为J第二个输入节点的位置为J+B=J+2^（L-1）；
	每级第二次调用WNP的蝶形结第一个输入节点的位置为J+2^L第