波形処理をハードウェアで行うかソフトウェアで行うか

組込み用の小型コンピュータが普及する以前は、このような波形に対する処理を行う場合には、電子回路的に波形のフィルタリングを行っていた。しかし電子回路的な方法では、回路の特性が変化してうまく処理ができなくなった場合に、回路の組み換えが発生するかもしれない。ただし回路の変更は基板の作り直しが必要であったりすることから、コストがかかる。

一方、A/D変換機能を内蔵した組込み用小型コンピュータも極めて安価になってきた。

• PIC(コントローラ)
• AVRマイコン / Arduino

こういったコンピュータの普及で、最近ではアナログ値をコンピュータに取り込んでデジタルの数値的な処理で余計な情報を取り除く。この方法であれば、プログラムを変更するだけなので、安価に変更が可能となる。ただし、アナログ値をA/D変換するのには時間がかかることから、周波数の高い信号には不向きである。

単純移動平均

前回説明を行った単純移動平均は、時刻ｔの平均を、その前後のデータで平均を求めた。この方式は、実際には与えられた波形のデータを全部記録した後に、単純移動平均をとる場合に有効である。

しかし、時々刻々変化する測定値の平均をその都度使うことを考えると、上記の方法は、未来の測定値を使っていることから、現実的ではない。

// 単純移動平均(未来の値も使う)
#define NS 3
int x[ SIZE ] ; // 入力値
int y[ SIZE ] ; // 出力値
for( int t = NS ; t < SIZE-NS ; t++ ) {
   int s = 0 ;
   for( int i = -NS ; i <= +NS ; i++ ) // 2*NS+1回の繰り返し
      s += x[t+i] ;
   y[t] = s / (2*NS + 1) ;
}

過去の値だけを使った移動平均

そこで、過去の値だけで移動平均をとることも考えられる。

この、単純移動平均と、過去の値だけを使う単純移動平均を、適当な測定値に対して適用した場合のグラフの変化を Excel によってシミュレーションした結果を以下に示す。

• 2018-06-19-mv-avg-excel

しかし、このグラフを見ると、波形後半の部分に注目するとよく分かるが、過去の値だけを使った移動平均では、測定値が立ち上がったのを追いかけて値が増えていく。これでは移動平均は時間的な遅れとなってしまう。

// 未来の値を使わない単純移動平均
for( int t = NS ; t < SIZE ; t++ ) {
   int s = 0 ;
   for( int i = 0 ; i <= NS ; i++ ) // NS+1回の繰り返し
      s += x[t-i] ;
   y[t] = s / (NS+1) ;
}

加重移動平均

過去の値を使った移動平均では遅れが発生する。でも、平均を取る際に、「n回前の値」と「現在の値」を考えた時、「その瞬間の平均値」は「現在の値」の方が近い値のはず。であれば、平均を取る時に、「ｎ回前の値は少なめ」「現在の値は多め」に比重をかけて加算する方法がある。

for( int t = 3 ; t < SIZE ; t++ ) {
   // 数個の移動平均だし、
   // ループを使わずに書いてみる。 
   int s = x[t]   * 3   // 現在の値は大きい重み
         + x[t-1] * 2   // 1つ前の値
         + x[t-2] * 1 ; // 2つ前の値(重みは最小)
   y[t] = s / (3+2+1) ;
}

この様に、過去に遡るにつれ、平均をとる比重を直線的に小さくしながら移動平均をとる方法は、加重移動平均と呼ばれる。以下にその変化をExcelでシミュレーションしたものを示す。

指数移動平均

ここまで説明してきた、単純移動平均や、加重移動平均は、平均をとる範囲の「過去の値」を記憶しておく必要がある。広い時間にわたる移動平均をとる場合は、それに応じてメモリも必要となる。これは、組み込み型の小型コンピュータであれば、メモリが足りず平均処理ができない場合もでてくる。

そこで、荷重移動平均の重みを、は、100%,は50%,は25%… というように、過去に遡るにつれ、半分にして平均をとる。

しかし、以降の項で、 を使うと以下のように書き換えることができる。

// 指数移動平均は、プログラムがシンプル
//  １つ前の平均y[t-1]を覚えるだけでいい。
for( int t = 1 ; t < SIZE ; t++ ) {
   y[t] = ( x[t] + y[t-1] ) / 2 ;
}

この方法であれば、直前の平均値を記録しておくだけで良い。このような移動平均を、指数移動平均と呼ぶ。

ここで示した指数移動平均は、過去を遡るにつれとなっているが、これをさらに一般化した指数移動平均は、以下の式で示される。前述の移動平均は、とみなすことができる。

#define ALPHA 0.5
for( int t = 1 ; t < SIZE ; t++ ) {
    y[t] = ALPHA * x[t] + (1.0 - ALPHA) * y[t-1] ;
}

以下のプログラムは、うまく動かない。理由を説明せよ。

#define RVA 4
for( int t = 1 ; t < SIZE ; t++ ) {
   // 以下はy[t]は全部ゼロになる。
   y[t] = 1/RVA * x[t] + (1.0 - 1/RVA) * y[t-1] ;

   // 以下は、整数型演算だけで、正しく動くだろう。
   // y[t] = ( x[t] + (RVA-1) * y[t-1] ) / RVA ;
}

理解度確認のための小レポート

上記の移動平均の理解のために、以下の資料(講義では印刷資料を配布)の表の中を、電卓などを使って計算せよ。
計算したら、その結果をグラフの中にプロットし、どういった波形となるか確認し、レポートとして提出すること。

• 2020-07-31-mov-avg

前回の授業で説明したようなA/D変換した数値データを読み取った場合、どのようなことが発生するか考える。

例えば、以下に示すような測定値があったとする。

• 2018-06-05-wave.csv

このデータの一部をグラフ化してみると、次のような波形であった。

この波形をみると、大きく見ればsinカーブだが、細かい点を見るとデータにブレがある。

誤差の原因

このような測定結果が得られた場合、本来コンピュータで処理したいデータは何であろうか？

原因は様々なものが考えられるが、

1. 回路のノイズ対策が不十分で、外部の電気的な影響が混入。
   オシロスコープで周期を図ると、60Hz なら、交流電源だったり…
2. D/A 変換を行う場合には、量子化誤差かもしれない。

例えば、最初の波形が、加速度センサーの値であったとして、船の上で揺れているために、大きな周期で加速度が変化しているかもしれない。一方で、船自体がエンジンによる揺れで加速度が変化しているかもしれない。

船の中で波の揺れと、エンジンの揺れが観測されている加速度センサーの情報で、船の揺れの大きさ・揺れの周期を知りたい場合、どうすればいいだろうか？

移動平均

このデータを見ると、10個のデータまでの間で、波形が上下に変動している。船の揺れとエンジンの揺れが原因であれば、10個ぐらいのデータのゆらぎが、エンジンによる揺れと考えられる。では、この10個ぐらいの範囲で値が上下の影響を減らしたければ、どうすればいいか？一番簡単な方法は、前後10個のデータで平均を取ればいいだろう。増減する値を加えれば、プラスの部分とマイナスの部分の値が相殺されて０に近くはず。そこでは、Excel で前後データの平均をとってみよう。

Excelで前後11点の平均を求める式をセルに入れる

青線：元波形データ(B列)、赤線：前後11点の平均(C列)

このように、データの前後の決められた範囲の平均を平均する処理は、移動平均(単純移動平均)と呼ぶ。

時間tにおけるデータをとした場合、前後5点の移動平均は、以下のような式で表せるだろう。

移動平均のプログラム

Excel で計算と同じ処理をプログラムで行うと以下のようになるだろう。

// moving-average.c
#include <stdio.h>

#define WIDTH 5
double data[ 1000 ] ; // 元データ
double ans[ 1000 ] ;  // 平均後のデータ

int main() {
   int t , i , size ;
   // 最初に全部のデータを読み込む
   for( size = 0 ; size < 1000 ; size++ ) {
      int num ;
      // コンマ区切りのデータを読む
      // 2つのデータが読み込めない時は入力を終了
      if ( scanf( "%d,%lf" , &num , &data[size] ) != 2 )
         break ;
   }
   // 移動平均を求める
   for( t = WIDTH ; t < size - WIDTH ; t++ ) {
      // t番目のデータの前後WIDTH個の合計
      double sum = 0 ;
      for( i = -WIDTH ; i <= WIDTH ; i++ )
         sum += data[ t + i ] ;
      ans[ t ] = sum / (2*WIDTH + 1) ;
   }
   // 計算後のデータをコンマ区切りで出力
   for( t = 0 ; t < size ; t++ ) {
      printf( "%d, %10.6lf, %10.6lf\n" , t , data[ t ] , ans[ t ] ) ;
   }
   return 0 ;
}

このプログラムを動かすと、データ番号とデータ値をコンマ区切りで与えること。

入力リダイレクトと出力リダイレクト

上記のプログラムでは、キーボードからデータを入力しなくてはいけない。これでは入力が大変なので、保存したファイルを使ってプログラムにデータを与える。

上記のプログラムを、パソコンの Z:¥課題¥moving-average.c に保存したとする。このプログラムを「コンパイル＆実行」すれば、Z:¥課題¥moving-average.exe という実行プログラムが作られ、プログラムが起動する。このままでは、キーボードからデータを入力する必要がある。

(1) ファイルから入力した値を使って処理を行うのであれば、コマンドを起動。

タスクバー左側の検索バーに、cmd.exe と入力すれば、命令入力画面が表示される。

(2) コマンド画面で、以下のように入力し、moving-average.exe があるか確認する。

  C:¥WINDOWS¥System32> Z:   青は表示される部分、赤が入力
  Z:¥> cd Z:¥課題
  Z:¥課題¥> dir *.exe
  06/21/2019 12:30PM  12345 moving-average.exe

(3) 最初のデータの記録されたCSVファイルを Z:¥課題 に保存する。

(4) コマンド画面で、以下のようにプログラム名の後ろに “＜ ファイル名” をつけて起動すると、キーボード入力の代わりに、ファイルから読み込んでプログラムが動く。このような起動は、入力リダイレクトと呼ぶ。

　Z:¥課題¥> moving-average.exe < 2018-06-05-wave.csv
　xx, xxx.xxxx, xxx.xxxx ←結果が画面に表示される

(5) これでは、結果がよく分からないので、ファイルに保存し Excel でグラフ化する。コマンド画面で、以下のようにプログラム名の後ろに“＞ ファイル名” をつけて起動すると、結果を画面に出力する代わりに、ファイルに結果を保存してくれる。このような起動は、出力リダイレクトと呼ぶ。

　Z:¥課題¥> moving-average.exe < 2018-06-05-wave.csv > out.csv

出力された out.csv は、データがコンマ区切りなので、Excel でひらけば、結果を表として簡単に読み込める。後は、グラフ化したい範囲を、マウスでドラッグ(もしくはシフトキーを押しながらカーソル移動)し、[挿入]-[グラフ]-[散布図]-[折れ線グラフ] でグラフ化すればいい。

自宅学習の課題

表計算ソフトで、移動平均を計算させてみよう。　　※

• 元波形
• 前後5点で移動平均
• 前後11点で移動平均
• 前後51点で移動平均

をとるような表計算の式を書き込んで、その結果の波形がどんなグラフになるのか確認しておくこと。

実数型(float / double)

実数型は、単精度実数(float型)と、倍精度実数(double型)があり、それぞれ32bit,64bitでデータを扱う。

単精度型(float)では、符号1bit,指数部8bit,仮数部23bitで値を覚え、数値としては、以下の値を意味する。(精度が低いので普通のコンピュータではあまり使われることはない)

符号✕ 1.仮数部 ✕ 2^{指数部-127}

倍精度型(double)では、符号1bit,指数部11bit,仮数部52bitで値を覚え、数値としては、以下の意味を持つ。

符号✕ 1.仮数部 ✕ 2^{指数部-1023}

これらの実数で計算を行うときには、0.00000001011×2¹⁰といった値の時に、仮数部に0が並んだ状態を覚えると、計算の精度が低くなるので、1.01100000000×2²のように指数部の値を調整して小数点の位置を補正しながら行われる。このため、float型であれば、23bit=10進数7桁相当、double型なら52bit=10進数16桁相当の精度がある。

倍精度型を使えば、正しく計算できるようになるかもしれないが、実数型はただの加算でも仮数部の小数点の位置を合わせたりする処理が必要で、浮動小数点専用の計算機能を持っていないような、ワンチップコンピュータでは整数型にくらべると10倍以上遅い場合もある。

実数の注意点

C言語でプログラムを作成していて、簡単な数値計算のプログラムでも動かないと悩んだことはないだろうか？解らなくて友達のプログラムを真似したら動いたけど、なぜ自分のプログラムは動かなかったのか深く考えたことはあるだろうか？

単純な合計と平均

整数を入力し、最後に合計と平均を出力するプログラムを以下に示す。
しかし、C言語でこのプログラムを動かすと、10,10,20,-1 と入力すると、合計(sum)40,件数(cnt)3で、平均は13と表示され、13.33333 とはならない。

小数点以下も正しく表示するには、どうすればいいだろうか？
ただし、変数の型宣言を “double data,sum,cnt ;” に変更しないものとする。

// 入力値の合計と平均を求める。
#include <stdio.h>

int main() {
   int data ;
   int sum = 0 ;
   int cnt = 0 ;
   for(;;) {
      printf( "数字を入力せよ。-1で終了¥n" ) ;
      scanf( "%d" , &data ) ;
      if ( data < 0 )
         break ;
      cnt = cnt + 1 ;
      sum = sum + data ;
   }
   printf( "合計 %d¥n" , sum ) ;
   printf( "平均 %d¥n" , sum / cnt ) ;
}

C言語では、int型のsum / int型のcnt の計算は、int 型で計算を行う(小数点以下は切り捨てられる)。このため、割り算だけ実数で行いたい場合は、以下のように書かないといけない。

   printf( "平均 %lf¥n" , (double)sum / (double)cnt ) ;
   // (double)式 は、sum を一時的に実数型にするための型キャスト

まずは動く例

以下のプログラムは、見れば判るけど、th を 0度〜360度まで5度刻みで変化させながら、y = sin(th) の値を表示するプログラム。

// sin の値を出力
#include <stdio.h>
#include <math.h>

int main() {
    double th , y ;
    for( th = 0.0 ; th <= 360.0 ; th += 5.0 ) {
        y = sin( th / 180.0 * 3.1415926535 ) ;
        printf( "%lf %lf¥n" , th , y ) ;
    }
    return 0 ;
}

動かないプログラム

では、以下のプログラムはどうだろうか？

// case-1 ---- プログラムが止まらない
#define PI 3.1415926535
int main() {
    double th , y ;
    // 0〜πまで100分割でsinを求める
    for( th = 0.0 ; th != PI ; th += PI / 100.0 ) {
        y = sin( th ) ;
        printf( "%lf %lf¥n" , th , y ) ;
    }
    return 0 ;
}

// case-2 ---- y の値が全てゼロ
int main() {
    int    th ;
    double y ;
    for( th = 0 ; th <= 360 ; th += 5 ) {
        y = sin( th / 180 * 3.1415926535 ) ;
        printf( "%d %lf¥n" , th , y ) ;
    }
    return 0 ;
}

どちらも、何気なく読んでいると、動かない理由が判らないと思う。そして、元のプログラムと見比べながら、case-1 では、「!=」を「<=」に書き換えたり、case-2 では、「int th ;」を「double th ;」に書き換えたら動き出す。

では何が悪かったのか…
回答編

---

表計算ソフトの使い方

情報制御基礎では、プログラムで計算する所を、Excel のような表計算ソフトを用いて検証してもらったりする予定なので、Excel で計算式を使う方法を説明する。

セルの場所と簡単な式

簡単な、品名・単価・個数・価格の表を考える。以下の表のように、列の名前と、品名・単価・個数まで入力した後、単価と個数をかけた価格を求めるとする。

Excel では、表の列には左から、A,B,C,D… , 表の行には上から1,2,3,4,5 と番号が振られていて、特定の列・特定の行のデータを表す時には、列行を組み合わせ、A1に品名、B3に￥80、C5に4 が入っている。

例えば、D2 に、ノート単価120円、ノート個数3個をかけた値を入れたい場合は、D2の場所に、

=B2*C2

を書き込めば、その場所には360が表示される。

Excelでは、入力する文字列の先頭が”=”の場合は、残り部分は計算式として扱われる。

D3には、”=B3*C3″を入力すれば、160 が表示される。しかし、この様な式を何度も入力するのは面倒である。

この場合、セル・カーソルを、D2 に合わせ、[右ボタン]-[コピー]を行い、D3 で[右ボタン]-[貼り付けオプション]-[貼り付け]を行えば、”=B3*C3″が入力される。

ここで注意しないといけないのが、式を張り付ける場合には、貼り付け先のセルの場所が一つ下の行なので、行番号を表す2の部分が1つ下の行番号3に書き換えられて、貼り付けが行われる。(相対参照)

関数式

例えば、下左図のような、数字とその平方根の表を作る場合、A2 に 1、B2に =sqrt( A2 ) を入力、A3 に =A2+1 を入力したあと、B2の式をB3にコピー＆ペーストし、A3,B3　を A4～A6にペーストすればいい。

B2に入力したような、sqrt( A2 ) のようなものは、関数式と呼ばれる。

また、A3,B3 といった複数の行・列をまとめた範囲を示す時は、A3:B3 といった表記方法であらわす。

絶対参照と相対参照

最初の例に戻って、単価と個数の積で今度は税率を加えて計算する例を考える。また、税率は後で変化するかもしれないので、B1 のセルに税率を記入しておく場合を考える。

この場合、D3 には、” =B3*C3*(1+B1) ” を入力すればいい。

ただ、このように式を入力すると、D3 の計算式を、D4,D5,D6 にコピーすると、セル D4 には =B4*C4*(1+B2) が入力されてしまい、B2 には単価という文字が記載されているため、正しい結果が求まらない。

こういった場合には、絶対参照を用いる。D3 に記入する式を

=B3*C3*(1+$B$2)

とし、この D3 の式を D4 にコピー＆ペーストすると、列記号、行番号の前に$がついた部分の式は、貼り付け場所に応じて変化しない。

このような、$B$2 といったセルの参照は、絶対参照と呼ぶ。これに対し、B2 といったセル参照は、貼り付け場所に応じて書き換えられるので、相対参照と呼ぶ。

絶対参照と相対参照が混ざった、$B2, B$2 といった書き方もある。
式の入力時にF4ボタンを押す度に、B2→$B$2→B$2→$B2→B2 と変化する

$B2 は、式をコピーすると列部分はBのまま、行部分は場所に合わせて変化する。

B$2 は、式をコピーすると列部分は場所に合わせて変化し、行部分は2のままとなる。

レポート課題(第4回)

Excel で、xを0〜180度まで変化させたときのsin(x),位相をyとした時のsin(x+y)の値の表を作り、グラフ機能で表示せよ。この時、計算式の入力をどのように行なったのか(相対参照や絶対参照をどのように使ったのか)説明を、グラフの下に入力欄を設け記入せよ。

そして出来上がった Excel のファイルを、指定されたフォルダに提出せよ。