信号の包絡線を計算する

統計以外の様々なことにもRを使用し始めています。今日はインパルス応答のピーク検出をするために、信号の包絡線（エンベロープ、envelope）を計算するプログラムを書きました。以前Matlabで書いたものをそのまま移植しようとしたら、ヒルベルト変換の関数がなかったので、そこから作っています。

## ヒルベルト変換で解析信号を計算する (Hilbert Transform)
hilbert <- function(x) {
  # 信号長が奇数の時には末尾にゼロを加えて計算
  if (length(x) %% 2 == 1) {
    is.odd.flag <- TRUE
    x <- c(x, 0)
  } else {
    is.odd.flag <- FALSE
  }

  # ヒルベルト変換はフーリエ変換→負領域をゼロに→逆変換で計算
  X <- fft(x, inverse=FALSE)
  H <- c(1, rep(2, length.out=length(X)/2-1),
         1, rep(0, length.out=length(X)/2-1))
  Y <- X * H
  y <- fft(Y, inverse=TRUE) / length(Y)

  # 信号長が奇数の時には末尾につけたゼロの分を削除
  if (is.odd.flag == TRUE) {
    return(y[1:(length(y)-1)])
  } else {
    return(y)
  }
}


## 信号包絡線 (Signal Envelope)
envelope <- function(x) {
  return(abs(hilbert(x)))
}

テストしてみると、以下のようになります。

> x <- c(1, 2, 3, 4)
> hilbert(x)
[1] 1+1i 2-1i 3-1i 4+1i

> envelope(x)
[1] 1.414214 2.236068 3.162278 4.123106

> t <- seq(from=-2*pi, to=3*pi, length.out=1000)
> x <- sin(t)
> y <- envelope(x)
> plot(t, x, type="l", ylim=c(-1, +1.5))
> lines(t, y, col="red")
> quartz.save(file="~/20131104_envelope.png", type="png")