1. 程式人生 > >DSP中的定點運算

DSP中的定點運算

DSP中的定點算數運算
1 數的定標
    在定點DSP晶片中,採用定點數進行數值運算,其運算元一般採用整型數來表示。一個整型數的最大表示範圍取決於DSP晶片所給定的字長,一般為16位或24位。顯然,字長越長,所能表示的數的範圍越大,精度也越高。如無特別說明,本書均以16位字長為例。
DSP晶片的數以2的補碼形式表示。每個16位數用一個符號位來表示數的正負,0表示數值為正,l則表示數值為負。其餘15位表示數值的大小。因此,
      二進位制數0010000000000011b=8195
      二進位制數1111111111111100b= -4
    對DSP晶片而言,參與數值運算的數就是16位的整型數。但在許多情況下,數學運算過程中的數不一定都是整數。那麼,DSP晶片是如何處理小數的呢?應該說,DSP晶片本身無能為力。那麼是不是說DSP晶片就不能處理各種小數呢?當然不是。這其中的關鍵就是由程式設計師來確定一個數的小數點處於16位中的哪一位。這就是數的定標。
通過設定小數點在16位數中的不同位置,就可以表示不同大小和不同精度的小數了。數的定標有Q表示法和S表示法兩種。表1.1列出了一個16位數的16種Q表示、S表示及它們所能表示的十進位制數值範圍。
    從表1.1可以看出,同樣一個16位數,若小數點設定的位置不同,它所表示的數也就不同。例如,
         16進位制數2000H=8192,用Q0表示
         16進位制數2000H=0.25,用Q15表示
但對於DSP晶片來說,處理方法是完全相同的。
    從表1.1還可以看出,不同的Q所表示的數不僅範圍不同,而且精度也不相同。Q越大,數值範圍越小,但精度越高;相反,Q越小,數值範圍越大,但精度就越低。例如,Q0 的數值範圍是一32768到+32767,其精度為1,而Q15的數值範圍為-1到0.9999695,精度為1/32768=0.00003051。因此,對定點數而言,數值範圍與精度是一對矛盾,一個變數要想能夠表示比較大的數值範圍,必須以犧牲精度為代價;而想精度提高,則數的表示範圍就相應地減小。在實際的定點演算法中,為了達到最佳的效能,必須充分考慮到這一點。
浮點數與定點數的轉換關係可表示為:
        浮點數(x)轉換為定點數(xq):xq=(int)x* 2Q
        定點數(xq)轉換為浮點數(x):x=(float)xq*2-Q
    例如,浮點數x=0.5,定標Q=15,則定點數xq=L0.5*32768J=16384,式中LJ表示下取整。反之,一個用Q=15表示的定點數16384,其浮點數為163幼*2-15=16384/32768=0.5。浮點數轉換為定點數時,為了降低截尾誤差,在取整前可以先加上0.5。

表1.1    Q表示、S表示及數值範圍
Q表示    S表示    十進位制數表示範圍
Q15    S0.15    -1≤x≤0.9999695
Q14    S1.14    -2≤x≤1.9999390
Q13    S2.13    -4≤x≤3.9998779
Q12    S3.12    -8≤x≤7.9997559
Q11    S4.11    -16≤x≤15.9995117
Q10    S5.10    -32≤x≤31.9990234
Q9     S6.9     -64≤x≤63.9980469
Q8     S7.8     -128≤x≤127.9960938
Q7     S8.7     -256≤x≤255.9921875
Q6     S9.6     -512≤x≤511.9804375
Q5     S10.5    -1024≤x≤1023.96875
Q4     S11.4    -2048≤x≤2047.9375
Q3     S12.3    -4096≤x≤4095.875
Q2     S13.2    -8192≤x≤8191.75
Q1     S14.1    -16384≤x≤16383.5
Q0     S15.0    -32768≤x≤32767

2 高階語言:從浮點到定點
    我們在編寫DSP模擬演算法時,為了方便,一般都是採用高階語言(如C語言)來編寫模擬程式。程式中所用的變數一般既有整型數,又有浮點數。如例1.1程式中的變數i是整型數,而pi是浮點數,hamwindow則是浮點陣列。
例1.1 256點漢明窗計算
int i;+
float pi=3.14l59;
float hamwindow[256];
for(i=0;i<256;i++) hamwindow=0.54-0.46*cos(2.0*pi*i/255);
    如果我們要將上述程式用某種足點DSP晶片來實現,則需將上述程式改寫為DSP晶片的組合語言程式。為了DSP程式除錯的方便及模擬定點DSP實現時的演算法效能,在編寫DSP彙編程式之前一般需將高階語言浮點演算法改寫為高階語言定點演算法。下面我們討論基本算術運算的定點實現方法。
2.1 加法/減法運算的C語言定點摸擬
設浮點加法運算的表示式為:
float x,y,z;
z=x+y;
將浮點加法/減法轉化為定點加法/減法時最重要的一點就是必須保證兩個運算元的定標
temp=x+temp; 
z=temp>>(Qx-Qz),若Qx>=Qz
z=temp<<(Qz-Qx),若Qx<=Qz
例1.4結果超過16位的定點加法
設x=l5000,y=20000,則浮點運算值為z=x+y=35000,顯然z>32767,因此
Qx=1,Qy=0,Qz=0,則定點加法為:
x=30000;y=20000;
temp=20000<<1=40000; 
temp=temp+x=40000+30000=70000;
z=70000L>>1=35000;
    因為z的Q值為0,所以定點值z=35000就是浮點值,這裡z是一個長整型數。當加法或加法的結果超過16位表示範圍時,如果程式設計師事先能夠了解到這種情況,並且需要保持運算精度時,則必須保持32位結果。如果程式中是按照16位數進行運算的,則超過16位實際上就是出現了溢位。如果不採取適當的措施,則資料溢位會導致運算精度的嚴重惡化。一般的定點DSP晶片都沒有溢位保護功能,當溢位保護功能有效時,一旦出現溢位,則累加器ACC的結果為最大的飽和值(上溢為7FFFH,下溢為8001H),從而達到防止溢位引起精度嚴重惡化的目的。
2.2乘法運算的C語言定點模擬
設浮點乘法運算的表示式為:
float x,y,z;
z=xy; 
假設經過統計後x的定標值為Qx,y的定標值為Qy,乘積z的定標值為Qz,則
z=xy
zq*2-Qx=xq*yq*2-(Qx+Qy)
zq=(xqyq)2Qz-(Qx+Qy)
所以定點表示的乘法為:
int x,y,z;
long temp;
temp=(long)x; 
z=(temp*y)>>(Qx+Qy-Qz);
例1.5定點乘法。
設x=18.4,y=36.8,則浮點運算值為=18.4*36.8=677.12;
根據上節,得Qx=10,Qy=9,Qz=5,所以
x=18841;y=18841;
temp=18841L;
z=(18841L*18841)>>(10+9-5)=354983281L>>14=21666;
因為z的定標值為5,故定點z=21666,即為浮點的z=21666/32=677.08。
2.3除法運算的C語言定點摸擬
設浮點除法運算的表示式為:
float x,y,z;
z=x/y;
假設經過統計後被除數x的定標值為Qx,除數y的定標值為Qy,商z的定標值為Qz,則
z=x/y
zq*2-Qz=(xq*2-Qx)/(yq*2-Qy)
zq=(xq*2(Qz-Qx+Qy))/yq
所以定點表示的除法為:
int x,y,z;
long temp;
temp=(long)x;
z=(temp<<(Qz-Qx+Qy))/y;
例1.6定點除法。
設x=18.4,y=36.8,浮點運算值為z=x/y=18.4/36.8=0.5;
根據上節,得Qx=10,Qy=9,Qz=15;所以有
z=18841,y=18841;
temp=(long)18841;
z=(18841L<<(15-10+9)/18841=3O8690944L/18841=16384;
因為商z的定標值為15,所以定點z=16384,即為浮點z=16384/215=0.5。
2.4程式變數的Q值確定
    在前面幾節介紹的例子中,由於x,y,z的值都是已知的,因此從浮點變為定點時Q值很好確定。在實際的DSP應用中,程式中參與運算的都是變數,那麼如何確定浮點程式中變數的Q值呢?從前面的分析可以知道,確定變數的Q值實際上就是確定變數的動態範圍,動態範圍確定了,則Q值也就確定了。
設變數的絕對值的最大值為|max|,注意|max|必須小於或等於32767。取一個整數n,使滿足
2n-1<|max|<2n
則有
2-Q=2-15*2n=2-(15-n)
Q=15-n
例如,某變數的值在-1至+1之間,即|max|<1,因此n=0,Q=15-n=15。
    既然確定了變數的|max|就可以確定其Q值,那麼變數的|max|又是如何確定的呢?一般來說,確定變數的|max|有兩種方法。一種是理論分析法,另一種是統計分析法。
1. 理論分析法
    有些變數的動態範圍通過理論分析是可以確定的。例如:
(1)三角函式。y=sin(x)或y=cos(x),由三角函式知識可知,|y|<=1。
(2)漢明窗。y(n)=0.54一0.46cos[nπn/(N-1)],0<=n<=N-1。因為-1<=cos[2πn/(N-1)]<=1,所以0.08<=y(n)<=1.0。
(3)FIR卷積。y(n)=∑h(k)x(n-k),設∑|h(k)|=1.0,且x(n)是模擬訊號12位量化值,即有|x(n)|<=211,則|y(n)|<=211。
(4)理論已經證明,在自相關線性預測編碼(LPC)的程式設計中,反射係數ki滿足下列不等式:|ki|<1.0,i=1,2,...,p,p為LPC的階數。
2. 統計分析法
    對於理論上無法確定範圍的變數,一般採用統計分析的方法來確定其動態範圍。所謂統計分析,就是用足夠多的輸入訊號樣值來確定程式中變數的動態範圍,這裡輸入訊號一方面要有一定的數量,另一方面必須儘可能地涉及各種情況。例如,在語音訊號分析中,統計分析時就必須來集足夠多的語音訊號樣值,並且在所採集的語音樣值中,應儘可能地包含各種情況。如音量的大小,聲音的種類(男聲、女聲等)。只有這樣,統計出來的結果才能具有典型性。
    當然,統計分析畢竟不可能涉及所有可能發生的情況,因此,對統計得出的結果在程式設計時可採取一些保護措施,如適當犧牲一些精度,Q值取比統計值稍大些,使用DSP晶片提供的溢位保護功能等。
2.5浮點至定點變換的C程式舉例
    本節我們通過一個例子來說明C程式從浮點變換至定點的方法。這是一個對語音訊號(0.3~3.4kHz)進行低通濾波的C語言程式,低通濾波的截止頻率為800Hz,濾波器採用19點的有限衝擊響應FIR濾波。語音訊號的取樣頻率為8kHz,每個語音樣值按16位整型數存放在insp.dat檔案中。
例1.7語音訊號800Hz 19點FIR低通濾波C語言浮點程式。
#i nclude <stdio.h>
const int length=180/*語音幀長為180點=22.5ms@8kHz取樣*/
void filter(int xin[],int xout[],int n,float h[]);/*濾波子程式說明*/
/*19點濾波器係數*/
static float h[19]= 
{0.01218354,-0.009012882,-0.02881839,-0.04743239,-0.04584568,
-0.008692503,0.06446265,0.1544655,0.2289794,0.257883,
0.2289794,0.1544655,0.06446265,-0.008692503,-0.04584568,
-0.04743239,-0.02881839,-0.009012882,O.01218354};
static int xl[length+20];
/*低通濾波浮點子程式*/
void filter(int xin[],int xout[],int n,float h[])
{
int i,j;
float sum;
for(i=0;i<length;i++)x1[n+i-1]=xin;
for(i=0;i<length;i++)
{
sum=0.0;
for(j=0;j<n;j++)sum+=h[j]*x1[i-j+n-1];
xout=(int)sum;
for(i=0;i<(n-l);i++)x1[n-i-2]=xin[length-1-i];

/*主程式*/
void main()
FILE *fp1,*fp2;
int ,indata[length],outdata[length];
fp1=fopen(insp.dat,"rb");/* 輸入語音檔案*/
fp2=fopen(Outsp.dat,"wb");/* 濾波後語音檔案*/
=0;
while(feof(fp1) ==0) 
{
++;
printf(“=%d\n”,);
for(i=0;i<length;i++)indata=getw(fp1); /*取一幀語音資料*/
filter(indata,outdata,19,h);/*呼叫低通濾波子程式*/
for(i=0;i<length;i++)putw(outdata,fp2);/*將濾波後的樣值寫入檔案*/

fcloseall();/*關閉檔案*/
return(0);
}
例1.8語音訊號800Hz l9點FIR低通濾波C語言定點程式。
#i nclude <stdio.h>
const int length=180; 
void filter (int xin[],int xout[],int n,int h[]);
static int h[19]={399,-296,-945,-1555,-1503,-285,2112,5061,7503,8450,
7503,5061,2112,-285,-1503,-1555,-945,-296,399};/*Q15*/
static int x1[length+20];
/*低通濾波定點子程式*/
void filter(int xin[],int xout[],int n,int h[])
int i,j;
long sum;
for(i=0;i<length;i++)x1[n+i-111=xin];
for(i=0;i<1ength;i++)
sum=0;
for(j=0;j<n;j++)sum+=(long)h[j]*x1[i-j+n-1];
xout=sum>>15;
for(i=0;i<(n-1);i++)x1[n-i-2]=xin[length-i-1];
}
主程式與浮點的完全一樣。“
3 DSP定點算術運算
    定點DSP晶片的數值表示基於2的補碼錶示形式。每個16位數用l個符號位、i個整數位和15-i個小數位來表示。因此:
00000010.10100000
表示的值為:
21+2-1+2-3=2.625
     這個數可用Q8格式(8個小數位)來表示,其表示的數值範圍為-128至+l27.996,一個Q8定點數的小數精度為1/256=0.004。
    雖然特殊情況(如動態範圍和精度要求)必須使用混合表示法。但是,更通常的是全部以Q15格式表示的小數或以Q0格式表示的整數來工作。這一點對於主要是乘法和累加的訊號處理演算法特別現實,小數乘以小數得小數,整數乘以整數得整數。當然,乘積累加時可能會出現溢位現象,在這種情況下,程式設計師應當瞭解數學裡面的物理過程以注意可能的溢位情況。下面我們來討論乘法、加法和除法的DSP定點運算,彙編程式以TMS320C25為例。
3.1定點乘法
    兩個定點數相乘時可以分為下列三種情況:
1. 小數乘小數
例1.9 Q15*Q15=Q30
0.5*0.5=0.25
0.100000000000000;Q15
* 0.100000000000000;Q15
--------------------------------------------
00.010000000000000000000000000000=0.25;Q30
    兩個Q15的小數相乘後得到一個Q30的小數,即有兩個符號位。一般情況下相乘後得到的滿精度數不必全部保留,而只需保留16位單精度數。由於相乘後得到的高16位不滿15位的小資料度,為了達到15位精度,可將乘積左移一位,下面是上述乘法的TMS320C25程式:
LT OP1;OP1=4000H(0.5/Q15)
MPY OP2;oP2=4000H(0.5/Ql5)
PAC
SACH ANS,1;ANS=2000H(0.25/Q15)
2. 整數乘整數
例1.10 Q0*Q0=Q0
17*(-5)=-85
0000000000010001=l7
*1111111111111011=-5
-------------------------------------------
11111111111111111111111110101011=-85
3. 混合表示法
    許多情況下,運算過程中為了既滿足數值的動態範圍又保證一定的精度,就必須採用Q0與Q15之間的表示法。比如,數值1.2345,顯然Q15無法表示,而若用Q0表示,則最接近的數是1,精度無法保證。因此,數1.2345最佳的表示法是Q14。
例1.11 1.5*0.75= 1.125
01.10000000000000=1.5;Q14
*00.11000000000000=0.75;Q14
---------------------------------------
0001.0010000000000000000000000000=1.125 Q28
    Q14的最大值不大於2,因此,兩個Q14數相乘得到的乘積不大於4。
一般地,若一個數的整數位為i位,小數位為j位,另一個數的整數位為m位,小數位為n位,則這兩個數的乘積為(i+m)位整數位和(j+n)位小數位。這個乘積的最高16位可能的精度為(i+m)整數位和(15- i- m)小數位。
    但是,若事先了解數的動態範圍,就可以增加數的精度。例如,程式設計師瞭解到上述乘積不會大於1.8,就可以用Q14數表示乘積,而不是理論上的最佳情況Q13。例3.11的TMS320C25程式如下:
LT OP1;OP1 = 6000H(1.5/Ql4)
MPY OP2;OP2 = 3000H(0.75/Q14)
PAC
SACH ANS,1;ANS=2400H(1.125/Q13)
    上述方法,為了精度均對乘的結果舍位,結果所產生的誤差相當於減去一個LSB(最低位)。採用下面簡單的舍人方法,可使誤差減少二分之一。
LT OP1
MPY OP2
PAC
ADD ONE,14(上舍入)
SACH ANS,1
    上述程式說明,不管ANS為正或負,所產生的誤差是l/2 LSB,其中儲存單元ONE的值為1。
3.2定點加法
    乘的過程中,程式設計師可不考慮溢位而只需調整運算中的小數點。而加法則是一個更加複雜的過程。首先,加法運算必須用相同的Q點表示,其次,程式設計師或者允許其結果有足夠的高位以適應位的增長,或者必須準備解決溢位問題。如果運算元僅為16位長,其結果可用雙精度數表示。下面舉例說明16位數相加的兩種途徑。
1.保留32位結果
LAC OP1;(Q15)
ADD OP2;(Ql5)
SACH ANSHI ;(高16位結果)
SACL ANSLO :(低16位結果)
2.調整小數點保留16位結果
LAC OP1,15;(Q14數用ACCH表示)
ADD OP2,15;(Q14數用ACCH表示)
SACH ANS;(Q14)
    加法運算最可能出現的問題是運算結果溢位。TMS320提供了檢查溢位的專用指令BV,此外,使用溢位保護功能可使累加結果溢位時累加器飽和為最大的整數或負數。當然,即使如此,運算精度還是大大降低。因此,最好的方法是完全理解基本的物理過程並注意選擇數的表達方式。
3.3定點除法
    在通用DSP晶片中,一般不提供單週期的除法指令,為此必須採用除法子程式來實現。二進位制除法是乘法的逆運算。乘法包括一系列的移位和加法,而除法可分解為一系列的減法和移位。下面我們來說明除法的實現過程。
    設累加器為8位,且除法運算為10除以3。除的過程包括與被除法有關的除數逐步移位,在每一步進行減法運算,如果能減則將位插入商中。
(1)除數的最低有效位對齊被除數的最高有效位。
0000l0l0
- 00011000
--------------------------------------
11110010
(2)由於減法結果為負,放棄減法結果,將被除數左移一位,再減。
00010100
- 00011000
----------------------------------------
11111000
(3)結果仍為負,放棄減法結果,被除數左移一位,再減。
00101000
-   00011000
------------------------------------------
00010000
(4)結果為正,將減法結果左移一位後加1,作最後一次減。
00100001
-   00011000
----------------------------------------
00001001
(5)結果為正,將結果左移一位加1 得最後結果。高4位代表餘數,低4位表示商。
00010011
即,商為0011= 3.餘數為0001= 1。
    TMS320沒有專門的除法指令,但使用條件減指令SUBC可以完成有效靈活的除法功能。使用這一指令的唯一限制是兩個運算元必須為正。程式設計師必須事先了解其可能的運算數的特性,如其商是否可以用小數表示及商的精度是否可被計算出來。這裡每一種考慮可影響如何使用SUBC指令的問題。下面我們給出兩種不同情況下的TMS320C25除法程式。
(1)分子小於分母
DIV_A:
LT NUMERA
MPY DENOM
PAC
SACH TEMSGN;取商的符號
LAC DENOM
ABS
SACL DENOM;使分母為正
ZALH NUMERA; 分子為正
ABS
RPTK 14
SUBC DENOM;除迴圈15次
SACL QUOT
LAC TEMSGN
BGEZ A1;若符號為正,則完成
ZAC
SUB QUOT
SACL QUOT;若為負,則商為負
A1: RET
    這個程式中,分子在NUMERA中,分母在DENOM中,商存在QUOT中,TEMSGN為暫存單元。
(2)規定商的精度
DIV_B: 
LT NUMERA
MPY DENOM
PAC 
SACH TEMSGN;取商的符號
LAC DENOM
ABS
SACL DENOM; 使分母為正
LACK 15
ADD FRAC
SACL FRAC;計算迴圈計數器
LAC NUMERA
ABS ; 使分子為正
RPT FRAC
SUBC DENOM; 除迴圈16+FRAC次
SACL QUOT
LAC TEMSGN
BGEZ B1;若符號為正,則完成
ZAC
SUB QUOT
SACL QUOT;若為負,則商為負
B1: RET
    與DIV_A相同,這個程式中,分子在NUMERA中,分母在DENOM中,商存在QUOT中,TEMSGN為暫存單元。FRAC中規定商的精度,如商的精度為Q13,則呼叫程式前FRAC單元中的值應為13。
4 非線性運算的定點快速實現
    在數值運算中,除基本的加減乘除運算外,還有其它許多非線性運算,如,對數運算,開方運算,指數運算,三角函式運算等,實現這些非線性運算的方法一般有:(1)呼叫DSP編譯系統的庫函式;(2)查表法;(3)混合法。下面我們分別介紹這三種方法。
1.呼叫DSP編譯系統的庫函式
    TMS320C2X/C5X的C編譯器提供了比較豐富的執行支援庫函式。在這些庫函式中,包含了諸如對數、開方、三角函式、指數等常用的非線性函式。在C程式中(也可在彙編程式中)只要採用與庫函式相同的變數定義,就可以直接呼叫。例如,在庫函式中,定義了以10為底的常用對數log10():
#i nclude<math.h>
double,log10(double x);
在C程式中按如下方式呼叫:
float x,y;
X=10.0;
y=log10(x);
    從上例可以看出,庫函式中的常用對數log10()要求的輸入值為浮點數,返回值也為浮點數,運算的精度完全可以保證。直接呼叫庫函式非常方便,但由於運算量大,很難在實時DSP中得到應用。
2.查表法
    在實時DSP應用中實現非線性運算,一般都採取適當降低運算精度來提高程式的運算速度。查表法是快速實現非線性運算最常用的方法。採用這種方法必須根據自變數的範圍和精度要求製作一張表格。顯然輸人的範圍越大,精度要求越高,則所需的表格就越大,即儲存量也越大。查表法求值所需的計算就是根據輸入值確定表的地址,根據地址就可得到相應的值,因而運算量較小。查表法比較適合於非線性函式是周期函式或已知非線性函式輸入值範圍這兩種情況、例1.12和例1. 13分別說明這兩種情況。
例1.12 已知正弦函式y=cos(x),製作一個512點表格,並說明查表方法。由於正弦函式是周期函式,函式值在-1至+1之間,用查表法比較合適。由於Q15的表示範圍為1-至32767/32768之間,原則上講-1至+1的範圍必須用Q14表示。但一般從方便和總體精度考慮,類似情況仍用Q15表示,此時+1用32767來表示。
(1)產生5l2點值的C語言程式如下所示。
#define N 512
#define pi 3.14l59
int sin_tab[5l2];
void main()
{
int i;
for(i=0;i<N;i++)sin_tab=(int)(32767*sin(2*pi*i/N));
(2)查表
    查表實際上就是根據輸人值確定表的地址。設輸入x在0~2π之間,則x對應於512點表的地址為:index=(int)(512*x/2π),則y=sin(x)=sin_tab[index]如果x用Q12定點數表示,將512/2π用Q8表示為20861,則計算正弦表的地址的公式為。
index=(x*20861L)>>20;
例1.12用查表法求以2為底的對數,已知自變數值範圍為0.5-1,要求將自變數範圍均勻劃分為10等分。試製作這個表格並說明查表方法。
(1)作表:
    y=log2(x),由於x在0.5到1之間,因此y在-1到0之間,x和y均可用Q15表示。由於對x均勻劃分為10段,因此,10段對應於輸入x的範圍如表3.2所示。若每一段的對數值都取第一點的對數值,則表中第一段的對數值為y0(Q15)=(int)(log(O.5)*32768),第二段的對數值為y1(Q15)=(int)(log2(0.55)*32768),依次類推,如表3.2所示。
(2)查表:
    查表時,先根據輸人值計算表的地址,計算方法為:
index=((x-16384)*20)>>15;
式中, index就是查表用的地址。例如,已知輸人x=26869,則index=6,因此,y= -10549。

表1.2 logtab0 10點對數表
地址    輸入值    對數值(Q15)
0    0.50-0.55    -32768
1    0.55-0.60    -28262
2    0.60-0.65    -24149
3    0.65-0.70    -20365
4    0.70-0.75    -16862
5    0.75-0.80    -13600
6    0.80-0.85    -10549
7    0.85-0.90    -7683
8    0.90-0.95    -4981
9    0.95-1.00    -2425

3.混合法
(1)提高查表法的精度
    上述方法查表所得結果的精度隨表的大小而變化,表越大,則精度越高,但儲存量也越大。當系統的儲存量有限而精度要求也較高時,查表法就不太適合。那麼能否在適當增加運算量的情況下提高非線性運算的精度呢?下面介紹一種查表結合少量運算來計算非線性函式的混合法,這種方法適用於在輸入變數的範圍內函式呈單調變化的情形。混合法是在查表的基礎上來用計算的方法以提高當輸入值處於表格兩點之間時的精度。提高精度的一個簡便方法是採用折線近似法,如圖1.1所示。
圖1.1提高精度的折線近似法”
    仍以求以2為底的對數為例(例1.12)。設輸入值為x,則精確的對數值為y,在表格值的兩點之間作一直線,用y'作為y的近似值,則有:
      y'=y0+△y
其中y0由查表求得。現在只需在查表求得y0的基礎上增加△y既可。△y的計算方法如下:                   △y=(△x/△x0)△y=△x(△y0/△x0)
    其中△y0/△x0對每一段來說是一個恆定值,可作一個表格直接查得。此外計算此時需用到每段橫座標的起始值,這個值也可作一個表格。這佯共有三個大小均為10的表格,分別為儲存每段起點對數值的表logtab0、儲存每段△y0/△x0值的表logtab1和儲存每段輸入起始值x0的表logtab2,表logtab1和表logtab2可用下列兩個陣列表示。
int logtab1[10]={22529,20567,18920,17517,16308,
15255,14330,13511,12780,12124};/*△y0/△x0:Q13*/
int logtab2[10]={16384,18022,19660,21299,22938,
24576,26214,27853,29491,31130};/*x0:Q15*/
綜上所述,採用混合法計算對數值的方法可歸納為:
(1)根據輸人值,計算查表地址:index=((x-16384)*20)>>15;
(2)查表得y0=logtab0[index];
(3)計算△x=x-logtab2[index];
(4)計算△y=(△x*logtab1[index])>>13;
(5)計算得結果y=y0+△y。
例1.13已知x=0.54,求log2(x)。
0.54的精確對數值為y=log2(0.54)=-0.889。
    混合法求對數值的過程為:
(1)定標Q15,定標值x=0.54*32768=17694;
(2)表地址index=((x-16384)*20)>>15=0;
(3)查表得y0=logtab0[0]=-32768;
(4)計算△x=x-logtab2[0]=17694-16384=1310;
(5)計算△y=(△xlogtab1[0]>>13=(13l0*22529L)>>13=3602 
(6)計算結果y=y0+△y=-32768+3602=-29166。
    結果y為Q15定標,析算成浮點數為-29166/32768=-0.89,可見精度較高。
(2)擴大自變數範圍
    如上所述,查表法比較適用於周期函式或自變數的動態範圍不是太大的情形。對於像對數這樣的非線性函式,輸入值和函式值的變化範圍都很大。如果輸入值的變化範圍很大,則作表就比較困難。那麼能否比較好地解決這個問題,即不便表格太大,又能得到比較高的精度呢?下面我們來討論一種切實可行的方法。
設x是一個大於0.5的數,則x可以表示為下列形式:
x=m*2e 
式中,0.5<=m<=1.0,e為整數。則求x的對數可以表示為:
    log2(x)=log2(m*2e)=log2(m)+log2(2e)=e+log2(m)
也就是說,求x的對數實際上只要求m的對數就可以了,而由於m的數值在0.5和1.0之間,用上面介紹的方法是完全可以實現的。例如:
    log2(10000)=log2(0.61035*214)=log2(0.61035)+14 =13.2877
可見,如果一個數可以用比較簡便的方法表示為上面的形式,則求任意大小數的對數也比較方便的。TMS320C2X/C5X指令集提供了一條用於對ACC中的數進行規格化的指令NORM,該指令的作用就是使累加器中的數左移,直至數的最高位被移至累加器的第30位。例如,對數值10000進行規格化的TMS320C25程式為。
LAC #10000
SACL TEMP
ZALH TEMP
LAR AR1,#0FH
RPT 14
NORM * -
    上述程式執行後,AR1=#0eH,ACCH=2000(10進位制)。對一個16位整數x進行上述程式處理實際上就是作這樣一個等效變換:
    x=[(x*2e)/32768]*215-Q
其中,暫存器AR1包含的值為15-Q累加器ACC高16位包含的值為x.2Q,其數值在16384至32768之間。
例1.14實現以2為底的對數的C定點模擬程式。
int logtab0[10]={-32768,-28262,-24149,-20365,-16862,
-13600),-1O549,-7683,-4981,-2425};/*Q15*/
int logtab1[10]={122529,20567,18920,175l7,16308,
15255,14330,13511,12780,12124};/*Q13*/
int logtab2[10]={16384,l8022,19660,21299,22938,
24576,26214,27853,29491,31130};/*Q15*/
int log2_fast(int Am)
{
int point,point1;
int index,x0,dx,dy,y;
point=0;
while(Am<16384){point++;Am=Am<<1;}/*對Am進行規格化*/
point1=(15-point-4)*512;/*輸入為Q4,輸出為Q9*/
index=((Am-16384)*20L)>>15;/*求查表地址*/
dx=Am-logtab2[index];
dy=((long)dx*logtab1[index])>>13;
y=(dy+longtab0[index])>>6;/*Q9*/
y=point1+y;
return(y);
}
    上述程式中,輸入值Am採用Q4表示,輸出採用Q9表示,如果輸入輸出的Q值與上面程式中的不同,則應作相應的修改。
    以上討論了DSP晶片進行定點運算所涉及的一些基本問題,這些問題包括:數的定標,DSP程式的定點模擬,DSP晶片的足點運算以及定點實現非線性函式的快速實現方法等。充分理解這些問題對於用定點晶片實現DSP演算法具有非常重要的作用。