real_sequence_fft.c

这是一个对实数序列进行FFT变换求频谱的C语言程序

/*程序说明：
  如果序列x(n)是实数，那么其傅立叶变换X（k）一般是复数，但其实部是偶对称，虚部是奇对称，即
X（k）具有如下共轭对称性：X（0）和X（N/2）都是实数，且有X（k）=X*(N-k),  1<=k<=N/2-1
在计算离散傅立叶变换时，利用这种共轭对称性，我们就可以不必计算与存储X(k)(N/2+1<=k<=N-1)
以及X（0）和X（N/2）的虚部，而仅需计算X（0）到X（N/2）即可。此处我们选择的是计算
X（0）到X（N/4）和X（N/2）到X(3N/4),这样做可以恰好利用复序列FFT算法的前（N/4）+1个复数蝶形。
这就是按时间抽取的基2实序列FFT算法，比复序列FFT算法大约可减少一半的运算量和存储量。

*/

/*形参说明:
  x----双精度实型一维数组，长度为n。开始时存放要变换的实数据，最后存放变换结果的前n/2+1个值，其存储顺序为
		［Re(0),Re(1),...,Re(n/2),Im(n/2-1),...,Im(1)］。其中Re(0)=X(0),Re(n/2)=X(n/2)(X(1)的实部为x(1),虚部为x(n-1),
		可推出结果的实部与虚部的数据位置下标之和为n)。
		根据X（k）的共轭对称性，很容易写出后半部分的值。
		程序为：
		printf("%10.7f+J%10.7f\n",x[n/2-1],-x[n/2+1]);//这是X[n/2＋1]的值。
			因X[n/2+1]关于共偶对称的是X[n/2-1]，而X[n/2-1]的实部存在x[n/2－1]，虚部存在x[n/2+1]从而对实部取偶，虚部取奇
			得X[n/2+1]
		for(i=2;i<n/2;i+=2)
			{	printf("%10.7f+J%10.7f\n",x[n/2-i],-x[n/2+i]);//这是X[n/2＋i]的值。
				printf("%10.7f+J%10.7f\n",x[n/2-i-1],-x[n/2+i+1]);//这是X[n/2＋i+1]的值。
			}
  n-----整型变量。数据长度，必须是2的整数次幂，即n=2*m
*/
#include "stdlib.h"
#include "math.h"
void rfft(x,n)
int n;
float      x[];
{
int i,j,k,m,i1,i2,i3,i4,n1,n2,n4;
float       a,e,cc,ss,xt,t1,t2;
puts("进行FFT变换......\n");
for(j=1,i=1;i<16;i++)
		{
			m=i;
			j=2*j;
			if(j==n) break;
		}
n1=n-1;
	for(j=0,i=0;i<n1;i++)
		{
			if(i<j)
			{
				xt=x[j];
				x[j]=x[i];
				x[i]=xt;
			}
			k=n/2;
			while(k<(j+1))
				{
					j=j-k;
					k=k/2;
			
					}
			j=j+k;
		 }
	for(i=0;i<n;i+=2)
		{xt=x[i];
		x[i]=xt+x[i+1];	
		x[i+1]=xt-x[i+1];	
		}
	n2=1;
	for(k=2;k<=m;k++)
		{
			n4=n2;
			n2=2*n4;
			n1=2*n2;
			e=6.28318530718/n1;
			for(i=0;i<n;i+=n1)
				{	xt=x[i];
					x[i]=xt+x[i+n2];
					x[i+n2]=xt-x[i+n2];
					x[i+n2+n4]=-x[i+n2+n4];
					a=e;
						for(j=1;j<=(n4-1);j++)
							{i1=i+j;
							i2=i-j+n2;
							i3=i+j+n2;
							i4=i-j+n1;
							cc=cos(a);
							ss=sin(a);
							a=a+e;
							t1=cc*x[i3]+ss*x[i4];
							t2=ss*x[i3]-cc*x[i4];
							x[i4]=x[i2]-t2;
							x[i3]=-x[i2]-t2;
							x[i2]=x[i1]-t1;
						    x[i1]=x[i1]+t1;
						
							}
				
					}
		}
puts("FFT变换完毕。\n");
}