一種三角函數(shù)的近似計算方法
【專利摘要】本發(fā)明涉及一種三角函數(shù)的近似計算方法，該方法包括：構建任意角P的近似三角函數(shù)C的計算公式為：C=(215-1)×cos(P)=(215-1)×cos(m×2π/1024)-n×2π/65536×(215-1)×sin(m×2π/1024)以m為地址從存儲單元中查找(215-1)×cos(m×2π/1024)的數(shù)值A并以n為地址從存儲單元中查找(215-1)×sin(m×2π/1024)的數(shù)值B；將n×2π/65536與數(shù)值B相乘得到乘積；以及將數(shù)值A減去所述乘積得到近似三角函數(shù)C。本發(fā)明所述方法優(yōu)點在于：可減小粗大偏差，降低存儲器使用量；采用簡單的加法和乘法計算，減少運行時間，在保持高精度的基礎上，提高運算速度，本方法可在短時間內將正、余弦計算精度達到16bit以上；本方法還可以擴展到計算坐標旋轉、雙曲函數(shù)等復雜函數(shù)的計算。
【專利說明】一種三角函數(shù)的近似計算方法
【技術領域】
[0001]本發(fā)明涉及一種計算的方法，特別是一種三角函數(shù)的近似計算方法。
【背景技術】
[0002]目前，F(xiàn)PGA(在線可編程邏輯陣列)已經廣泛應用于各種數(shù)字信號處理領域中，用于各種數(shù)學計算、信號濾波計算等，其中三角函數(shù)計算是數(shù)學計算中的一個重要方面。由于FPGA的內部存儲器件數(shù)量和數(shù)據(jù)總線寬度的限制，一些數(shù)值或數(shù)值的計算結果均往往只能以近似形式表示，結果準確性較差。比如常見的三角函數(shù)計算，由于受到存儲器容量和數(shù)據(jù)總線寬度的限制，通常只能用IObit數(shù)據(jù)線寬來近似表示，近似的結果與實際的準確值有2-10(約為0.1% )左右的相對偏差，這樣的偏差在很多場合(比如調頻波的計算)是不能接受的。如果想采用較寬的數(shù)據(jù)線寬表示，則又可能面臨內部存儲空間不夠的問題，因為計算時存儲器件數(shù)量大致是以線寬的2的指數(shù)變化，即線寬每增加一位，存儲器增加一倍，數(shù)據(jù)精度也提高一倍。
[0003]目前三角函數(shù)計算方案大致有兩種:一種方案采用直接數(shù)據(jù)查表的方法，將三角函數(shù)值的結果存儲于FPGA內部存儲器件上。地址對應角度，存儲內容代表三角函數(shù)值，這種方式速度快，只需一個時鐘脈沖即可得到結果，但精度不高；另一種方案是CORDIC算法，這種方案使用內存較少，計算結果精度相對較高一些，但需要很多個時鐘脈沖，重復迭代計算，精度越高，需要的時鐘脈沖越多。在對計算 時間要求較高的情況下，很難滿足要求。

【發(fā)明內容】

[0004]本發(fā)明要解決的技術問題是提供一種三角函數(shù)的近似計算方法，以解決數(shù)字信號處理技術中用FPGA實現(xiàn)三角函數(shù)計算的精度低和速度慢的問題。
[0005]為解決上述問題本發(fā)明提供一種三角函數(shù)的近似計算方法，對于角度P，
[0006]P=mX 2 π /1024+ηΧ 2 π /65536,其中 m=0 ~1023, n=0 ~63;
[0007]該角度P的近似三角函數(shù)C的計算公式為:
[0008]C=(215-l) Xcos(P) = (215-1) Xcos(mX2 π /1024) -ηΧ2 π /65536 X (215_1) Xsin(mX2 3? /1024)
[0009]該方法包括以下步驟:
[0010]以m為地址從存儲單元中查找(215-1) Xcos(mX2 31 /1024)的數(shù)值A并以η為地址從存儲單元中查找(215-1) Xsin(mX2Ji/1024)的數(shù)值B ;
[0011]將η X 2 31 /65536與數(shù)值B相乘得到乘積；以及
[0012]將數(shù)值A減去所述乘積得到近似三角函數(shù)C。
[0013]優(yōu)選的，該方法以一個時鐘脈沖完成數(shù)值查找步驟，以及以一個時鐘脈沖完成乘法計算步驟。
[0014]一種用于近似計算三角函數(shù)的FPGA，該FPGA包括
[0015]函數(shù)表儲存模塊，包括1024X16bit存儲區(qū)，用于存儲(215-1) X cos (mX 2 3i /1024)的數(shù)值 A 和(215_1) X sin (mX 2 π/1024)的數(shù)值 B，其中 m=0 ~1023 ；
[0016]數(shù)據(jù)讀取單元，以m為地址從所述存儲模塊中查找(215_1) X cos (mX 2 π /1024)的數(shù)值A并以η為地址從存儲單元中查找(215-1) Xsin(mX2 3i/1024)的數(shù)值B，其中n=0~63
[0017]乘法器，用于將nX2 /65536與數(shù)值B相乘得到乘積；
[0018]加法器，用于將數(shù)值A減去所述乘積得到近似三角函數(shù)C，
[0019]C=(215_l) Xcos(P)，=(215_1) Xcos (mX2 π /1024) -ηΧ2 π /65536 X (215_1) Xsin(mX2 3? /1024)
[0020]其中，P=mX23? /1024+nX2 π /65536。
[0021]優(yōu)選的，該FPGA包括輸入和輸出線寬均為16bit的數(shù)據(jù)線。
[0022]優(yōu)選的，該所述乘法器的輸出數(shù)據(jù)線線寬為16bit。
[0023]優(yōu)選的，所述數(shù)據(jù)讀取單元以一個時鐘脈沖讀取數(shù)據(jù)；所述乘法器以一個時鐘脈沖計算乘積。根據(jù)本發(fā)明所述方法優(yōu)點在于:
[0024]1、可減小粗大偏差，降低存儲器使用量；
[0025]2、采用簡單的加法和乘法計算，減少`行時間，在保持高精度的基礎上，提高運算速度，本方法可在短時間內將正、余弦計算精度達到16bit以上；
[0026]3、本方法還可以擴展到計算坐標旋轉、雙曲函數(shù)等復雜函數(shù)的計算。
【專利附圖】

【附圖說明】
[0027]圖1示為一種三角函數(shù)的近似計算方法示意圖；
[0028]圖2示為一種近似計算三角函數(shù)的小型FPGA的示意圖。
[0029]1、函數(shù)表儲存模塊，2、數(shù)據(jù)讀取單元，3、乘法器，4、加法器。
【具體實施方式】
[0030]本發(fā)明提供一種三角函數(shù)的近似計算方法，該方法包括以下步驟:對于角度P，P=mX2 31 /1024+nX2 31 /65536，其中m=0~1023，n=0~63;該角度P的近似三角函數(shù)C的計算公式為:
[0031]C=(215-l) Xcos(P) (SI) ;=(215_1) Xcos(mX2 π /1024)-ηX2 π /65536X (215_1)X sin (mX 2 3? /1024)
[0032]以m為地址從存儲單元中查找(215-1) Xcos(mX2 π/1024)的數(shù)值A并以η為地址從存儲單元中查找(215-1) Xsin(mX2 /1024)的數(shù)值B (S2);將nX2 π/65536與數(shù)值B相乘得到乘積，以及將數(shù)值A減去所述乘積得到近似三角函數(shù)C (S3)。
[0033]本發(fā)明提供一種計算三角函數(shù)的FPGA，該FPGA包括函數(shù)表儲存模塊1、數(shù)據(jù)讀取單元2、乘法器3和加法器4 ;所述函數(shù)表儲存模塊I用于查詢并調取正、余弦值；所述乘法器3用于通過數(shù)據(jù)讀取單元2從所述函數(shù)表儲存模塊中調取正弦值與偏移量相乘；所述加法器4用于通過數(shù)據(jù)讀取單元2從所述函數(shù)表儲存模塊中調取余弦值減去正弦值與ηΧ2 31 /65536 的乘積。
[0034]下面結合附圖對本發(fā)明做進一步描述。[0035]本發(fā)明所述的一種計算三角函數(shù)的小型FPGA包括函數(shù)表儲存模塊1、數(shù)據(jù)讀取單兀2、乘法器3和加法器4,輸入輸出均米用16bit表不。以下將以計算一個O~2π角度的余弦值為例對本發(fā)明作進一步描述。
[0036]函數(shù)表儲存模塊I是一個1024X16bit的存儲區(qū)，地址用IObit數(shù)據(jù)線寬，余弦值用16bit表示,地址為a的位置,余弦值為b,那么a、b的關系為:
[0037]b=(215_l) Xcos (a/1024X2 π )，a=0, I, 2，…，1023 ;這樣就得到一個具有 1024個存儲單元的基本正、余弦表，利用這個表也可以查出正、余弦值。這個表的最小相位分辨力只有2 /1024rad.，為了提高計算精度，本發(fā)明將原來的相位最小分辨量2 π /1024rad再平均分成64 (即26)等分，這時O~2JI角度就需要用16bit (1024X64=210X26=216)表示，分辨力就可達到2 31 /65536 rad.帶來的效果是所表示相位的準確度由2 π /1024rad提高到2 π/65536 rad，即準確度提高64倍。而這個改進需要的硬件只是增加一個6bit (64=26)的存儲器，用來存儲16bit相位量值中的最低的6bit。用來參與下面的計算之用。
[0038]此時需要計算的任意一個相角P都可以表示為:
[0039]P=mX2 31 /1024+ηΧ2 π /65536，m=0, I, 2，…，1023 ；n=0, I, 2，...，63 ;即將相角表示為一個2 31 /1024倍數(shù)的整數(shù)部分的角度mX 2 31 /1024和小于2 31 /1024的小數(shù)部分的角度nX2 /65536的和。而整數(shù)部分的角度mX 2 π /1024的正、余弦可以利用前面提到的1024個三角函數(shù)存儲單元I得到；小數(shù)部分的角度可以利用泰勒公式在整數(shù)部分的角度附近展開得到。因此，任意相角P的三角函數(shù)C的計算公式可表示為:C=(215-1)XCoMP) = (215-1) X cos (mX 2 31 /1024) -η X 2 π /65536 X (215_1) X sin (mX 2 π /1024)。實際取值中，根據(jù)要獲取的精度要求僅取展開式的前兩項即可，隨后的所有項對最終結果不會產生嚴重影響。
[0040]其中(215-1)X cos (`mX 2 31 /1024)和(215_1) Xsin(mX2 π /1024)可以利用數(shù)據(jù)讀取單元2從函數(shù)表存儲器模塊中分別獲得數(shù)值A和數(shù)值B。這樣三角函數(shù)C的結果就可以利用數(shù)值A減去數(shù)值B與nX 2 /65536的乘積獲得。全部過程包括一次查表、一次相乘、一次相加，分別需要I個時鐘脈沖、I個時鐘脈沖、O個時鐘脈沖?？偣残枰?個時鐘脈沖即可完成，相位精度達到2Ji/65536rad.，計算結果的精度達到2/65536以上(含舍去的高階項)。表示相位的數(shù)據(jù)線寬度和計算結果使用的數(shù)據(jù)線寬度為16bit ;過程中間乘法計算線寬24bit，得到最終結果時略去低Sbit (等效數(shù)學計算時略去小數(shù)位一樣)，占用存儲區(qū):1024X16bit加上少量乘法計算存儲區(qū)。如果完全采用查表方式，則需要I個時鐘完成，相位精度達到2 Ji/65536rad，計算結果的精度達到1/65536以上，使用線寬16bit，占用存儲區(qū):65536X16bit0
[0041]如果采用CORDIC算法，存儲器用量少，要達到計算結果的精度達到1/65536以上，理論上至少需要16個時鐘脈沖，實驗中發(fā)現(xiàn)需要27個時鐘脈沖，速度大為降低。并且中間的乘法運算也需要非常寬(32bit以上)的數(shù)據(jù)線寬，這樣的數(shù)據(jù)線寬在小型FPGA中也是難以實現(xiàn)的，因此，利用現(xiàn)有方法難以在小型FPGA中實現(xiàn)65536父1613^大小的存儲和較寬的數(shù)據(jù)線寬。而本發(fā)明所述方法及計算單元可大幅度提高余弦函數(shù)的計算精度和計算速度，特別適合在小型FPGA上實現(xiàn)高精度快速計算。
[0042]實際調頻時，如果相對頻偏量為500時，需要計算cos [500*COS (P)]，上例中查表所帶來的偏差被放大500倍，再進行一次余弦計算，所得結果的偏差將不可接受，而本發(fā)明通過泰勒級數(shù)展開的方法，結合低精度的正(余)弦表，將查表和計算結合起來實現(xiàn)高精度的正、余計算，既具備有查表的速度又有計算帶來的簡便及少量存儲器占用。這樣整個計算既提高了計算精度又保持了計算速度，同時還節(jié)省了存儲器。實施例1
[0043]相角P=(2ji/1024)*128.25rad,如果直接查具有1024個存儲單元的基本余弦表時，因小于2Ji/1024rad的相位無法表示，首先必將相位截為(2 π/1024) *128rad，再查表得到 65536*cos (P) =46340(而真實值應為 46269)；
[0044]而采用本發(fā)明時，相位值可準確的表示為P= (2 31 /1024) *128.4= (2 /1024) *128+(2 31 /65536) *16 ;計算 65536*cos (P) =cos [ (2 π /1024) *128+ (2 π /65536) *16] ^ cos [ (2 π/1024)*128]-16X2 3? /65536*sin[ (2 π /1024)*128] =46340-16X2 π /65536*46340=46268 (而真實值應為46269)，結果要精確許多。
[0045]綜上所述，利用本發(fā)明所述方法可減小粗大偏差，降低存儲器使用量；基于簡單的加法和乘法計算，減少運行時間，在保持高精度的基礎上，提高運算速度，同時還可在短時間內將正、余弦計算精度達到16bit以上；本發(fā)明不僅限于計算三角函數(shù)還可以擴展到計算坐標旋轉、雙曲函數(shù)等復雜函數(shù)的計算。
[0046]以上內容是結合具體的優(yōu)選實施方式對本發(fā)明所作的進一步詳細描述說明，不能認定本發(fā)明的具 體實施只局限于這些說明。對于本發(fā)明所述【技術領域】的普通技術人員來說在不脫離本發(fā)明構思的前提下，還可以做出若干簡單推演或替換，都應當視為屬于本發(fā)明的保護范圍。
【權利要求】
1.一種三角函數(shù)的近似計算方法，其特征在于:對于角度P，
P=mX2 31 /1024+nX2 31 /65536，其中 m=0 ~1023，n=0 ~63; 該角度P的近似三角函數(shù)C的計算公式為:
C=(215-l) Xcos(P) = (215-1) X cos (mX 2 3i /1024) -nX2 π /65536 X (215_1) Xsin(mX23i /1024) 該方法包括以下步驟: 以m為地址從存儲單元中查找(215-1) Xcos(mX2 /1024)的數(shù)值A并以η為地址從存儲單元中查找(215-1) X sin (mX 2 31 /1024)的數(shù)值B ； 將nX 2 /65536與數(shù)值B相乘得到乘積；以及 將數(shù)值A減去所述乘積得到近似三角函數(shù)C。
2.根據(jù)權利要求1所述的三角函數(shù)的近似計算方法，其特征在于:該方法以一個時鐘脈沖完成數(shù)值查找步驟，以及以一個時鐘脈沖完成乘法計算步驟。
3.一種用于近似計算三角函數(shù)的FPGA，其特征在于:該FPGA包括 函數(shù)表儲存模塊，包括1024X16bit存儲區(qū)，用于存儲(215-1) X cos (mX 2 Ji/1024)的數(shù)值 A 和(215-1) X sin (mX 2 31 /1024)的數(shù)值 B，其中 m=0 ~1023 ； 數(shù)據(jù)讀取單元，以m為地址從所述存儲模塊中查找(215-1) XC0S(mX2 π /1024)的數(shù)值A并以η為地址從存儲單元中查找(215-1) Xsin(mX2 3i/1024)的數(shù)值B，其中n=0~63乘法器，用于將nX2 /65536與數(shù)值B相乘得到乘積； 加法器，用于將數(shù)值A減去所述乘積得到近似三角函數(shù)C，
C=(215-l) Xcos(P)，=(215_1) Xcos (mX2 π /1024) -ηΧ2 π /65536 X (215_1) Xsin (mX2 3 /1024) 其中，P=mX2 3? /1024+nX2 3i /65536。
4.根據(jù)權利要求3所述的FPGA，其特征在于:該FPGA包括輸入和輸出線寬均為16bit的數(shù)據(jù)線。
5.根據(jù)權利要求4所述的FPGA，其特征在于:該所述乘法器的輸出數(shù)據(jù)線線寬為16bit0
6.根據(jù)權利要求3所述的FPGA，其特征在于: 所述數(shù)據(jù)讀取單元以一個時鐘脈沖讀取數(shù)據(jù)； 所述乘法器以一個時鐘脈沖計算乘積。
【文檔編號】G06F17/17GK103699518SQ201310740098
【公開日】2014年4月2日 申請日期:2013年12月26日 優(yōu)先權日:2013年12月26日
【發(fā)明者】吳遠伍, 武衛(wèi)平, 李文意, 任士卿, 韓蕊, 程春悅, 陳晉龍 申請人:北京無線電計量測試研究所