のろのろのろ雑記

_ [ColdFire] ColdFire×P2P地震情報(A) PWMと圧電ブザーで音を鳴らす

　Interface 9月号付録基板を使い，P2P地震情報のサービスを作ってしまう勝手な企画です．が，番号がアルファベットのものはおまけです．

二番煎じに関する注意事項

　これは新適当マイコン電子工作研究所で行われていることの二番煎じです．が，PWMの出力をクリスタルイヤホンで…　などと見ていると，「これだ」と思ってそのまま突っ走ってしまったわけです．お付き合いください．

経緯

1. 「PWMの出力を…　クリスタルイヤホン？」というわけで，検索する．
2. 「ふむふむ，圧電効果とかいうのか…」
3. 「圧電？　圧電ブザーってあるね？　つーか同じじゃね？」
4. 「じゃあ圧電ブザー繋げて同じことすれば音鳴るんじゃね？」

材料

• ColdFire基板 (Interface 9月号付録，RJ-45とDCジャックを実装済み)
• 何に使ったか覚えていない導線
• 古いマシンから切り取った圧電ブザー
  • "HYCOM HY-05"と書かれていて，どうやら3.0 V〜8.0 Vで動作するらしい．
• 私の時間

手順1: 本当に鳴るのか

　まずは，MCF52233付録基板 - SilentCが使っていないモジュールを探せ(4)に記載されているとおりにやってみます．

　(PDF)リファレンスマニュアルや(PDF)回路図を見ると，PNQ[1]＝IRQ1のようなので，次のように接続します．

　すると，確かに「プツッ，プツッ，プツッ…」と3.6 Hzくらいで鳴っているのが聞こえます．いけるようです．

手順2: プリスケーラ？　スケーラ？　デューティ？　周波数を設定するには…

　圧電ブザーは高速でON/OFFを繰り返すパルス信号によって，その信号の周波数の音を出すようです．ですから，この周波数さえ設定できれば音の高低を自由自在に操れるわけです．

　が，WindowsのBeep関数のように「周波数」という設定項目はありません．あるのは，プリスケーラ，スケーラ，デューティ，意味の分からない項目ばかりです．じっくり向き合う必要があるようです．

手順3: ColdFireの"Pulse Width Modulation(PWM)"って？

　と言いつつ，PWMの解説は他所に譲ります．ここでは，ColdFireで用意されている設定項目について見ていきます．リファレンスマニュアルをご用意ください．

PWM Clock Select Register (PWMCLK) - 29.2.3/534ページ

PWMのクロック源を選ぶ．AやBを指定するとプリスケーラ(PWMPRCLK)のみ，SAやSBを指定するとプリスケーラとスケーラ(PWMSCLA/PWMSLCB)の両方を通る．

PWM Prescale Clock Select Register (PWMPRCLK) - 29.2.4/535ページ

プリスケーラの分周比を選ぶ．簡単に言うと，内部バスクロックの何分の1のクロックで動くようにするかを指定する．ただし，分母は2の0乗,1乗,2乗,…,7乗．

PWM Scale A/B Register (PWMSCLA/B) - 29.2.7〜/538ページ〜

プリスケーラクロックの何分の1で動くようにするかを指定する．ただし，分母は2,4,6,…,512．

PWM Channel Period Registers (PWMPERn) - 29.2.10/540ページ

PWMCLKなどで指定したクロックを元にカウントを行うが，そのカウンタの上限値，1サイクルに掛かるクロック数．

PWM Channel Duty Registers (PWMDTYn) - 29.2.11/541ページ

1サイクルに掛かるクロック数に対し，信号がONになるクロック数．PWMPERnと同じにすれば常にONだし，半分ならON/OFFがピッタリ交互だし，0ならOFFのまま．

　プリスケーラ，スケーラ，ピリオド（デューティ）を組み合わせることで，やっと周波数を指定出来るというわけです．これはなかなかしんどい．

　スケーラを通る（PWMCLKにSA/SBを指定した）場合の周波数の計算式を置いておきます．ただし，デューティ比（信号ONの割合）は50 %とします．

=   60*1000*1000   / ( 2 ^ PWMPRCLK ) / ( 2 * PWMSCLA/B ) / PWMPERn
  内部バスクロック     プリスケーラ          スケーラ       ピリオド

手順4: 音階に対応する組み合わせ

　早速組み合わせてみました．平均律を用いています．0.6 Hzくらいズレているところもありますが，お遊びということで．

　ポイントは，「プリスケーラを1減らすと，周波数が2倍＝1オクターブ上になる」こと，そして逆もまたあるということです．上記のような1オクターブ分の音階表を持っておけば，あとはプリスケーラの値だけで対応できます．

　プリスケーラは0〜7まで指定できるので，この場合だとG#2〜G 8までは自由に出せることになります．この範囲を外れてしまうと所々音が鳴らなくなります．

手順5: 実装

mel(){
  // char *meloは可読性のため改行を入れたので，実行時には削除すること
  char *melo="5E-85D#85E-85D#85E-84B-85D-85C-84A-83E-83A-84C-84E-84A-84B-83E-8
              3G#84E-84G#84B-85C-83E-83A-84E-85E-85D#85E-85D#85E-84B-85D-85C-8
              4A-83E-83A-84C-84E-84A-84B-83E-83G#84D-85C-84B-84A-83E-83A-8SS-8
              EE";
  char *pt=melo; int bt,slp;
  bt=10000/(16000/60); // テンポを指定し，4分音符の長さを10 ms単位で計算
                       // "10000 / (テンポ*10 / 60)"
                       // 小数計算が出来ないため簡易精度確保

  int *pnq = 0x40100068;  // PNQ[1]
  char *pwm = 0x401b0000; // PWM

  *pnq |= 12; pwm[2]=2;   // PNQ[1]とPWM1を接続，PWM1を有効に

  #stop 0
  for(;;){
    if(Getc(0)!=0)break;
    if(pt[0]==69)break;

    if(pt[1]==67&&pt[2]==45){pwm[3]=4;pwm[8]=256;pwm[21]=14;pwm[29]=7;} //C
    if(pt[1]==67&&pt[2]==35){pwm[3]=3;pwm[8]=199;pwm[21]=34;pwm[29]=17;} //C#
    if(pt[1]==68&&pt[2]==45){pwm[3]=3;pwm[8]=103;pwm[21]=62;pwm[29]=31;} //D
    if(pt[1]==68&&pt[2]==35){pwm[3]=3;pwm[8]=137;pwm[21]=44;pwm[29]=22;} //D#
    if(pt[1]==69&&pt[2]==45){pwm[3]=3;pwm[8]=39;pwm[21]=146;pwm[29]=73;} //E
    if(pt[1]==70&&pt[2]==45){pwm[3]=3;pwm[8]=179;pwm[21]=30;pwm[29]=15;} //F
    if(pt[1]==70&&pt[2]==35){pwm[3]=3;pwm[8]=181;pwm[21]=28;pwm[29]=14;} //F#
    if(pt[1]==71&&pt[2]==45){pwm[3]=5;pwm[8]=46;pwm[21]=26;pwm[29]=13;} //G
    if(pt[1]==71&&pt[2]==35){pwm[3]=2;pwm[8]=215;pwm[21]=42;pwm[29]=21;} //G#
    if(pt[1]==65&&pt[2]==45){pwm[3]=2;pwm[8]=147;pwm[21]=58;pwm[29]=29;} //A
    if(pt[1]==65&&pt[2]==35){pwm[3]=2;pwm[8]=149;pwm[21]=54;pwm[29]=27;} //A#
    if(pt[1]==66&&pt[2]==45){pwm[3]=2;pwm[8]=211;pwm[21]=36;pwm[29]=18;} //B
    if(pt[1]==83){*pwm=0;}else{*pwm=2;}
    pwm[3]-=pt[0]-51; //所詮ビープなので，指定より2オクターブ高く出して聞こえやすく

    if(pt[3]==52){slp=bt;}
    if(pt[3]==56){slp=bt/2;}
    if(pt[3]==57){slp=bt/3;}
    if(pt[3]==50){slp=bt*2;}
    if(pt[3]==49){slp=bt*4;}

    Sleep(slp);

    pt+=4;
  }
  *pwm=0;
  PrStr("E.\r\n");
}

　char *meloでメロディを作ります．4バイトで1個の音を奏でる独自形式です．

[オクターブ][音階][シャープ][長さ]
 ※音階はA〜GまたはS（休み）で，シャープは#または-（ありまたはなし）
 ※長さは 1（全音符），2（2分音符），4，8，9（3連8分音符）のみ

　forループ内で次のような処理を繰り返します．

1. オクターブが E であれば，ループを抜ける．
2. 音階とシャープを見て，プリスケーラなどを設定
3. 休みであればPWMオフ，そうでなければPWMオン
4. オクターブを見て，プリスケーラを加減
5. 長さを見て，演奏時間を計算してスリープ
6. 次の位置へポインタをセット

　若干順番が乱れていますが，気にしないことにします．かなり強引にやってしまいましたが，これを実行すると「エリーゼのために」が流れます．

　制限として，約256バイトに収める必要があるようです．それを超えたり，終了命令（オクターブに E）を忘れたりすると，範囲外にまで演奏を続行して変な音になります．

手順6: 実演

• ニコニコ動画: ColdFire基板と圧電ブザーで音を鳴らしてみた(YouTube版)

　SilentCからここまで操作できちゃうとは．プログラミングガイドを読んで「こんなものか」と思っていたのがバカみたいです．スバラシイ．

参考

• http://noritan-micon.blog.so-net.ne.jp/
• http://www.freescale.com/files/32bit/doc/ref_manual/MCF52235RM.pdf
• http://www.cqpub.co.jp/interface/download/2008/09/kairozu.pdf
• http://www.apache-band.com/apchhp/music_theory/frequency_temperament.html

_ [ColdFire] ColdFire×P2P地震情報(7) "33.3 Hzの壁"を超える

　Interface 9月号付録基板を使い，P2P地震情報のサービスを作ってしまう勝手な企画です．過去の記事は，タイトルのColdFireカテゴリからどうぞ．

前回のまとめ

• 計測震度とフィルタとの対応を見ました．
• 静止状態は震度3相当でした…

"33.3 Hzの壁"

　静止状態で震度3になってしまったということは，すなわち震度4以上でないと働かないということであって，それを震度計・地震計とは呼べないということです．

　何とかしてノイジーな状態を解消しなければなりません．ひとまず，"オーバーサンプリング"を考えてみます．オーバーサンプリングの効果は私も詳しく知りませんが，とりあえず震度3という表示にビックリしたので何でもやってみます．

loop_axis(){
 char t1; int i,cr;
 t1=CreateTimer(0,30000,0);
 InitAd(0x70);
 for(i=0;i<1000;i++){
  GetAd(4);GetAd(5);GetAd(6);
 }
 cr=GetTimerCount(t1);
 PrNum(cr);PrStr("\r\n");
 CreateTimer(t1,0,0);
}

Ax_Timer::loop_axis
26999

Ax_Timer::loop_axis
27000

　上記のコードは，3軸を1000回サンプリングした時の所要時間を計測するものです．タイマは10 ms単位ですから， "(30000 - RESULT) / 100" と計算すれば所要秒数が分かります．

　計算してみると，ちょうど30 s掛かっていることが分かります．サンプリングのためにGetAd関数を3軸×1000回＝3000回呼び出していますから，「もしかして，GetAd関数はキッチリ10 ms掛かるようになっているんじゃないか？」，つまり3軸取得時には『33.3 Hzの壁』があるのではないか，という疑問が出てきます．

"33.3 Hzの壁"を超える

　この壁を超えます．新適当マイコン電子工作研究所あたりを見ていると，どうやらメモリマップドI/Oのようです．ということは，A/Dコンバータ(ADC)とてメモリアクセスで簡単に出来るはず…

　リファレンスマニュアルを調べると，"Chapter 28 Analog-to-Digital Converter (ADC)"，493ページ以降に書かれています．

• 8チャネルの入力と，2つのA/Dコンバータ
  • AN0〜3はコンバータAに，AN4〜7はコンバータBに接続
• サンプリングするチャネルを指定し，シーケンシャルまたはパラレルサンプリング
  • 1回のサンプリング命令で最大8サンプル（チャネル）をサンプリングしデータを格納
• 2チャネルの差分値，オフセットなどの設定が可能

　GetAd関数からは想像もつかない作りになっています（？）．この辺りを自由自在に操れると，SilentCでもまだまだ頑張れる気がしてきます．

　実際にこれを操作する方法も，リファレンスマニュアルに載っています．ここでは簡単に示します．

1. ADLST1とADLST2の設定(28.4.4/500ページ)．リザルトレジスタが8サンプル分用意されているので，どのリザルトレジスタにどのチャネルからのサンプリングデータを注ぎ込むか指定する．
   • パラレルサンプリング実行時にはチャネル指定に制限があるので注意する．
2. 必要であれば，ADOFSnでオフセットを設定(28.4.11/509ページ)．
3. 必要であれば，CTRL1でサンプリング方法を指定．
4. CTRL1を設定してサンプリングし，ADRSLTnを読み取ることの繰り返し．
   • ADRSLTnはシフトされていることに注意する．

　これを使って全チャネル(AN0〜7)を読み取ってみましょう．

adc(){
  int i;
  char *ctr; ctr=0x40190000; //ADC CTRL1
  int *dat; dat=0x40190012; //ADC ADRSLT
  InitAd(0x70);
  #stop 0
  for(;;) {
    if(Getc(0)!=0)break;
    ctr[0]=32; //サンプリング
    for(i=0;i<8;i++){
      PrNum(dat[i]>>3);PrStr(", "); //AN0-AN7
    }
  PrStr("\r\n");
  SystemSleep();
}

2548, 2352, 2256, 2211, 1921, 2238, 2771, 2275,
2550, 2335, 2158, 2138, 1923, 2238, 2789, 2195,
2544, 2306, 1972, 1990, 1924, 2250, 2802, 2272,
2548, 2304, 2029, 2181, 1930, 2249, 2801, 2470,
2551, 2305, 2288, 2388, 1910, 2210, 2790, 2306,

OK

　AN4,5,6がX,Y,Z軸です．ちゃんと読み取れました．サンプリングにはある程度時間が掛かりますが，SilentCがコードを読む間に終わっているようなのでウェイトは入れません．さて，これでどれだけ早くなるでしょうか．

loop_axis(){
  char t1, *k; int *a,i,cr;
  k=0x40190000; a=0x40190012;
  t1=CreateTimer(0,30000,0);
  InitAd(0x70);
  for(i=0;i<1000;i++){
    k[0]=32;
    cr=a[4] >> 3;cr=a[5] >> 3;cr=a[6] >> 3;
  }
  cr=GetTimerCount(t1);
  PrNum(cr);PrStr("\r\n");
  CreateTimer(t1,0,0);
}

Ax_Timer::loop_axis
29710

Ax_Timer::loop_axis
29705

　ビットシフトして値を出すところまで行いました．およそ3 s，GetAd関数を呼ぶ時よりも10倍早くなっています！　移動平均などの計算時間があるのでストレートに10倍のサンプリング周波数は難しいですが，GetAd関数での待ち時間がなくなるのは大きいです．

参考

• http://noritan-micon.blog.so-net.ne.jp/
• http://www.freescale.com/files/32bit/doc/ref_manual/MCF52235RM.pdf
• http://www.cqpub.co.jp/interface/download/2008/09/kairozu.pdf

_ [ColdFire] ColdFire×P2P地震情報(6) 計測震度との対応を見る

　Interface 9月号付録基板を使い，P2P地震情報のサービスを作ってしまう勝手な企画です．過去の記事は，タイトルのColdFireカテゴリからどうぞ．

前回までのまとめ

• 「ズレ」「ノイズ」「サンプリング間隔」の問題をクリアし，"地震計としてのColdFire基板"が出来上がりました．
• これから，いよいよ"ColdFire基板 地震計"を制作していきます。

計測震度との対応を見る

　揺れの大きさとして「震度」が一般的であることや，震度階級の元となる計測震度が揺れの周期などを十分考慮して決定されていることを考えれば，ColdFire基板でもおおよそ相当する震度を出したいものです．

　そこで，これまでズレやノイズを解消するために作った『フィルタ』を通した値が，気象庁の計測震度とどう関係しているかを見ていきます．

　幸い，気象庁には100 Hzサンプリングの加速度データと計測震度が一部公開されています．20 Hzサンプリングにするため5個に1個の間隔でデータを取得してフィルタを通し，ベクトルを用いて加速度の大きさを求めます．

　フィルタを通さない場合の3方向，通した場合の3方向と2方向，3ケースを比べました．フィルタを通すことにより，計測震度との対応が良くなっていることが分かります．±0.5くらい誤差がありますが，ひとますこれでやってみましょう．

まずはシミュレーション

　まずは，表計算ソフトを用いてシミュレーションしてみます．ExcelやOpenOffice.org Calcを例にすると，フィルタはAVERAGE，ベクトルは各軸の2乗を足したものをSQRTすればよいので，プログラムと全く同等のシミュレーションができます．

　震度の算出については，Excelが "y = 0.7968 Ln(x) + 0.9717" と言っているのでこれを用います．最初に，静止した状態のノイズレベルがどの程度の震度になるのかシミュレーションしてみました．

　なんということでしょう！　静止しているのに震度3だなんて…

　つづく…

参考

• http://www.seisvol.kishou.go.jp/eq/kyoshin/jishin/index.html

_ [ColdFire] ColdFire×P2P地震情報(5) ノイズを除去する

　Interface 9月号付録基板を使い，P2P地震情報のサービスを作ってしまう勝手な企画です．過去の記事は，タイトルのColdFireカテゴリからどうぞ．

(2)〜(4) のおさらい

• 地震計として使うためには，「ズレ」「ノイズ」「サンプリング間隔」の3点を考える必要があるようです．
  • タイマを使うことで，サンプリング回数を保つことがそれなりにできました．
  • "0 gの電圧差""傾いた設置"の問題は，静止した状態を"地震計にとっての0 g"と考えることでクリアしました．
  • 時間の経過による傾きの変化などは，移動平均を用いることでクリアしました．

1 gあたりの感度差は，今回はそのままで

　"(2) 地震計としてのColdFire基板"において，『加速度センサによる1 gあたりの感度差は約±75 Galの差として現れることがあります．』と記しました．この差を補正するには，それぞれの軸で0 gの状態と1 gの状態を調べ，正しい値となるよう常に一定倍率を掛ける方法が考えられます．

　しかし，一般的な環境で各軸に0 gと1 gの状態を作り出すのは難しく，正確な補正は期待できません．誤った補正が掛かる可能性もあります．さらに，この感度差による誤差はそれほど大きいものではありません．

　感度差が 740 mV/g 〜 800 mV/g (Typical) 〜 860 mV/g の間で変動するということは， 92.5 % 〜 107.5 % の間で変動する，つまり誤差は±7.5 %ということになります．

　気になると言えば気になりますが，そこまで問題となるレベルではありませんので今回はそのままにします．

"移動平均" でノイズ問題もクリア

　これまで様々な問題を取り上げてきましたが，最後に『ノイズがあります．』と記した部分を解決しましょう．前回オフセットの算出に移動平均を使いましたが，これをデータにも使います．「データをならす」といったイメージで，細かなノイズを取り除きます．

　これも一種のフィルタで，「非常にすばやい変化には追従できない」という副作用があります．以下に数パターンの特性を，気象庁が計測震度の計算に用いるフィルタの特性とともに示しました．凡例は "サンプリング周波数 / 移動平均サンプル数" を表しています．

　高周波数帯（1 Hz以上）の特性が変わり，ローパスフィルタのように作用しています．気象庁フィルタと同じとはいきませんが，1 Hz〜3 Hz辺りは"20 Hz / 4 Samples"が近い値になっています．

　実は，これだけではノイズは完全に取れていません．フーリエ変換すると分かりますが，静止状態のノイズはかなり幅広く分布しています．ただ，ノイズ除去を厳しくしすぎると地震波特有の波形を見失う恐れがあるので，これ以上フィルタを掛けるのは避けます．

　それでは恒例になりましたサンプルプログラム，前回のオフセットと今回のノイズ除去を組み合わせてみました．サンプリング間隔が未調整で25 Hzくらいで出てきますが，あくまでサンプルです．

main(){
  int i,index,ind_n;
  int *sample=MemoryAlloc(120); int *nr=MemoryAlloc(6); long *offset=MemoryAlloc(12);
  // init
  InitAd(0x70); index=0; ind_n=17;
  for(i=0;i<3;i++){offset[i]=0;nr[i]=0;}
  for(i=0;i<60;i++)sample[i]=0;
  // loop
  #stop 0
  for(;;) {
    if(Getc(0)!=0)break;
    index=(index+1)%20;ind_n=(ind_n+1)%20;
    for(i=0;i<3;i++){
      offset[i]-=sample[i*20+index];nr[i]-=sample[i*20+ind_n];
      sample[i*20+index]=GetAd(i+4);
      offset[i]+=sample[i*20+index];nr[i]+=sample[i*20+index];
      PrNum(nr[i]/3-offset[i]/20);PrStr(", ");
    }
    PrStr("\r\n");SystemSleep();
  }
  MemoryFree(sample);MemoryFree(nr);MemoryFree(offset);
}

2種類の"移動平均"で，気象庁フィルタに近づく？

　サンプリング周波数20 Hzを基本とし，オフセット（基準点）に最新20サンプルの移動平均を用いるケース，データに最新3サンプルの移動平均を用いるケース，それらを組み合わせたケース，そして気象庁が計測震度の計算に用いるフィルタ，それぞれの特性を示しました．

　当初目標としていた『5HzPGA』，0.1 Hz〜5 Hzまでのバンドパスフィルタを掛けることと比べると，0.1 Hzや5 Hzで信号強度が約0.3 倍に落ち込んでいることから，『5HzPGA』として用いるには不適切な特性であることが分かります．まっとうにバンドパスしようとすると複雑な処理が必要になるので，無理があります．

　一方，気象庁フィルタと比べると多くの部分が近似しています．気象庁の計測震度「っぽいもの」を目指す可能性が急浮上してきました．がんばってみましょう．

今回のまとめ

• ノイズについては，移動平均を使うことで緩和しました．
• 何故か気象庁フィルタに近い特性を示したので，その方向でがんばることにしました

次回以降の予定

　これで全問題をクリアし，"地震計としてのColdFire基板"が一応出来上がりました．次回以降は，以下のような点を考えて"ColdFire基板 地震計"を制作していきます．

• 震度の計算が必要ですが，気象庁の方法を実装するのは無理があります．また，『5HzPGA』と別物になるため，5HzPGAのデータに基づく算出も出来ません．
  • 計測震度と，その元になった波形データをフィルタに通した値との相関を見出せるかもしれません．
  • ここから最大加速度だけでがんばれるかもしれません．
• 通信仕様を全く考えていません．

参考

• http://www.seisvol.kishou.go.jp/eq/kyoshin/kaisetsu/calc_sindo.htm
• http://www.jma.go.jp/jma/kishou/know/shindo/keisoku.html

_ [ColdFire] ColdFire×P2P地震情報(4) オフセット（基準点）を求める

　Interface 9月号付録基板を使い，P2P地震情報のサービスを作ってしまう勝手な企画です．過去の記事は，タイトルのColdFireカテゴリからどうぞ．

前回・前々回のおさらい

• 地震計として使うためには，「ズレ」「ノイズ」「サンプリング間隔」の3点を考える必要があるようです．
  • タイマを使うことで，サンプリング回数を保つことがそれなりにできました．

"設置の傾き"と"0 gの電圧差"の問題をまとめてクリアする

　前々回，『加速度センサによる0 gの電圧差は約±200 Galの差として現れることがあります．』と記しました．これは，オフセットを求めればすぐ解決します．電圧差が短時間で変化しないことが前提になりますが，「ここが0 gなんだよ」と教えてやれば良いのです．

　しかし，もうちょっとひねりましょう．基板が傾いて設置されているかもしれませんし，時間の経過とともに傾きが変わるかもしれません．さらに，求めたいのは「重力を取り除いた地震による加速度」です．

　つまり，傾きがどうであれ「"各軸に重力が掛かって静止している状態"が"地震計にとっての0 g"」，求めるべきは静止状態のオフセット，ということです．そうすれば，0 gの電圧差や傾いた設置の問題をまとめてクリアできます．

静止状態のオフセット

　静止状態のオフセットを出来るだけ正確に求めるには，静止したまま何サンプルか取って平均値を求めるのが良いでしょう．

calibration() {
  int i; long *avg;
  avg = MemoryAlloc(12);
  avg[0]=0;avg[1]=0;avg[2]=0;
  InitAd(0x70);
  for (i = 0; i < 30; i++) {
    avg[0] += GetAd(4); avg[1] += GetAd(5); avg[2] += GetAd(6);
    Sleep(10);
  }

  for (i = 0; i < 3; i++) {
    avg[i] = avg[i] / 30;
    PrNum(avg[i]); PrStr(",");
  }
  MemoryFree(avg);
}

2018,2058,2825,
OK

　ポインタ配列が登場していますが，その話はここではしません．30サンプルを足して30で割ることで平均値を求めました．地震による加速度の変化を調べるには，"最新サンプルの値 - 平均値"を計算することで実現できます。

　しかし，時間の経過によって基板の傾きが変わったり，地震によって傾いたまま戻らないケースも考えられます．これに対応するには、オフセットとなる平均値をサンプリング毎に再計算すると良いでしょう．

main(){
  int i,index;
  int *sample=MemoryAlloc(60); long *offset=MemoryAlloc(12);
  // init
  InitAd(0x70); index=0;
  for(i=0;i<3;i++)offset[i]=0;
  for(i=0;i<30;i++)sample[i]=0;
  // loop
  #stop 0
  for(;;) {
    if(Getc(0)!=0)break;
    index=(index+1)%10;
    for(i=0;i<3;i++){
      offset[i]-=sample[i*10+index];
      sample[i*10+index]=GetAd(i+4);
      offset[i]+=sample[i*10+index];
      PrNum(sample[i*10+index]-offset[i]/10);PrStr("(avg:");
      PrNum(offset[i]/10);PrStr("), ");
    }
    PrStr("\r\n");
  }
  MemoryFree(sample);MemoryFree(offset);
}

4(avg:2006), 19(avg:2072), 8(avg:2812),
0(avg:2006), 18(avg:2075), -17(avg:2811),
6(avg:2007), 5(avg:2075), -1(avg:2808),
-2(avg:2007), -18(avg:2076), 9(avg:2810),
7(avg:2008), 4(avg:2076), 8(avg:2810),

OK

　例として，最新10サンプルの平均値をオフセットにしました．シェイクすればきちんと変動しますし，傾きを変えても適応していきます．このような平均値の求め方を『移動平均』と呼びます．

　移動平均をオフセットにすると，副作用として「非常にゆっくりとした加速度の変化は，オフセットが追従するために値として現れなくなる」ということがあります．一種のフィルタとして作用するわけです．

　数パターンの特性と，気象庁が計測震度の計算に用いているフィルタの特性を上記に示します．凡例は "サンプリング周波数 / オフセットに使う移動平均サンプル数" を表しています．

　これを見ると，"10Hz / 10Samples"や"20Hz / 20Samples"の低周波数帯（1Hz以下）の特性が気象庁フィルタの特性とほぼ一致していることが分かります．なかなか面白いことになってきました．

　特性はノイズ対策の辺りに回すとして，これで"0 gの電圧差"問題は解決としましょう．

今回のまとめ

• "0 gの電圧差""傾いた設置"の問題は，静止した状態を"地震計にとっての0 g"と考えることでクリアしました．
• 時間の経過による傾きの変化などは，移動平均を用いることでクリアしました．

参考

• http://www.seisvol.kishou.go.jp/eq/kyoshin/kaisetsu/calc_sindo.htm
• http://www.jma.go.jp/jma/kishou/know/shindo/keisoku.html