2014年9月23日火曜日

Sinテーブルを極限まで(?)小さくする。

なんとなくtiny13でどーにかしようとしているドM野郎の冨塚です。

要は「マルツtiny85だと秋月tiny13が4つ買える。」ってことです。
(メモリは8倍なので得なのだけど)

全然電子工作エリアではないのですが
 如何にゼニをかけないかという点で
  セコセコと最小構成で遊びたいわけです。

○三角関数系でいかにメモリを小さくするかで色々考えてたら思い出しました。

以下のネタは、1980年代の8bitパソコンで3次元ワイヤーフレームアニメーションをやるため、いかに小さなメモリで三角関数を計算するか考えてマシン語で組んでたのがルーツです。
(なんせFM-7のサブシステム共有メモリは128バイトしかなかったのです、RGBの一枚を捨てて、VRAM上でプログラムを走らせてました。)

☆テーブルは60バイトでなんとかなっちまう。
 60バイトってのが重要で、tiny13のEEPROMエリアが64バイト。
 つまりフラッシュメモリの1kbyteをプログラム領域としてフルに使用できるのね。

☆0度から360度の1度毎のdegで精度はほぼ16bitです。

○ 下準備
1、0度から90度までのSin値を32768倍したものを10度毎に
  20バイトのメモリ上に置きます。
  (なぜこいつだけ32768にするのかは後述)

2、0度から9度までのSin値を65536倍したものを1度毎に
  20バイトのメモリ上に置きます

3、0度から9度までのCos値を65536倍したものを1度毎に
  20バイトのメモリ上に置きます

C言語的に書くとこんなね。数値はホンマにつかえます。

uint16_t sin90[10] = { 0,5690,11207,16383,21062,25101,28377,30791,32270,32768};
uint16_t sin10[10] = { 0,1143,2287,3429,4571,5711,6850,7986,9120,10252};
uint16_t cos10[10] = { 0,65526,65496,65446,65376,65286,65176,65047,64898,64729};


 このunsigned16bitの30配列だけで0-360度の三角関数をなんとかします。
まずSin値を求めたい角度の1の位を捨ててφとし、1の位をθとします。

4、先ずSinφとCosφを求めます。
○10度毎のSinφテーブルは90度までしかないのですが
100-180度はSin(180-φ) 190-270度は-Sin(φ-180) 280-350度は-Sin(360-φ)
で求めます。

○10度毎のCosφは上記に90度足してあげれば求まります。

5、次に1の位のSinθとCosθを求めます。
これはそれぞれ0度から9度までしかないので、テーブルを参照すればOk

6、この公式を使います。
     Sin(φ+θ) = Sinφ×Cosθ+Cosφ×Sinθ
 これで0度から360度のSin値が求まるわけです。
 
☆なぜ 10度毎のSinφの値が32768倍なのか。
 つまるところ整数値なので
 (Sinφ×65536)×(Cosθ×65536)+(Cosφ×65536)×(Sinθ×65536) / 65536;
   でいいのですが、これやっちゃうと計算中に32ビットで桁溢れします。
 まあlong64でやればいいんですけど、8bitCPUで64ビット演算なんかあんまりしたくないしバイナリコードがでかくなったら本末転倒。
 精度はオチるがSinφのテーブルは32768倍にして
 (Sinφ×32768)×(Cosθ×65536)+(Cosφ×32768)×(Sinθ×65536) / 32768;
 とやることでint32で計算できるというわけです。

2π=360でなく2π=256でやるなら
上位5bit(16バイト)下位3bit(32バイト)計48バイトでいけるか?

サンプルはgccでコンパイルできる普通のcのソースです。
※tiny13用ではありません。

#include <stdio.h>
#include <math.h>
#define PI 3.141592

unsigned short  sin90[10] = {  0, 5690, 11207, 16383, 21062, 25101, 28377, 30791, 32270, 32768};
unsigned short  sin10[10] = { 0,1143,2287,3429,4571,5711,6850,7986,9120,10252};
unsigned short  cos10[10] = { 0,65526,65496,65446,65376,65286,65176,65047,64898,64729};

int sin36(int d ) {
    d = d % 36;
    if ( d <= 9 ) return sin90[d];
    if ( d <= 18 ) return sin90[18-d];
    if ( d <= 27 ) return -sin90[d-18];
    if ( d <= 36 ) return -sin90[36-d];
}

int cos36(int d) {
    return sin36(d+9);
}

int sin16(int d) {
    int d36,d10,value ;
    int s36,s010,c36,c010;
    d = d % 360;
    d10 = d % 10;
    d36 = ( d - d10 ) / 10 ;
    s36 = sin36(d36) ;
    if ( d10 == 0 ) return s36*2;
    c36 = cos36(d36) ;
    s010 = (int)sin10[d10] ;
    c010 = (int)cos10[d10] ;
    value = (s36 * c010 + c36 * s010) /32768 ;
    return (int) value ;
}

int cos16(int d) {
    return sin16(d+90);
}

int main(void)
{
    int d , s16 , c16;
    double deg , rad , rsin , tsin, rcos ,tcos;
    double dsin ,dcos,dmaxs,dmaxc ;
   
    /*どのくらいの精度が出るのか、math.hと比較*/
    dmaxs = dmaxc = 0.0;
    for ( d = 0 ; d < 360 ; d+=1 ) {
        deg = (double) d;
          rad = deg * PI / 180.0;
          rsin = sin(rad);
        s16 = sin16(d) ;
        tsin = (double)s16 / 65536.0 ;
          rcos = cos(rad);
        c16 = cos16(d) ;
        tcos = (double)c16 / 65536.0 ;
        dsin = rsin - tsin ;
        if ( dsin < 0 ) dsin = -dsin ;
        if ( dmaxs < dsin ) dmaxs = dsin ;
        dcos = rcos - tcos ;
        if ( dcos < 0 ) dcos = -dcos ;
        if ( dmaxc < dcos ) dmaxc = dcos ;
        printf("sin(%.3d)  math.h[%1.6f]table[%1.6f]Gosa[%1.6f] cos(%.3d) math.h[%1.6f]table[%1.6f]Gosa[%1.6f]\n"

        ,d,rsin,tsin,dsin,d,rcos,tcos,dcos);
    }
    printf("max Gosa sin:%1.6f  cos:%1.6f\n",dmaxs,dmaxc);
      return 0;
}


3 件のコメント:

  1. 見ただけではどういう理屈か読みきれないですが
    Sin波生成はこういう方法もあるのか

    ダメ元でVisual Studioでビルドしましたが
    動きました



    返信削除
  2. 今ではあまり使いどころのない小手先技です。
    VSでもいけましたかー。

    速度はオチますが他のでもメモリが不足している時にはテーブル節約に使えると思います。
    応用範囲としてはもっと細かくしたいとか…遅い周波数の波形を奇麗にするために0.1度毎のSinを求めたい場合など、0度から9度のSinCosテーブルを0.1度刻みに100配列づつとすれば、合計420バイトのテーブルで賄えるわけです。

    昔はSRAMが高くて、名古屋大須のアメ横で6264という8kbyteのメモリを買うために数日分のバイト代をつぎ込んだものです。

    さらには8bitのCPUは加算減算はできても頑張らないと乗算はできず、除算は少なくともビット数分ループせにゃならんので時間がかかるという代物でした。富士通FM-7で採用されていたモトローラの6809という8bitCPU(68000の前身)だけは8x8bit乗算が1命令でできたので、この方式でワイアフレームアニメーションができたのです。

    返信削除
  3. 調子づいて同じアルゴリズムで作った、円周を1024等分、組み込み向け2π=1024版です。
    一般には2π=256でテーブルが作られますがその4倍精度の16bitSinです。
    参照テーブルの配列はshortの(17+16+16)で49配列98バイト。
    (考えてみればテーブル先頭0を省略してもいけますので、46配列92バイトにできるか。)

    遅い周波数の波形も滑らかに…出るか?

    #include
    #include
    #define PI 3.141592

    unsigned short s64[17]={ 0, 3211, 6392, 9512, 12539, 15446, 18204, 20787,
    23170, 25330, 27245, 28898, 30273, 31357, 32138, 32610, 32767};
    unsigned short s16[16]={ 0, 402, 804, 1206, 1608, 2010, 2412, 2814,
    3215, 3617, 4018, 4420, 4821, 5222, 5622, 6023};
    unsigned short c16[16]={ 0, 65534, 65531, 65524, 65516, 65505, 65491, 65475,
    65457, 65436, 65412, 65386, 65358, 65327, 65294, 65258};


    int sin64(int d ) {
    d = d % 64;
    if ( d <= 16 ) return s64[d];
    if ( d <= 32 ) return s64[32-d];
    if ( d <= 48 ) return -s64[d-32];
    if ( d <= 64 ) return -s64[64-d];
    }

    int cos64(int d) {
    return sin64(d+16);
    }

    int sin16(int d) {
    int d64,d16,value ;
    int s64,s016,c64,c016;
    d = d % 1024;
    d16 = d % 16;
    d64 = ( d - d16 ) / 16 ;
    s64 = sin64(d64) ;
    if ( d16 == 0 ) return s64*2;
    c64 = cos64(d64) ;
    s016 = (int)s16[d16] ;
    c016 = (int)c16[d16] ;
    value = (s64 * c016 + c64 * s016) /32768 ;
    return (int) value ;
    }

    int cos16(int d) {
    return sin16(d+256);
    }

    int main(void)
    {
    int d , s16 , c16;
    double deg , rad , rsin , tsin, rcos ,tcos;
    double dsin ,dcos,dmaxs,dmaxc ;

    /*どのくらいの精度が出るのか、math.hのSinと比較*/
    dmaxs = dmaxc = 0.0;
    for ( d = 0 ; d < 1024 ; d+=1 ) {
    deg = (double) d;
    rad = deg * PI / 512.0;
    rsin = sin(rad);
    s16 = sin16(d) ;
    tsin = (double)s16 / 65536.0 ;
    rcos = cos(rad);
    c16 = cos16(d) ;
    tcos = (double)c16 / 65536.0 ;
    dsin = rsin - tsin ;
    if ( dsin < 0 ) dsin = -dsin ;
    if ( dmaxs < dsin ) dmaxs = dsin ;
    dcos = rcos - tcos ;
    if ( dcos < 0 ) dcos = -dcos ;
    if ( dmaxc < dcos ) dmaxc = dcos ;
    printf("sin(%.3d) math.h[%1.6f]table[%1.6f]Gosa[%1.6f] cos(%.3d) math.h[%1.6f]table[%1.6f]Gosa[%1.6f]\n",d,rsin,tsin,dsin,d,rcos,tcos,dcos);
    }
    printf("max Gosa sin:%1.6f cos:%1.6f\n",dmaxs,dmaxc);
    return 0;
    }


    返信削除