fft6.c

fft和ifft的c程序代码

#include <math.h>
#include <stdio.h>
#define PI 3.14159265
#define NN 32
static float x[32]={ -1, 1, 1,-1,-1, 1,-1,
-1,-1, 1,-1, 1, 1, 1, 1,
-1,-1, 1,-1, 1, 1,-1, -1,
1,-1, 1,-1,-1,-1, 1, 1, 1,};

struct compx {float real,imag;};    /*定义一个复数结构*/
struct compx EE(struct compx,struct compx); /*定义复数相乘结构*/
void ifft(float *);
void ditfft(struct compx * , int);          /*FFT函数定义*/

struct compx EE(struct compx b1,struct compx b2)
{
    struct compx b3;
    b3.real=b1.real*b2.real-b1.imag*b2.imag;
    b3.imag=b1.real*b2.imag+b1.imag*b2.real;
    return(b3);
}

void ifft(float *a)
{
    int p;
    struct compx s[NN+1];
    for(p = 1; p <= NN; p++)
    {
        s[p].real = a[2*p-2];		/*a的数据格式是偶数位为实部，奇数位为虚部*/
        s[p].imag = -a[2*p-1];		/*这里取模，做ifft*/
    }

    ditfft(s, NN);
	
    for(p = 1; p <= NN; p++)
    {
        a[2*p-2] = s[p].real/NN;
        a[2*p-1] = -s[p].imag/NN;	/*再取一次模，并除已NN，做ifft*/
    }
}

void ditfft(struct compx *xin,int N)
{
	int le,lei,ip;  /*le DFT组间距离 lei 每组DFT内的间距*/
	float pi;
    int f,m,nv2,nml,i,k,j=1,l;
    struct compx v,w,t;
    nv2 = N/2;
    f = N;
    for(m = 1; (f = f/2) != 1; m++) {;}
    nml = N-1;
	
    /*变址运算*/
    for(i = 1; i <= nml; i++)
    {
        if(i<j) {
			t=xin[j];
			xin[j]=xin[i];
			xin[i]=t;
		}
        k=nv2;
        while(k<j)  {
			j=j-k;
			k=k/2;
		}
        j=j+k;
    }
	
    /*fft*/
	for(l=1;l<=m;l++)
	{
		le=pow(2,l);
		lei=le/2;
		pi=PI;
		/*pi=3.14159265;*/
		v.real=1.0;
		v.imag=0.0;
		w.real=cos(pi/lei);
		w.imag=-sin(pi/lei);
		for(j=1;j<=lei;j++) /*j 每组DFT的个数*/
		{
			for(i=j;i<=N;i=i+le)    /*i i各个组相对应的位置上的DFT运算*/
			{
				ip=i+lei;
				t=EE(xin[ip],v);
				xin[ip].real=xin[i].real-t.real;
				xin[ip].imag=xin[i].imag-t.imag;
				xin[i].real=xin[i].real+t.real;
				xin[i].imag=xin[i].imag+t.imag;
			}
			v=EE(v,w);
		}
	}
}	

main()
{
	int i, j = 0;  
	float a[NN*2];   
	for(i = 0; i < NN * 2; i += 2)  
	{
		a[i] = x[j];
		j++;
		//a[i]=1.0/(((i-NN-1.0)/NN)*((i-NN-1.0)/NN)+1.0);
		a[i + 1] = 0.0; 
		printf("%f	+i	%f	", a[i],  a[i + 1]);
		if((i+1) % 2 == 0)	
			printf("\n");
	}
	printf("\n");
	ifft(a);
	printf("\n");
	for(i = 0; i < NN; i += 2)
	{
		printf("%f	+i	%f	", a[i], a[i + 1]);
		if((i+1) % 2 == 0)	
			printf("\n");
	}
}