fft.c

快速傅立叶变化的fft算法的C语言实现。

#include "fft.h"

static double cosine[8]=				/*	cosine[k] = cos(PI/(2^k)) */
{
	-1.0, 0.0, 0.7071068, 0.9238795, 0.9807853, 0.9951847, 0.9987955, 0.9996988
};


static  double sine[8]=					/*	sine[k] = sin(PI/(2^k))  */
{
     0.0, 1.0, 0.7071068, 0.3826834, 0.1950903, 0.0980171, 0.0490677, 0.0245412
};

/*
此程序是假设一个周期最多采64点的fft变换子程序,如果需要更多点，可以继续添加系数表。

入口：n－－采样点数
	  a－－a是指针，指向采样点构成的数组首地址，采样点存在数组的复数成员的实部。
	  
出口：a－－由于a是指针，可以改变其指向的储存单元的数值，实现双向传递。因此fft运算后的
		   结果保存在原采样点所在的数组，这也就是原位运算。
*/


void fft(int n, complex *a) 
{
	unsigned char le, le1, ip;
 	unsigned char i,j,k;
 	complex	t, u, w;	/*用于计算角度*/
 	unsigned char l;	/*矩形窗的长度*/
	int n1;

	l = 0;
	n1 = n;
	while ( n1 > 1 )	   /*求矩形窗的长度*/
	{
		n1 >>= 1;
		l ++;
	}

	/**********  码位倒置************/
	/*例：10011变为11001      */
  	for(j = 1, i = 1; i <= n - 1; i++) 
  	{
   		if(i < j) 
   		{
    			t.r = a[j-1].r; 
    			t.i = a[j-1].i;
    			a[j-1].r = a[i-1].r; 
    			a[j-1].i = a[i-1].i;
    			a[i-1].r = t.r; 
    			a[i-1].i = t.i;
   		}
   		
   		k = n / 2;
   		while (k < j) 
   		{
			j -= k;
			k >>= 1;
   		}
   		j += k;
  	}
			
 
	for (k = 1; k <= l; k++)
	{ 

  		le = 1 << k;
  		le1 = le / 2;

  		u.r = 1.0; 
  		u.i = 0.0;


  		w.r = cosine[k-1]; 
  		w.i = sine[k-1];

  		for(j = 1; j <= le1; j++) 
  		{
   			for (i = j; i <= n; i += le) 
   			{
    				ip = i + le1;
    				t.r = a[ip-1].r * u.r - a[ip-1].i * u.i;
    				t.i = a[ip-1].r * u.i + a[ip-1].i * u.r;

    				a[ip-1].r = a[i-1].r - t.r; 
    				a[ip-1].i = a[i-1].i - t.i;
	 
    				a[i-1].r = a[i-1].r + t.r; 
    				a[i-1].i = a[i-1].i + t.i;
   			}

   			t.r = u.r * w.r + u.i * w.i;
   			t.i = -u.r * w.i + u.i * w.r;

   			u.r = t.r; 
   			u.i = t.i;
  		}
 	}
}