前回の授業で説明したような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点で移動平均

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

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

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

一方、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
         + x[t-2] * 1 ; // 数個の移動平均だし、
   y[t] = s / (3+2+1) ; //  ループを使わずに書いてみる。
}

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

指数移動平均

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

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

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

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 ;
}

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

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

• 2019-07-04-mov-avg.xlsx

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

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

一方、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
         + x[t-2] * 1 ; // 数個の移動平均だし、
   y[t] = s / (3+2+1) ; //  ループを使わずに書いてみる。
}

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

指数移動平均

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

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

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

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 ;
}