2015年10月8日木曜日

固定小数点の四則演算のテスト

波形に振幅変調をかけようと思うと小数点演算が必要になりそうなので、Visual Studioで固定小数点演算のテストをしてみた。「音遊び!BlackfinDSP基板でディジタル信号処理初体験 2015年 04月号[雑誌]インターフェース増刊」の記事を参考にした。Amazonではまた転売価格になってますが、マルツにはまだ在庫ありそうです。

<fixedpoint.h>
 #ifndef _FIXEDPOINT_H_  
 #define _FIXEDPOINT_H_  
 #include <stdint.h>  
 typedef int32_t fp32;  
 typedef int64_t fp64;  
 #define FIXQ                    (24u)  
 #define double_to_fp32(x)     ((fp32)((1u<<FIXQ)*((double)x)))  
 #define fp32_to_double(x)     (((double)(x))/(1u<<FIXQ))  
 // 四則演算  
 // オーバーフローに注意  
 // 特に除算の除数  
 #define fp32_add(x,y)          ((x)+(y))       
 #define fp32_sub(x,y)          ((x)-(y))  
 #define fp32_mul(x,y)          ((fp32)(((fp64)(x)*(fp64)(y))>>FIXQ))  
 #define fp32_div(x,y)          ((fp32)(((fp64)(x)<<FIXQ)/(y)))  
 #endif //_FIXEDPOINT_H_  

PSoCでは関数呼び出しのコストが高そうなので、マクロ関数で定義してみた。

テストプログラム
 #include "stdafx.h"  
 #include <stdio.h>  
 #include <stdlib.h>  
 #include "fixedpoint.h"  
 #define CALC_DIV     1  
 #define CALC_REPEAT     0xffff     // Excelの最大行  
 #define     FP32_MAX     (1.0f)  
 double dbl_rand()  
 {  
      double dr = (double)rand() / RAND_MAX;  
      return (dr - 0.5) * FP32_MAX;  
 }  
 void calc(double dx, double dy)  
 {  
      fp32 fpx = double_to_fp32(dx);  
      fp32 fpy = double_to_fp32(dy);  
      double x_v = fp32_to_double(fpx);  
      double y_v = fp32_to_double(fpy);  
      double add_v = fp32_to_double(fp32_add(fpx, fpy));  
      double sub_v = fp32_to_double(fp32_sub(fpx, fpy));  
      double mul_v = fp32_to_double(fp32_mul(fpx, fpy));  
 #if CALC_DIV  
      double div_v = fp32_to_double(fp32_div(fpx, fpy));  
 #endif  
      printf("x\t%lf\t%lf\t%lf\n", dx, x_v, dx - x_v);  
      printf("y\t%lf\t%lf\t%lf\n", dy, y_v, dy - y_v);  
      printf("x+y\t%lf\t%lf\t%lf\n", dx + dy, add_v, (dx + dy) - add_v);  
      printf("x-y\t%lf\t%lf\t%lf\n", dx - dy, sub_v, (dx - dy) - sub_v);  
      printf("x*y\t%lf\t%lf\t%lf\n", dx * dy, mul_v, (dx * dy) - mul_v);  
 #if CALC_DIV  
      printf("x/y\t%lf\t%lf\t%lf\n", dx / dy, div_v, (dx / dy) - div_v);  
 #endif  
 }  
 int _tmain(int argc, _TCHAR* argv[])  
 {  
      printf("expr\tdouble\tfp32\terr\n");  
      calc(1.0f, -1.0f);  
      for (int i = 0; i < CALC_REPEAT; i++) {  
           double dx = dbl_rand();  
           double dy = dbl_rand();  
           // 除算エラー回避  
           if (dy == 0.0f)  
                continue;  
           calc(dx, dy);  
      }  
      return 0;  
 }  

Q8.24の固定小数点と浮動小数点で、-1.0~1.0の範囲で乱数出力して誤差を出力した。

除算あり

E列、F列はExcelの関数で誤差の最大値と最小値を抽出した。かなりごつい誤差がでているが、これは除算の時に除数が小さい時に発生する。

除算なし

プログラムを修正して除算を省いてみると誤差はなくなった。ただし出力桁数が少ないのでもっと細かい誤差はあると思うが、DDSで使うならこれぐらいの精度で十分だと思う。

Visual Studioのプロジェクト
https://github.com/ryood/FixedPoint_Test