fd_sp.cpp

时间抽取(DIT)基2FFT算法: 输入:先输入数组长度(8,16,32...),再输入数组内容. 输出:FFT变换结果.

/***********       FFT     *************/

#include     <stdio.h>
#include     <math.h>
#include     <stdlib.h>

#define      N    25            
#define      Pi   3.14159    
#define      M    3              

struct  complex        
{
  double x;
  double y;
} ;
typedef struct complex COMPLEX          ;
int            L                        ;

COMPLEX     mul(COMPLEX a,COMPLEX b)    ;            /* 复数相乘       */
COMPLEX     add(COMPLEX a,COMPLEX b)    ;            /* 复数相加       */
COMPLEX     sub(COMPLEX a,COMPLEX b)    ;            /* 复数相减       */
int         back_n(int)                 ;            /* 对下标实现倒序 */
COMPLEX *   back_x(double a[],int)      ;            /* 实现倒序列     */
COMPLEX *   comput_x(COMPLEX X[],int)   ;            /* 计算蝶形       */
void        output(COMPLEX X[],int)     ;            /* 输出           */


int main(void)
{
	int n, i;
	double	in[N];
	COMPLEX  *q;
	printf("Please input the n:");       
	scanf ("%d$\n",&n);
    L=int(ceil(log(n)/log(2)));                      
    printf("Now,please input X[n].\n");
	for (i=0;i<n;++i)
	    scanf("%lf$",&in[i]);
	q=back_x(in,n);
	q=comput_x(q,n);
	output(q,n);	
    return 0;
}


int   back_n(int u)
{
	int k[M], i, sum=0, temp;
	for (i=M-1; i>=0 ;--i)
	{
		k[i]=u%2;
		u=u/2;
	}

    for (i=0;i<M/2;i++)
	{
		temp=k[i];
		k[i]=k[M-1-i];
        k[M-1-i]=temp;
	}

	for (i=0;i<M;i++)
    {
		sum*=2;
		sum+=k[i];
	}
	return sum;
}


COMPLEX * back_x(double a[],int n)                    
{	
	int i, j;
	double p[N];
	static COMPLEX TEM[N];
	for (i=0;i<n;++i)                 
	{
		j=back_n(i);
		p[i]=a[j];
	}
    for (i=0;i<n;++i)
	{
		TEM[i].x=p[i];
		TEM[i].y=0;
	}
	printf("\nThe backward is as follows:\n");	
	for (i=0;i<n;++i)
	printf("%lf\n",p[i]);
	return TEM;
}

COMPLEX mul(COMPLEX a,COMPLEX b)
{
	COMPLEX c;
	c.x=a.x*b.x-a.y*b.y;
	c.y=a.x*b.y+a.y*b.x;
	return c;
}

COMPLEX add(COMPLEX a,COMPLEX b)
{
	COMPLEX c  ;
	c.x=a.x+b.x;
	c.y=a.y+b.y;
	return c   ;
}


COMPLEX sub(COMPLEX a,COMPLEX b)
{
	COMPLEX c  ;
	c.x=a.x-b.x;
	c.y=a.y-b.y;
	return    c;
}

COMPLEX * comput_x(COMPLEX X[],int n)
{
	int le, k, x, bty, r, m, g=7;                     /* le为级，bty为组，k为蝶中X[n]变化 */
	COMPLEX E, x1, x2, mult;
	for (le=0;le< L;++le)                                                  
	{
		k=0;
		x=1<<(le);
		
		for(bty=0; bty<(1<<(L-le-1));                  
		    bty++,x+=(1<<(le+1)),k+=(1<<(le)))
					
			for ( ;k<x; k++)                              
			{
				r=(k<<(L-le-1))&g;         
				m=1<<(le);
				E.x=cos((2*Pi/n)*r);                      
				E.y=-sin((2*Pi/n)*r);
				x1=X[k];
				x2=X[k+m];
				mult=mul(x2,E);
                X[k]=add(x1,mult);
				X[k+m]=sub(x1,mult);
				                                
			}
	
	}
	return X;
}

void output(COMPLEX X[],int n)
{
	int i;
	printf("The result X[n] is as follows:\n");
    for( i=0; i<n; ++i)
		if(X[i].y==0)
	        printf("%lf\n",X[i].x);
		else if(X[i].y<0)
			printf("%lf %lfj\n",X[i].x,X[i].y);
		else
			printf("%lf +%lfj\n",X[i].x,X[i].y);
}