浮點轉定點運算
阿新 • • 發佈:2018-12-08
一 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[i]=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語言浮點程式。
#include <stdio.h>
const int length=180
void filter(int xin[],int xout[],int n,float h[]);
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[i];
for(i=0;i<length;i++)
{
sum=0.0;
for(j=0;j<n;j++)sum+=h[j]*x1[i-j+n-1];
xout[i]=(int)sum;
for(i=0;i<(n-l);i++)x1[n-i-2]=xin[length-1-i];
}
void main()
FILE *fp1,*fp2;
int frame,indata[length],outdata[length];
fp1=fopen(insp.dat,"rb");
fp2=fopen(Outsp.dat,"wb");
frame=0;
while(feof(fp1) ==0)
{
frame++;
printf(“frame=%d\n”,frame);
for(i=0;i<length;i++)indata[i]=getw(fp1);
filter(indata,outdata,19,h);
for(i=0;i<length;i++)putw(outdata[i],fp2);
}
fcloseall();
return(0);
}
例1.8語音訊號800Hz l9點FIR低通濾波C語言定點程式。
#include <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};
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][i];
for(i=0;i<1ength;i++)
sum=0;
for(j=0;j<n;j++)sum+=(long)h[j]*x1[i-j+n-1];
xout[i]=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)
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[i]=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語言浮點程式。
#include <stdio.h>
const int length=180
void filter(int xin[],int xout[],int n,float h[]);
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[i];
for(i=0;i<length;i++)
{
sum=0.0;
for(j=0;j<n;j++)sum+=h[j]*x1[i-j+n-1];
xout[i]=(int)sum;
for(i=0;i<(n-l);i++)x1[n-i-2]=xin[length-1-i];
}
void main()
FILE *fp1,*fp2;
int frame,indata[length],outdata[length];
fp1=fopen(insp.dat,"rb");
fp2=fopen(Outsp.dat,"wb");
frame=0;
while(feof(fp1) ==0)
{
frame++;
printf(“frame=%d\n”,frame);
for(i=0;i<length;i++)indata[i]=getw(fp1);
filter(indata,outdata,19,h);
for(i=0;i<length;i++)putw(outdata[i],fp2);
}
fcloseall();
return(0);
}
例1.8語音訊號800Hz l9點FIR低通濾波C語言定點程式。
#include <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};
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][i];
for(i=0;i<1ength;i++)
sum=0;
for(j=0;j<n;j++)sum+=(long)h[j]*x1[i-j+n-1];
xout[i]=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)