AWK で十分ですよ - その4 - 移動平均

2012年 12月 1日掲載

「AWK，それは~~最後の~~今でもフロンティア」

AWK でいろいろする第4回である．今回は，移動平均によるデータの平滑化を行ってみる．

まずは，サンプルデータ作りである．振幅 2 のサインカーブのデータ data0.txt を作って，これに乱数で誤差を加え data1.txt を作る．

$ awk 'BEGIN{for(x=0;x<=6.0;x+=0.25)print 2*sin(x)}' > data0.txt
$ awk '{print $1+rand()-0.5}' data0.txt > data1.txt

さて，できあがったデータは，どんな風になっているだろうか．グラフを描いて見てみたい．

ということで，簡単なグラフを描く AWKスクリプトを作ってみた．

作成する画像ファイルは，PBMフォーマットとした．このフォーマットはアスキーファイルなので AWK でもなんとかできる． 1次元のデータを読みこんで， 320x240 のサイズを目いっぱい使って点をプロットする AWKスクリプトが，以下の plot1pbm.awk である．

BEGIN{
  width = 320
  height = 240

  for(i = 0; i < width * height; i++){
    p[i] = 0;
  }
}
{
  d[NR] = $1
}
END{
  if (NR <= 1){
    printf "P1\n1 1\n0\n"
    exit
  }

  dmin = dmax = d[1]

  for(i in d){
    if(d[i] < dmin) dmin = d[i]
    if(d[i] > dmax) dmax = d[i]
  }

  if(dmin == dmax) dmax = dmin + 1

  for(i in d){
    ix = int((width - 1) * (i - 1) / (NR - 1))
    iy = int((height - 1) * (d[i] - dmin) / (dmax - dmin))

    p[ix + (height - 1 - iy) * width] = 1
  }

  print "P1"
  print width, height

  for(i = 0; i < width * height; i++){
    printf "%d ", p[i]
  }

  printf "\n"
}

さっそく，data1.txt を図化してみる．

$ awk -f plot1pbm.awk data1.txt > data1.pbm

data1.pbm

う〜ん，よく見るとたしかに点はプロットされてるけど，やっぱり見にくいなあ．ちょっとがんばって，データを線で結んでみよう（線分の描画には，ブレゼンハムのアルゴリズムを使う）．

BEGIN{
  width = 320
  height = 240

  for(i = 0; i < width * height; i++){
    p[i] = 0;
  }
}
{
  d[NR] = $1
}
END{
  if (NR <= 1){
    printf "P1\n1 1\n0\n"
    exit
  }

  dmin = dmax = d[1]

  for(i in d){
    if(d[i] < dmin) dmin = d[i]
    if(d[i] > dmax) dmax = d[i]
  }

  if(dmin == dmax) dmax = dmin + 1

  ix1 = 0
  iy1 = int((height - 1) * (d[1] - dmin) / (dmax - dmin))

  for(i = 2; i <= NR; i++){
    ix2 = int((width - 1) * (i - 1) / (NR - 1))
    iy2 = int((height - 1) * (d[i] - dmin) / (dmax - dmin))

    line(ix1, iy1, ix2, iy2)

    ix1 = ix2
    iy1 = iy2
  }

  print "P1"
  print width, height

  for(i = 0; i < width * height; i++){
    printf "%d ", p[i]
  }

  printf "\n"
}

function line(x1, y1, x2, y2)
{
  dx = x2 - x1
  dy = y2 - y1

  sx = 1
  sy = 1

  if (dx < 0){
    dx *= -1
    sx = -1
  }

  if (dy < 0){
    dy *= -1
    sy = -1
  }

  p[x1 + (height - 1 - y1) * width] = 1
  
  if (dx > dy){
    de = dx / 2
    for(j = dx; j > 0; j--){
      x1 += sx
      de += dy
      if (de > dx){
        de -= dx
        y1 += sy
      }
      p[x1 + (height - 1 - y1) * width] = 1
    }
  }else{
    de = dy / 2
    for(j = dy; j > 0; j--){
      y1 += sy
      de += dx
      if (de > dy){
        de -= dy
        x1 += sx
      }
      p[x1 + (height - 1 - y1) * width] = 1
    }
  }
}

$ awk -f plot1pbm.awk data1.txt > data1a.pbm

data1a.pbm

うん，これなら見やすくて，いいんじゃないかな．スクリプトは，ごちゃごちゃして読みづらくなっちゃったけど（なんか BASIC時代のプログラムみたい…）．

さて，本題に入ろう．上のばらつきのあるデータ data1.txt を，移動平均を使って平滑化してみる．

N個のデータ y₁, y₂, ... ,y_N を，隣接する 3点を使って移動平均した値 y~ は，

y~_j = ( y_j-1 + y_j + y_j+1 ) / 3
ただし，j = 2, 3, ... , N-1

である．これくらい，AWK で計算するのは簡単である．ささっと書いた AWKスクリプト moving_average.awk は以下のとおり．

NR==1{
  y1 = $1
}
NR==2{
  y2 = $1
}
NR>2{
  y3 = $1

  print (y1 + y2 + y3) / 3

  y1 = y2
  y2 = y3
}

これに，data1.txt をつっこむと，

$ awk -f moving_average.awk data1.txt > data2.txt
$ awk -f plot1pbm.awk data2.txt > data2.pbm

data2.pbm

もとの data1.txt と比べると，移動平均した data2.txt は，なめらかになっているのがわかる．

が，このスクリプトには大きな問題がある．データの数が減ってしまうのだ．

$ wc -l data1.txt data2.txt
      25 data1.txt
      23 data2.txt
      48 total

data1.txt では 25個あったデータが，data2.txt では 23個に減っている．これは，y~_j = ( y_j-1 + y_j + y_j+1 ) / 3 では，両端の y~₁ と y~_N を計算することができないからである．この，データ列の両端の値を求める方法は，次回に解説することにする．

今回の AWKスクリプト：

過去の記事：

[ KDO ホームにもどる]