Board logo

标题: C6701 的DSP 按频率抽取FFT程序,仿真通过 [打印本页]

作者: 苹果也疯狂    时间: 2014-5-25 18:02     标题: C6701 的DSP 按频率抽取FFT程序,仿真通过

首先说下为什么要FFT吧。主要进行的是时域和频域的变换。比如在设计滤波器的时候需要用到大量的卷积。这时候在时域上做卷积通过FFT转换到频域的时候,只需要对频域做乘加操作。再通过IFFT变换回时域即可。


#include "math.h"


#define PI  
3.14159265358979


#define numlength  


2048


void radix2(int n,short *xy,short *w);


void bitrev_cplx(int *x,short *index,int nx);


void digitrev_index(short *index, int n, int radix);




static short index1[64], w1[numlength];


static short Input[numlength*2];


static int Output[numlength];


static int radix;


int Num;


int i;





void main()


{


radix=2;


Num=numlength;



for (i=0;i


{


w1[2*i] = 32767*(-cos((double)i*2*PI/Num));


w1[2*i+1] = 32767*(-sin((double)i*2*PI/Num));


}



for (i=0;i


{


Input[2*i] = (short)((cos(PI*i/20.0)+cos(PI*i/10.0)+cos(PI*i/5.0))*0x80);


Input[2*i+1] = 0;


}



radix2(Num, Input, w1);


for (i=0; i


{


Output = sqrt(Input[2*i]*Input[2*i]+Input[2*i+1]*Input[2*i+1]);  


}






digitrev_index(index1,Num,radix);


bitrev_cplx(Output, index1, Num);


}  

































//频率抽取FFT算法函数体



void radix2(int Num, short xy[], short w[])


{


short n1,n2,ie,ia,i,j,k,l;


short xt,yt,c,s;


n2 = Num;


ie = 1;


for (k=Num;k > 1;k = (k >> 1


) )


{


n1 = n2;


n2 = n2>>1;


ia = 0;


for (j=0; j < n2; j++)


{


c = w[2*ia];


s = w[2*ia+1];


ia = ia + ie;


for (i=j; i < Num; i += n1)


{


l = i + n2;


xt = xy[2*l] - xy[2*i];


xy[2*i] = xy[2*i] + xy[2*l];


yt = xy[2*l+1] - xy[2*i+1];


xy[2*i+1] = xy[2*i+1] + xy[2*l+1];



xy[2*l]  
= (c*xt + s*yt)>>15;


xy[2*l+1] = (c*yt - s*xt)>>15;


}


}


ie = ie<<1;


}


}





//生成倒序表子程序


void digitrev_index(short *index, int n,int radix)


{


int  
i, j, k;


short nbits, nbot, ntop, ndiff, n2, raddiv2;


nbits = 0;


i = n;


while (i > 1)


{


i = i >> 1;


nbits++;


}


raddiv2 = radix >> 1;


nbot = nbits >> raddiv2;


nbot = nbot << raddiv2-1;



ndiff  
= nbits & raddiv2;


ntop  
= nbot + ndiff;


n2  

= 1 << ntop;



index[0] = 0;


for(i=1,j=n2/radix+1;i


{


index=j-1;


for (k = n2/radix; k*(radix-1) < j; k /= radix)


j -= k*(radix-1);


j += k;


}


index[n2-1] = n2 - 1;


}







//倒序子程序


void bitrev_cplx(int *x, short *index, int nx)


{


int  

i;


short  


i0, i1, i3;


short  


j0, j1, j3;


int  

xi0, xi1, xi3;


int  

xj0, xj1, xj3;


short  


t;


int  

a, b, ia, ib, ibs;


int  

mask;


int  

nbits, nbot, ntop, ndiff, n2, halfn;


nbits = 0;


i = nx;


while (i > 1)


{


i = i >> 1;


nbits++;


}


nbot  
= nbits >> 1;


ndiff  
= nbits & 1;


ntop  
= nbot + ndiff;


n2  

= 1 << ntop;


mask  
= n2 - 1;


halfn  
= nx >> 1;


for (i0 = 0; i0 < halfn; i0 += 2)


{


b  
= i0 & mask;


a  
= i0 >> nbot;


if (!b) ia  = index[a];


ib  = index;


ibs = ib << nbot;


j0  = ibs + ia;


t  
= i0 < j0;


xi0 = x[i0];


xj0 = x[j0];


if (t){x[i0] = xj0;


x[j0] = xi0;}


i1  = i0 + 1;


j1  = j0 + halfn;


xi1 = x[i1];


xj1 = x[j1];


x[i1] = xj1;


x[j1] = xi1;


i3  = i1 + halfn;


j3  = j1 + 1;


xi3 = x[i3];


xj3 = x[j3];


if (t){x[i3] = xj3;


x[j3] = xi3;}


}


}







欢迎光临 电子技术论坛_中国专业的电子工程师学习交流社区-中电网技术论坛 (http://bbs.eccn.com/) Powered by Discuz! 7.0.0