FFT在DM642中的運行報告
1 FFT編程思想程序清單見附錄）
本算法采用了運算量較大但效果較好的動態溢出檢查。也就是說在計算每一級之前，對所有值進行遍歷，如果發現有超過溢出允許輸入值的變量，則對整個數組進行右移縮小處理。這樣就保證了數值較小的輸入不會被過度縮放，而數值很大的數據也能夠通過FFT運算。由於對所有的數組值進行縮放，因此不會造成失真，但是幅度大小會有變化。因此可以理解為定點算法只求解相對值。下面詳細介紹縮放方法：
一個乘加運算可以理解為：
X1+X2*WN
我們需要找出在什麼情況下輸出會最大，並在設定輸出在此數據類型允許值下時，輸入最大取多少。由於X1，X2，WN都為復數，且abs(X2)=abs(X2*WN)，因此最壞的情況是，在模一定的情況下，原來X2的實部和虛部的數值經過旋轉因子，被全部轉移到了實部，而虛部為零。或者反之。而由於當X2的實部和虛部相等時，它的模值最大，假設全部轉移到實部，在加上X1的實部的值就構成了最大的輸出。我們另這個最大的輸出等於類型最大允許值8bit時為0x7f，16bit時為0x7fff），由於我們需要對所有輸入做統一的限定值以便程序判斷，所以只能有X1.real=X2.real<max，此值得到的最壞情況下輸出為
Sqrt(X2.real^2+X2.imag^2)+X1.real= Sqrt(2*X2.real^2)+X1.real=(sqrt(2)+1)X1.real
我們將此輸出值限定在類型數據允許的最大整數上，則有
(sqrt(2)+1)X1.real<MAX
其中(sqrt(2)+1)為常數=2.4142得到X1.real最大不能超過MAX/2.4142，虛部一樣。繼續觀察運算我們發現，2.4142為大於2小於4的數字，因此，最多可能在同一級上溢出兩個bit，我們動態的調整縮放的量級則會改善上述情況。也就是在溢出兩位時除以4，溢出一位時除以2，經過計算機模擬，發現此時只要滿足單極不溢出的條件如8bit時為52）則總體就不會溢出，但是，付出的代價是兩個bit的精度被處理掉了。因此整體的精度為7bit-2bit8bit-符號位）=5bit，當然誤差可能由於多次這樣的操作而累積。程序中，實現時在每次進入新的一級運算之前，將所有值遍歷，如果超出2倍的ovs值則除以4，超出一倍除以2，否則不處理。這樣，乘加運算在數據類型內部不會造成溢出了。
2 定浮點FFT在DM642中的運行結果運行時間和精度比較）
定浮點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%.
圖1 浮點程序運行結果
圖2 定點程序運行結果
附：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
