//**************************************************
//           128點浮點FFT子程序(Mega128)
//                Passion_FFT
//                 2007.8.18
//**************************************************
#include <avr/io.h>
#include <math.h>
//**************************************************
float FFT_Re[128]={82,88,108,140,164,177,169,148,114,92,81,90,111,141,165,177,169,144,115,91,82,89,112,141,
166,178,168,144,114,90,80,91,112,143,167,178,166,144,112,89,80,91,112,143,168,177,168,142,111,89,82,91,115,
144,168,177,166,142,110,89,81,92,115,146,168,178,165,139,109,88,81,93,117,145,169,177,165,139,109,86,82,93,
117,148,170,176,164,138,106,88,81,94,118,149,170,177,163,137,106,87,80,94,121,149,171,177,163,106,86,82,96,
120,150,171,178,161,136,104,85,82,98,120,152,172,175,161,134,104};    //輸入序列實部，初值為20KHz正弦波AD采樣得到 
float FFT_Im[128]={0.0};     //虛部分存，因為是序列變換因此初值為零
float Wr_Buff[64]={1,0.9988,0.9952,0.9892,0.9808,0.9700,0.9569,0.9415,0.9239,0.9040,
0.8819,0.8577,0.8315,0.8032,0.7730,0.7410,0.7071,0.6716,0.6344,0.5957,0.5556,0.5141,
0.4714,0.4276,0.3827,0.3369,0.2903,0.2430,0.1951,0.1467,0.0980,0.0491,0,-0.0491,
-0.0980,-0.1467,-0.1951,-0.2430,-0.2903,-0.3369,-0.3827,-0.4276,-0.4714,-0.5141,
-0.5556,-0.5957,-0.6344,-0.6716,-0.7071,-0.7410,-0.7730,-0.8032,-0.8315,-0.8577,
-0.8819,-0.9040,-0.9239,-0.9415,-0.9569,-0.9700,-0.9808,-0.9892,-0.9952,-0.9988}; //旋轉因子實部
float Wi_Buff[64]={0,-0.0491,-0.0980,-0.1467,-0.1951,-0.2430,-0.2903,-0.3369,-0.3827,
-0.4276,-0.4714,-0.5141,-0.5556,-0.5957,-0.6344,-0.6716,-0.7071,-0.7410,-0.7730,
-0.8032,-0.8315,-0.8577,-0.8819,-0.9040,-0.9239,-0.9415,-0.9569,-0.9700,-0.9808,
-0.9892,-0.9952,-0.9988, -1,-0.9988,-0.9952,-0.9892,-0.9808,-0.9700,-0.9569,-0.9415,
-0.9239,-0.9040,-0.8819,-0.8577,-0.8315,-0.8032,-0.7730,-0.7410,-0.7071,-0.6716,
-0.6344,-0.5957,-0.5556,-0.5141,-0.4714,-0.4276,-0.3827,-0.3369,-0.2903,-0.2430,
-0.1951,-0.1467,-0.0980,-0.0491};         //旋轉因子虛部
float Window[128]={0.0806,0.0822,0.0850,0.0888,0.0938,0.0998,0.1069,0.1150,0.1242,
0.1343,0.1454,0.1575,0.1705,0.1844,0.1992,0.2147,0.2311,0.2482,0.2660,0.2844,0.3035,
0.3232,0.3433,0.3640,0.3850,0.4065,0.4282,0.4503,0.4725,0.4949,0.5174,0.5400,0.5626,
0.5851,0.6075,0.6297,0.6518,0.6735,0.6950,0.7160,0.7367,0.7568,0.7765,0.7956,0.8140,
0.8318,0.8489,0.8653,0.8808,0.8956,0.9095,0.9225,0.9346,0.9457,0.9558,0.9650,0.9731,
0.9802,0.9862,0.9912,0.9950,0.9978,0.9994, 1,0.9994,0.9978,0.9950,0.9912,0.9862,
0.9802,0.9731,0.9650,0.9558,0.9457,0.9346,0.9225,0.9095,0.8956,0.8808,0.8653,0.8489,
0.8318,0.8140,0.7956,0.7765,0.7568,0.7367,0.7160,0.6950,0.6735,0.6518,0.6297,0.6075,
0.5851,0.5626,0.5400,0.5174,0.4949,0.4725,0.4503,0.4282,0.4065,0.3850,0.3640,0.3433,
0.3232,0.3035,0.2844,0.2660,0.2482,0.2311,0.2147,0.1992,0.1844,0.1705,0.1575,0.1454,
0.1343,0.1242,0.1150,0.1069,0.0998,0.0938,0.0888,0.0850,0.0822,0.0806,0.0800};   //窗函數，改進漢寧窗
float Res[128];             //轉換結果
float Result1_r;
float Result1_i;
float Result2_r;
float Result2_i;
//*************************************************
void FFT_Proc(void);
//*************************************************
int main(void)
{
    unsigned char i;
	for(i=0;i<128;i++)
	    FFT_Re[i]*=Window[i];
    FFT_Proc();
	while(1);
	return 1;
}
//************************************************* 
/*------------------------------------------------
FFT_Proc ：計算128點FFT
    信息 ：標準C浮點運算，可移植
------------------------------------------------*/
void FFT_Proc(void)
{
    unsigned char i;
    unsigned char j;
    unsigned char a;
    unsigned char b;
    unsigned char c;
    unsigned char k=0;
    unsigned char temp=0;
    unsigned char revers=0;
    float FFT_Re1[128];
    for(i=0;i<128;i++)
	    FFT_Re1[i]=FFT_Re[i];
    for(i=0;i<128;i++)
    {
        for(j=0;j<7;j++)
        { 
            a=(i&(0x01<<j));
            b=(6-(j<<1));
            c=6-2*(6-j);
            if(j<3)
	       revers+=a<<b;
            else if(j>3)
	       revers+=a>>c;
            else
               revers+=a;
        }
        FFT_Re[i]=FFT_Re1[revers];	    
        revers=0;
    }                               //對128點輸入序列進行位倒序
    for(i=0;i<7;i++)                //第一重循環控制蝶形的級數
    {
        for(j=0;(j+(1<<i))<128;j++) //控制每級的蝶形和旋轉因子
        {
             if(k<64)
             {
                c=k;
                k+=(1<<(6-i));
             }
             else
             k=0;
             if(temp<(1<<i))
             {
                 temp++;
                 a=j;
                 b=(j+(1<<i));
                 Result1_r=FFT_Re[a]+Wr_Buff[c]*FFT_Re[b]-Wi_Buff[c]*FFT_Im[b];
	             Result1_i=FFT_Im[a]+Wi_Buff[c]*FFT_Re[b]+Wr_Buff[c]*FFT_Im[b];
	             Result2_r=FFT_Re[a]-Wr_Buff[c]*FFT_Re[b]+Wi_Buff[c]*FFT_Im[b];
	             Result2_i=FFT_Im[a]-Wi_Buff[c]*FFT_Re[b]-Wr_Buff[c]*FFT_Im[b];  //蝶形單元
                 FFT_Re[a]=Result1_r;
                 FFT_Im[a]=Result1_i;
                 FFT_Re[b]=Result2_r;
                 FFT_Im[b]=Result2_i;                                              
             }
             else
             {
                 j=(j+(1<<i)-1);
                 temp=0;
             }
		    
         }
         k=0;
         temp=0;
     }
     for(i=0;i<128;i++)
     Res[i]=sqrt(FFT_Re[i]*FFT_Re[i]+FFT_Im[i]*FFT_Im[i]);                      //結果取模
}
//**************************************************************************