Oracle ORION IO 测试工具

[align=center]FFT[/b]在[/b]DM642[/b]中的运行报告[/b][/b][/align]

1 FFT[/b]编程思想（程序清单见附录）[/b][/b]
本算法采用了运算量较大但效果较好的动态溢出检查。也就是说在计算每一级之前，对所有值进行遍历，如果发现有超过溢出允许输入值的变量，则对整个数组进行右移缩小处理。这样就保证了数值较小的输入不会被过度缩放，而数值很大的数据也能够通过FFT运算。由于对所有的数组值进行缩放，因此不会造成失真，但是幅度大小会有变化。因此可以理解为定点算法只求解相对值。下面详细介绍缩放方法：
一个乘加运算可以理解为：
[align=center]X1+X2*WN[/align]
我们需要找出在什么情况下输出会最大，并在设定输出在此数据类型允许值下时，输入最大取多少。由于X1，X2，WN都为复数，且abs(X2)=abs(X2*WN)，因此最坏的情况是，在模一定的情况下，原来X2的实部和虚部的数值经过旋转因子，被全部转移到了实部，而虚部为零。或者反之。而由于当X2的实部和虚部相等时，它的模值最大，假设全部转移到实部，在加上X1的实部的值就构成了最大的输出。我们另这个最大的输出等于类型最大允许值（8bit时为0x7f，16bit时为0x7fff），由于我们需要对所有输入做统一的限定值以便程序判断，所以只能有X1.real=X2.real<max，此值得到的最坏情况下输出为
[align=center]Sqrt(X2.real^2+X2.imag^2)+X1.real= Sqrt(2*X2.real^2)+X1.real=(sqrt(2)+1)X1.real[/align]
我们将此输出值限定在类型数据允许的最大整数上，则有
[align=center](sqrt(2)+1)X1.real<MAX[/align]
其中(sqrt(2)+1)为常数=2.4142得到X1.real最大不能超过MAX/2.4142，虚部一样。继续观察运算我们发现，2.4142为大于2小于4的数字，因此，最多可能在同一级上溢出两个bit，我们动态的调整缩放的量级则会改善上述情况。也就是在溢出两位时除以4，溢出一位时除以2，经过计算机模拟，发现此时只要满足单极不溢出的条件（如8bit时为52）则总体就不会溢出，但是，付出的代价是两个bit的精度被处理掉了。因此整体的精度为7bit-2bit（8bit-符号位）=5bit，当然误差可能由于多次这样的操作而累积。程序中，实现时在每次进入新的一级运算之前，将所有值遍历，如果超出2倍的ovs值则除以4，超出一倍除以2，否则不处理。这样，乘加运算在数据类型内部不会造成溢出了。
[/b]
2 [/b]定浮点[/b]FFT[/b]在[/b]DM642[/b]中的运行结果（运行时间和精度比较）[/b][/b]
定浮点FFT程序在CCS中的运行结果如下图1，2。
经过在MATLAB中验证数据浮点程序运行结果正确。对于定点程序，我们看到标识位三个ov=1；故输入的八个数据总共被缩小2^3=8倍。所以，把所有的定点数据换算成浮点数据然后再乘以8即得到实际数据，经验算结果正确。
本例在CCS中设定的周期时间是40ns，观察两个程序的运行周期，浮点程序的运行周期是99510个，所以其运行时间是99510*40ns=3.9804ms。定点程序的运行周期是56143个，所以其运行时间是56143*40ns=2.24572ms。二者相差的时间数量级是1.7个ms，相差的相对百分比是（3.9804-2.24572）/9.9804=43.58%。
对于定点运算的精度，计算公式如下：n=m/32768*8.最后得到的16个数据为：
3.598877 0 -0.400391 0.965576 -0.399902 0.399902 -0.400391 0.165771
0.399658 0 -0.399902 0.165771 -0.399902 -0.399902 -0.399902 0.965576
与浮点所得结果差值的绝对值为：
0.001123 0 0.000391 0.000104 0.000098 0.000098 0.000391 0.000089
0.000342 0 0.000098 0.000089 0.000098 0.000098 0.000098 0.000104
其平均值为：0.000226 其差值的百分比为0.003612/9.06272=0.03986%.



[align=center]图1 浮点程序运行结果[/align]



[align=center]图2 定点程序运行结果[/align]

附：FFT定浮点程序清单
浮点程序：
#include<stdio.h>
struct complex
{
double real;
double imag;
};
struct complex x[8]={{0.1,0.0},{0.2,0.0},{0.3,0.0},{0.4,0.0},{0.5,0.0},{0.6,0.0},{0.7,0.0},{0.8,0.0}};
struct complex cheng,y[8];
double costable[7]={1.0000,1.0000,0.0000,1.0000,0.7071,0.0000,-0.7071},sintable[7]={0.0000,0.0000,1.0000,0.0000,0.7071,1.000,0.7071};
int reverse(int i)
{
int a[3],b,j=0;
for(b=0;b<3;b++)
{
a[b]=(1&i)<<(2-b);
j=j+a[b];
i/=2;
}
return j;
}
void paixu()
{
int i,j;
for(i=0;i<8;i++)
{
j=reverse(i);
y[j].real=x[i].real;
y[j].imag=x[i].imag;
}
}
void main()
{

int m=4,n=1,i,j,k,num1,num2;
paixu();
for(i=0;i<3;i++)
{
for(j=0;j<m;j++)
{
for(k=0;k<n;k++)
{
num1=2*n*j+k;
num2=num1+n;
cheng.real=y[num2].real*costable[(1<<i)+k-1]+y[num2].imag*sintable[(1<<i)+k-1];
cheng.imag=y[num2].imag*costable[(1<<i)+k-1]-y[num2].real*sintable[(1<<i)+k-1];
y[num2].real=y[num1].real-cheng.real;
y[num2].imag=y[num1].imag-cheng.imag;
y[num1].real=y[num1].real+cheng.real;
y[num1].imag=y[num1].imag+cheng.imag;
}
}
m/=2;
n*=2;
}
for(i=0;i<8;i++)
{
printf("%f\n%f\n",y[i].real,y[i].imag);
}

}

定点程序：
#include<stdio.h>
struct complex
{
int real;
int imag;
};
struct complex x[8]={{3276,0},{6553,0},{9830,0},{13107,0},{16384,0},{19660,0},{22937,0},{26214,0}};
struct complex cheng,y[8];
int costable[7]={32767,32767,0,32767,23170,0,-23170},sintable[7]={0,0,32767,0,23170,32767,23170};
int reverse(int i)
{
int a[3],b,j=0;
for(b=0;b<3;b++)
{
a[b]=(1&i)<<(2-b);
j=j+a[b];
i/=2;
}
return j;
}
void paixu()
{
int i,j;
for(i=0;i<8;i++)
{
j=reverse(i);
y[j].real=x[i].real;
y[j].imag=x[i].imag;
}
}
void main()
{
int p,i,j,k,num1,num2,m=4,n=1,ov=0,ovd=0;
paixu();
for(i=0;i<3;i++)
{

for(p=0;p<8;p++)
if (y[p].real>27146 || y[p].imag>27146 || y[p].real<-27146 || y[p].imag<-27146) //超出两倍允许值，所有值右移2个bit
{
ovd=1;
break;
}
else if (y[p].real>13573 || y[p].imag>13573 || y[p].real<-13573 || y[p].imag<-13573) //超出允许值，所有值右移1个bit
{
ov=1;
}
if (ovd==1)
for(p=0;p<8;p++)
{
y[p].real=y[p].real>>2;
y[p].imag=y[p].imag>>2;
}
else if(ov==1)
for(p=0;p<8;p++)
{
y[p].real=y[p].real>>1;
y[p].imag=y[p].imag>>1;
}
printf("%d\n%d\n",ovd,ov);
for(j=0;j<m;j++)
{
for(k=0;k<n;k++)
{

num1=2*n*j+k;
num2=num1+n;
cheng.real=((y[num2].real*costable[(1<<i)+k-1])>>15)+((y[num2].imag*sintable[(1<<i)+k-1])>>15);
cheng.imag=((y[num2].imag*costable[(1<<i)+k-1])>>15)-((y[num2].real*sintable[(1<<i)+k-1])>>15);
//cheng.real=x[num2].real*cos(pi*k*m/4)+x[num2].imag*sin(pi*k*m/4);
//cheng.imag=x[num2].imag*cos(pi*k*m/4)-x[num2].real*sin(pi*k*m/4);
y[num2].real=y[num1].real-cheng.real;
y[num2].imag=y[num1].imag-cheng.imag;
y[num1].real=y[num1].real+cheng.real;
y[num1].imag=y[num1].imag+cheng.imag;
}
}
m>>=1;
n<<=1;
}
for(i=0;i<8;i++)
{
printf("%d\n%d\n",y[i].real,y[i].imag);
}

}
本文出自 “嵌入式学习” 博客，请务必保留此出处http://mousimin.blog.51cto.com/1624003/319026