经过一段时间的学习，今天在CCS内成功仿真了基2频域抽取FFT(radix2 DIF FFT)。
环境：CCS 3.1 + C6201 Device Simulator
仿真结果如下图：
![]() ![]()
全部代码在下面，为了描述简单，把所有函数写在了一起，实际运用中不应该这样做。另外位反转索引表和旋转因子表都是在里面实时生成的，实际运用时，应该做好表格后嵌入，或在编译时生成而不是运行时生成。
感谢 pacificxu 的教程，他貌似曾经是闻亭的高手。
#include “math.h”
/*========================================*/
/* GLOBAL CONSTANTS */
/*========================================*/
#define PI 3.14159265358979
#define N 2048
/*========================================*/
/* GLOBAL VARIABLES */
/*========================================*/
static short index1[64], w1[N];
static short x1[N*2];
static int x2[N];
int i;
void DSP_radix2(int n, short *restrict xy, const short *restrict w)
{
short n1,n2,ie,ia,i,j,k,l;
short xt,yt,c,s;
n2 = n;
ie = 1;
for (k=n; 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 < n; 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 DSP_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;}
}
}
void bitrev_index(short *index, int n)
{
int i, j, k, radix = 2;
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 < n2 – 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;
}
/*===================================*/
/* MAIN() */
/*===================================*/
void main()
{
double delta;
/*======== init the index table for FFT bitrev =====*/
bitrev_index(index1, N);
/*======== init the FFT coeficients =====*/
delta = 2*PI/N;
for (i=0;i<N/2;i++)
{
w1[2*i] = 32767*(-cos((double)i*delta));
w1[2*i+1] = 32767*(-sin((double)i*delta));
}
/*======== init the input data =====*/
for (i=0;i<N;i++)
{
x1[2*i] = (short)((cos(PI*i/20.0)+cos(PI*i/10.0)+cos(PI*i/5.0))*0×80);
x1[2*i+1] = 0;
}
/*======== implement FFT =====*/
DSP_radix2(N, x1, w1);
for (i=0; i<N; i++)
{
x2 = sqrt(x1[2*i]*x1[2*i]+x1[2*i+1]*x1[2*i+1]); //—?
}
/*======== bitrev =====*/
DSP_bitrev_cplx(x2, index1, N); |
楼主
标题置顶
标题高亮
