深度學習FPGA實現基礎知識20(FPGA小數運算--DSP基礎“定點小數運算”)
需求說明：深度學習FPGA實現知識儲備
來自：第一部分 DSP基礎--定點小數運算
第二部分 dsp定點運算基本方法
整理來自：時間的詩
第一部分 DSP基礎--定點小數運算
在DSP世界中，由於DSP晶片的限制,經常使用定點小數運算。所謂定點小數，實際上就是用整數來進行小數運算。下面先介紹定點小數的一些理論知識，然後以C語言為例，介紹一下定點小數運算的方法。在TI C5000 DSP系列中使用16位元為最小的儲存單位，所以我們就用16位元的整數來進行定點小數運算。 先從整數開始，16位元的儲存單位最多可以表示0x0000到0xffff，65536種狀態，如果它表示C語言中的無符號整數的話，就是從0到65535。如果需要表示負數的話，那麼最高位就是符號位，而剩下的15位可以表示32768種狀態。這裡可以看出，對於計算機或者DSP晶片來說，符號並沒有什麼特殊的儲存方式，其實是和數字一起儲存的。為了使得無論是無符號數還是符號數，都可以使用同樣的加法減法規則，符號數中的負數用正數的補碼錶示。
我們都知道-1 + 1 =0，而0x0001表示1，那麼-1用什麼來表示才能使得-1 + 1 =0呢？答案很簡單：0xffff。現在就可以開啟Windows的計算器，用16進位制計算一下0xffff+0x0001，結果是0x10000。那麼0x10000和0x0000等價麼，我們剛才說過用16位元來表達整數，最高位的1是第17位，這一位是溢位位，在運算暫存器中沒有儲存這一位，所以結果是低16位，也就是0x0000。現在我們知道負數的表達方式了。舉個例子：-100。首先我們需要知道100的16進位制，用計算器轉換一下，可以知道是0x0064，那麼-100就是0x10000
- 0x0064，用計算器算一下得0xff9c。還有一種簡單的轉換符號的方法，就是取反加一：把數x寫成二進位制格式，每位0變1，1變0，最後把結果加1就是-x了。
好，複習了整數的相關知識之後，我們進入定點小數運算環節。所謂定點小數，就是小數點的位置是固定的。我們是要用整數來表示定點小數，由於小數點的位置是固定的，所以就沒有必要儲存它（如果儲存了小數點的位置，那就是浮點數了）。既然沒有儲存小數點的位置，那麼計算機當然就不知道小數點的位置，所以這個小數點的位置是我們寫程式的人自己需要牢記的。
先以10進製為例。如果我們能夠計算12+34=46的話，當然也就能夠計算1.2+3.4 或者 0.12+0.34了。所以定點小數的加減法和整數的相同，並且和小數點的位置無關。乘法就不同了。 12*34=408，而1.2*3.4=4.08。這裡1.2的小數點在第1位之前，而4.08的小數點在第2位之前，小數點發生了移動。所以在做乘法的時候，需要對小數點的位置進行調整？！可是既然我們是做定點小數運算，那就說小數點的位置不能動！！怎麼解決這個矛盾呢，那就是捨棄最低位。 也就說1.2*3.4=4.1，這樣我們就得到正確的定點運算的結果了。所以在做定點小數運算的時候不僅需要牢記小數點的位置，還需要記住表達定點小數的有效位數。上面這個例子中，有效位數為2，小數點之後有一位。
現在進入二進位制。我們的定點小數用16位二進位制表達，最高位是符號位，那麼有效位就是15位。小數點之後可以有0 - 15位。我們把小數點之後有n位叫做Qn，例如小數點之後有12位叫做Q12格式的定點小數，而Q0就是我們所說的整數。
Q12的正數的最大值是 0 111 . 111111111111，第一個0是符號位，後面的數都是1，那麼這個數是十進位制的多少呢，很好運算，就是 0x7fff / 2^12 = 7.999755859375。對於Qn格式的定點小數的表達的數值就它的整數值除以2^n。在計算機中還是以整數來運算，我們把它想象成實際所表達的值的時候，進行這個運算。
反過來把一個實際所要表達的值x轉換Qn型的定點小數的時候，就是x*2^n了。例如 0.2的Q12型定點小數為：0.2*2^12 = 819.2，由於這個數要用整數儲存， 所以是819 即 0x0333。因為捨棄了小數部分，所以0x0333不是精確的0.2，實際上它是819/2^12 =0.199951171875。
我們用數學表示式做一下總結：
x表示實際的數（*一個浮點數）， q表示它的Qn型定點小數（一個整數）。
q = (int) (x * 2^n)
x = (float)q/2^n
由以上公式我們可以很快得出定點小數的+-*/演算法：
假設q1，q2，q3表達的值分別為x1，x2，x3
q3 = q1 + q2 若 x3 = x1 + x2
q3 = q1 - q2 若 x3 = x1 - x2
q3 = q1 * q2 / 2^n若 x3 = x1 * x2
q3 = q1 * 2^n / q2若 x3 = x1 / x2
我們看到加減法和一般的整數運算相同，而乘除法的時候，為了使得結果的小數點位不移動，對數值進行了移動。
用c語言來寫定點小數的乘法就是：
short q1,q2,q3;
....
q3=((long q1) * (long q2)) >> n;
由於/ 2^n和* 2^n可以簡單的用移位來計算，所以定點小數的運算比浮點小數要快得多。下面我們用一個例子來驗證一下上面的公式：
用Q12來計算2.1 * 2.2，先把2.1 2.2轉換為Q12定點小數：
2.1 * 2^12 = 8601.6 = 8602
2.2 * 2^12 = 9011.2 = 9011
(8602 * 9011) >> 12 = 18923
18923的實際值是18923/2^12 = 4.619873046875 和實際的結果 4.62相差0.000126953125，對於一般的計算已經足夠精確了。
第二部分 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定點算術運算
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演算法具有非常重要的作用。