信号包絡線の計算 / Computing signal envelope

信号包絡線（エンベロープ）は、信号の上部だけに接するような曲線で、ピーク検出などに使用される。

たとえば下図のような黒い実線の信号があるとき、その信号のおおよそのレベル変動を知りたいとする。音響信号は正負どちらの方向にも変動するので、その絶対値をとると灰色の破線のようになる。そこでできた頂点をなめらかに結ぶような曲線が包絡線である。

plot of chunk unnamed-chunk-2

信号包絡線を計算するには、ヒルベルト変換を使って解析信号を作成する必要がある。少し探してみたが、Rには標準でヒルベルト変換が用意されていないようだったので、まず以下のように関数を準備した。

## ヒルベルト変換で解析信号を計算する（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)))
}

使用法としては、以下のようになる。

## WAVファイルを扱うパッケージを読み込む
library(tuneR)

## WAVファイルを読み込み、途中の一部分だけを取り出す
wobj <- readWave("file_in.wav")
fs <- wobj@samp.rate
smpl <- seq(from=fs, to=fs+0.03*fs)
x <- wobj@left[smpl]
tm <- smpl/fs

## 包絡線の計算
xe <- envelope(x)

## 原信号と包絡線を重ねてプロットする
plot(tm, x, type="l", col=1, xlab="Time (s)", ylab="Amplitude", ylim=c(-20000,+20000))
par(new=T)
plot(tm, xe, type="l", col=2, xlab="", ylab="", ylim=c(-20000,+20000))

plot of chunk unnamed-chunk-5

`seewave`の`hilbert()`を使う場合

（2017-06-24追記：R 3.4.0でインストールできるseewaveのバイナリパッケージがないため、seewaveを使った方法はエラーが出ています。バイナリパッケージが登場したら更新します。)

## パッケージを読み込む
library(tuneR)
library(seewave)

## 
## Attaching package: 'seewave'

## The following object is masked _by_ '.GlobalEnv':
## 
##     hilbert

## WAVファイルを読み込み、途中の一部分だけを取り出す
wobj <- readWave("file_in.wav")
wobj2 <- extractWave(from=1, to=1.03, xunit="Time")

## Error in is(object, "Wave"): argument "object" is missing, with no default

## 包絡線の計算（ここではseewaveのhilbertを明示的に呼び出している）
xe <- abs(seewave::hilbert(wobj2))

## Error in is.data.frame(wave): object 'wobj2' not found

## 包絡線をプロットする
plot(xe, type="l")

plot of chunk unnamed-chunk-6

[←戻る](index.html)

信号包絡線の計算 / Computing signal envelope

seewaveのhilbert()を使う場合

`seewave`の`hilbert()`を使う場合