R初心者がRMSをプロットする。 - Seeva's Hot Dog ScrapBox!

初投稿です。投稿テーマは「GNU RでRMSをプロットする」です。

初投稿がこんな一部の人向けの様な投稿になっていいのか？という感じですね。筆者はどこにでもいる芸術大学の大学生で、Rはこの4月に始めた本当の初心者です。

以前は同じGNUプロジェクトの中のOctaveと呼ばれる言語を音響処理用として少しだけ学んでいました。しかし、この春から統計心理にも少し足を踏み入れることになり、同じソフトウェア上で統計と音響処理が出来ればと思い、Rに挑戦することにしました。

さて、そもそもRって何？って方も多いと思います。

RというのはGNUという、Unix向けの自由ソフトウェア運動の名の下にソフトウェアを配布している団体によるソフトウェアの１つです。GNUは、自由に使用・改変・再配布等ができるフリーウェアを作り、日頃使うすべてのソフトウェアをそれらでまかなうのを目標としている団体です。

そして、肝心のRですが、そんな彼らが作るフリーウェアの一つで統計処理などに多く用いられています。Rで使用できるR言語というのは、AT&T ベル研究所で作られたS言語、そしてプログラム言語Lispの方言に当たる言語で、MITの学生が作ったSchemeという言語などを参考に作られたものです。

統計処理に向いた言語ですが、パッケージと呼ばれるR言語で書かれたプログラムのライブラリを入れることで音響処理や画像処理なども行う事が出来ます。また、パッケージはCRANと呼ばれるネットワークによって共有されており、インターネットに接続できる環境があれば、世界中で書かれたあらゆるパッケージを読み込むことが出来ます。

~~今回はRで音の信号処理をしていきますが、Rでそのまま音を読み込むことは出来ないので、tuneRと呼ばれるパッケージを使用します。~~

~~source("tuneR")~~

~~（tuneRをお持ちでない方は、GNU Rソフトウェアのパッケージインストーラー等からインストールしておいてください。）~~
今回は音の入力などを行っていなかったので、これは必要無いですね・・・^^;

というわけで本題です。

RMSはRootMeanSquareのそのままの意味で、二乗して平均を取ったものの平方根を取ることです。
当然、時間軸方向にある程度の範囲が必要なので、今回は100msの窓を10分の1の大きさで移動させながら計算値を出していきます。FFTなどで用いられる窓がけの様な感じです。

プログラムの大まかな流れとしては、
1. 窓の大きさを計算する。
2. 窓に合う長さに音の長さを調節する。（ゼロ詰め）
3. 計算回数を設定する。
4. RMSを計算する。
5. 計算したRMS値を、最大値を0dBとした対数の値に変換する。
6. 計算された物をプロットする。
という流れになります。

plot.rms <- function(x,fs,nbits){
  
  #buffer time(ms)
  buffer <- 100
  buffer.samp <- ((buffer/1000) * fs)
  
  #add zero to fit 
  x.length1 <- length(x)
  x.length2 <- buffer.samp - (x.length1 %% buffer.samp)
  if(x.length2 == buffer.samp){
    x.length <- x.length1
  }else{
    x[x.length1 : (x.length1 + x.length2)] <- 0
    x.length <- x.length1 + x.length2 
  }
  
  
  #number of trials
  x.end <- (x.length / (buffer.samp / 10)) - ( 2 * buffer / 10 ) 
  
  x.rms <- c(1:x.end)
  x.rms[1:x.end] <- 0 
  
  
  k <- 1
  l <- 0
  while(k < x.end){
    
    x.rms[k] <-sqrt( sum( x[ ( l + 1 ) : l + buffer.samp ] ^ 2 )  / buffer.samp)
    
    k <- k + 1
    l <- l + buffer.samp / 10
    
  }
  
  #x.rms change to dB
  x.rms <- x.rms / max(x.rms)
  x.rms <- 20 * log10(x.rms)
    
  time.length <- c(1:x.end)
  time.length <- time.length / x.end * x.length / fs
  
  plot(time.length,x.rms,type="l",xlim=c(0,max(time.length)),ylim=c(-80,1),xlab="Time(sec)",ylab="PowerSpectrum(dB)",col="blue")
  par(new=T)
  abline(h=c(-20,-40,-60,-80),v=c(0:max(time.length)),lty="dotted")
  
  
}

ゼロ詰めの必要性ですが、もしゼロ詰めをせずに進めてしまうと、最後の方の試行で計算結果がNAになってしまい、プロット出来なくなってしまいます。
また、time.lengthが出てくる直前までは、xの値などはLarge numericとして扱われています。つまりサンプル値として扱われているので、冒頭のbufferをサンプル値に変換したり、最後のプロット時には、分かりやすいように秒に直すなどの手間が必要です。

というわけで、これを使えばあなたもRMSがグラフにプロット出来るはずです。