fft.cpp

32位嵌入式系统实现

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

#define PI 3.14159265358979 
#define N  1024

struct comp_float
{
	float real;
	float img;
};

struct comp_int
{
	int real;
	int img;
};

struct comp_float  float_in[N], out[N];
struct comp_int  hex_in[N],hex_out[N];
struct comp_int  *w,in[N];

int size_x = 1024;
int m = 10;        /*级数*/

int beginFFT()
{
	int comput_m(int size_x);	
	void fft();
	void initw();
	void change();
	int i=0;
    FILE *fp_input_r,*fp_input_i,*fp_output_r,*fp_output_i;
    fp_input_r=fopen("RHr_infft.txt","r");
	fp_input_i=fopen("RHi_infft.txt","r");
	fp_output_r=fopen("RKr_fftout.txt","w");
	fp_output_i=fopen("RKi_fftout.txt","w");
    /*printf("Please input the size of x:\n");
    scanf("%d",&size_x);
    m=comput_m(size_x);*/	
	for(i=0;i<size_x;i++)
    {  	   
	  fscanf(fp_input_r,"%x",&hex_in[i].real);
	  fscanf(fp_input_i,"%x",&hex_in[i].img);
	  in[i].real = (int) (hex_in[i].real*4096);
	  in[i].img  = (int) (hex_in[i].img*4096);
    }	
	initw();
	fft();
	for(i=0;i<size_x;i++)
	{
		float_in[i].real = (float) ((float) in[i].real*256);
		float_in[i].img = (float) ((float) in[i].img*256);
		hex_out[i].real = (int) (float_in[i].real);
		hex_out[i].img = (int) (float_in[i].img);
		fprintf(fp_output_r,"%x\n",hex_in[i].real);
		fprintf(fp_output_i,"%x\n",hex_in[i].img);
		//printf("%d %d\n",in[i].real,in[i].img);	
		
	}
	fprintf(fp_output_r,"--------------------------------\n");
	fprintf(fp_output_i,"--------------------------------\n");
	printf("-------------------------------------------\n");
	for(i=0;i<size_x;i++)
	{				
		printf("%f %f\n",float_in[i].real,float_in[i].img);
	}

	fclose(fp_input_r);
	fclose(fp_input_i);
    fclose(fp_output_r);
	fclose(fp_output_i);
    return 0;
}


void initw()
{
   int i=0;
   w=(comp_int *)calloc( size_x/2, sizeof(comp_int) );
   for(i=0;i<size_x/2;i++)
   {   
       w[i].real = (int)(cos(2*(PI)/size_x*i)*32768);
       w[i].img = (int)(-1*sin(2*(PI)/size_x*i)*32768);
	   
   }         
}

void change(comp_int *z)
{  comp_int temp;
   int i=0,j=0,k=0,t;
   for(i=0;i<size_x;i++)
   {
     k=i;j=0;
     t=m;
     while(t--)
	 {
        j=j<<1;
        j|=(k & 1);
        k=k>>1;
     }
     if(j>i)
	 {
        temp=z[i];
        z[i]=z[j];
        z[j]=temp;
     }
   }
}

int comput_m(int n)
{
	int i=0,number;
	number=n;
	do
	{
		number>>=1;
		i++;
	}while(number!=1);
	return i;
}

void fft()
{
	void func_die(int ,int ,int , int ); 
	change(in);
	int  i=0,j=0,k=0;
	int  L,c_num,beg;   //组数,组内元素个数,组内计算时的起始位置

	for(i=0;i<m;i++)
	{
		L = pow(2,m-i-1);		    	

		for(j=0;j<L;j++)
		{
			c_num = pow(2,i+1);
			beg = j*c_num;
			for(k=0;k<c_num/2;k++)
			{
				
				func_die(beg,c_num,k,i);
			}
		}
	}
}

void func_die(int beg,int c_num,int k,int i)
{
	void add(comp_int ,comp_int ,comp_int *);
	void sub(comp_int ,comp_int ,comp_int *);
	void mul(comp_int ,comp_int ,comp_int *);

	int temp_real =0 ,temp_img = 0;
	struct comp_int temp1,up,down;
	
	int a,b,p;
	a = beg + k;
	b = beg + k + c_num/2;
	p = k*pow(2,m-1-i);
	mul(in[b], w[p], &up);
	add(in[a], up, &up);
	mul(in[b], w[p], &down);
	sub(in[a], down, &down);
	in[a]=up;
	in[b]=down;

	/*temp_real =  (((w[p].real*in[b].real)>>15)-((w[p].img*in[b].img)>>15))>>1;
	temp_img  =  (((w[p].real*in[b].img)>>15) + ((w[p].img*in[b].real)>>15))>>1;
	temp1.real =  (in[a].real + temp_real)>>1;
	temp1.img  =  (in[a].img +temp_img)>>1;
	temp_real =  (((w[p].real*in[b].real)>>15)-((w[p].img*in[b].img)>>15))>>1;
	temp_img  = (((w[p].real*in[b].img)>>15) + ((w[p].img*in[b].real)>>15))>>1;
	in[b].real = (in[a].real - temp_real)>>1;
	in[b].img  = (in[a].img - temp_img)>>1;
	in[a].real =temp1.real;
	in[a].img =temp1.img;*/	
	
}


void add(comp_int a,comp_int b,comp_int *c)//(16,15)+(16,15)->(16,14)
{
	c->real=((a.real+b.real)>>1);
	c->img=((a.img+b.img)>>1);
}

void mul(comp_int a,comp_int b,comp_int *c)//(16,15)*(16,15)->(16,15)
{
	c->real=((a.real*b.real)>>15) -((a.img*b.img)>>15);
	c->img= ((a.real*b.img)>>15) +((a.img*b.real)>>15);
}

void sub(comp_int a,comp_int b,comp_int *c)//(16,15)+(16,15)->(16,14)
{
	c->real=((a.real-b.real)>>1);
	c->img=((a.img-b.img)>>1);
}