統計以外の様々なことにも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")