Empirical Mode Decomposition（EMD、経験的モード分解）は、時系列データをヒルベルト変換で扱いやすい振動成分に分解する方法です。ヒルベルト変換については，以下の記事も参考にしてください．
【Rで時系列解析】ヒルベルト変換で振動のエンベロープを抽出 - ケィオスの時系列解析メモランダム
EMDの前段階の知識として、ヒルベルト変換に関する以下の記事にある動画を見ておくと、イメージがわくと思います。
【Rで時系列解析】ヒルベルト変換でエンベロープを計算するときの注意 - ケィオスの時系列解析メモランダム
今回はEMDの計算過程を一通り説明します。
EMDの計算手順
EMDの計算手順は、以下です。
1. 極大値と極小値の点を見つける
下の動画の赤点と青点ように、まず、入力された時系列データの中で、極大値と極小値をすべて見つけます。
Rで処理する場合は、以下ような関数を使います。
find_extrema <- function(x) { n <- length(x) if(n < 3) return(list(max_idx = integer(0), min_idx = integer(0))) # 内部の要素に対して、前後との比較を一括で実施 max_idx <- which((x[2:(n-1)] > x[1:(n-2)]) & (x[2:(n-1)] > x[3:n])) + 1 min_idx <- which((x[2:(n-1)] < x[1:(n-2)]) & (x[2:(n-1)] < x[3:n])) + 1 list(max_idx = max_idx, min_idx = min_idx) }
この関数を定義して、時系列のベクトルxに対し、
i.extrema <- find_extrema(x)
とすれば、極大値の配列番号がi.extrema$max_idx、極小値の配列番号がi.extrema$minax_idxに入ってます。
2. 上側の包絡線と下側の包絡線を求める
見つけた極大値をなめらかな曲線でつないで上側の包絡線（エンベロープ）を作ります。包絡線は、スプライン補間したり、線形補間したりして計算します。 極大値の包絡線と同様に、極小値をつないで下側の包絡線も作ります。
ここまでの処理をRで実行する例が以下のスクリプトです。
# サンプルデータの作成 N <- 501 # データ長 i <- 0:(N - 1) # インデックス（0からN-1まで） # 複数の正弦波を重ね合わせた信号の生成 x <- sin(2 * pi * i / 200) + 0.9 * sin(5 * pi * (i - 11) / 200) + 0.8 * sin(16 * pi * i / 200) # 局所的な極大値と極小値のインデックスを求める関数 find_extrema <- function(x) { n <- length(x) # 要素数が3未満の場合は極値は存在しないので、空のリストを返す if (n < 3) { return(list(max_idx = integer(0), min_idx = integer(0))) } # ベクトル演算による局所的な最大値の検出 # x[2:(n-1)] が前後の要素 x[1:(n-2)] と x[3:n] より大きい場合 max_idx <- which((x[2:(n - 1)] > x[1:(n - 2)]) & (x[2:(n - 1)] > x[3:n])) + 1 # ベクトル演算による局所的な最小値の検出 # x[2:(n-1)] が前後の要素 x[1:(n-2)] と x[3:n] より小さい場合 min_idx <- which((x[2:(n - 1)] < x[1:(n - 2)]) & (x[2:(n - 1)] < x[3:n])) + 1 return(list(max_idx = max_idx, min_idx = min_idx)) } # 局所的な極値のインデックスの抽出 i.extrema <- find_extrema(x) # 局所的な極大値と極小値のインデックス（信号全体での位置） i.max <- i[i.extrema$max_idx] i.min <- i[i.extrema$min_idx] # 極値の個数を取得 n.max <- length(i.max) n.min <- length(i.min) # 局所的な極大値と極小値における信号の値 x.max <- x[i.extrema$max_idx] x.min <- x[i.extrema$min_idx] # 信号の端点を極値リストに追加（両端を含める処理） # 先頭の値の処理 if (i.min[1] > i.max[1]) { # 最初に極大値が現れる場合：最初の極小値の前に端点を追加 i.min <- c(0, i.min) x.min <- c(x[1], x.min) i.max <- c(0, i.max) x.max <- c(x.max[1], x.max) } else { # 最初に極小値が現れる場合：最初の極大値の前に端点を追加 i.min <- c(0, i.min) x.min <- c(x.min[1], x.min) i.max <- c(0, i.max) x.max <- c(x[1], x.max) } # 終端の値の処理 if (i.min[n.min] < i.max[n.max]) { # 最後に極大値が現れる場合：最後の極小値の後に端点を追加 i.min <- c(i.min, N - 1) x.min <- c(x.min, x[N]) i.max <- c(i.max, N - 1) # x.max の最後の値を繰り返し利用（※注意：この実装は意図に依存） x.max <- c(x.max, x.max[n.max + 1]) } else { # 最後に極小値が現れる場合：最後の極大値の後に端点を追加 i.max <- c(i.max, N - 1) x.max <- c(x.max, x[N]) i.min <- c(i.min, N - 1) # x.min の最後の値を繰り返し利用（※注意：この実装は意図に依存） x.min <- c(x.min, x.min[n.min + 1]) } # 上側包絡線（upper envelope）の作成 # 極大値の点を基に自然スプライン補間を行う upper_envelope <- spline(i.max, x.max, xout = i, method = "natural")$y # 下側包絡線（lower envelope）の作成 # 極小値の点を基に自然スプライン補間を行う lower_envelope <- spline(i.min, x.min, xout = i, method = "natural")$y # 上下包絡線の平均を求め、信号のトレンドなどの推定に利用可能 mean_envelope <- (upper_envelope + lower_envelope) / 2 # プロットの設定（n.plot が未定義の場合は使用箇所をコメントアウトまたは定義する必要があります） # i.mean_envelope.plot <- round(seq(1, N, length.out = n.plot)) # 信号と包絡線のプロット plot(i, x, type = "l", col = 1, lwd = 2, ylab = "x[i]", xaxs = "i", xlab = "i") lines(i, upper_envelope, col = 2, lwd = 2) lines(i, lower_envelope, col = 4, lwd = 2)
3. 平均包絡線を計算して、元の信号から引く
上側の包絡線と下側の包絡線のちょうど真ん中の線を計算します。この線を、平均包絡線と呼びます。 もともとの信号から、この平均包絡線を引き算します。
4. 上の処理をしつこく繰り返してIMFを求める
上の1～3の手続きを、何度も繰り返して、平均包絡線がほとんど真っ平になるようにします。
そうして得られるのが、最初のIntrinsic Mode Function (IMF)と呼ばれる振動成分です。
5. 残りの信号に対してさらに繰り返し処理
IMFを取り出した後の残りのデータに対して、上の1～5の処理を繰り返します。
これを何度も繰り返して、すべての成分を分解するまで進めます。
Rスクリプトの例を以下にのせておきます。
# サンプル信号の作成 N <- 1001 t <- 1:N x0 <- sin(2 * pi * (t-1)/200) + 0.9 * sin(5 * pi * ((t-1)-11)/200) + 0.8 * sin(16 * pi * (t-1)/200) # ベクトル演算による局所極値検出関数 find_extrema <- function(x) { dx <- diff(x) sdx <- sign(dx) dsdx <- diff(sdx) max_idx <- which(dsdx < 0) + 1 # 上昇→下降: 局所最大 min_idx <- which(dsdx > 0) + 1 # 下降→上昇: 局所最小 return(list(max_idx = max_idx, min_idx = min_idx)) } # EMDの実行：各IMFを順次抽出 imfs <- list() # 抽出したIMFの保存用リスト residual <- x0 # 初期残差は元の信号 imf_index <- 1 repeat { extrema <- find_extrema(residual) if (length(extrema$max_idx) + length(extrema$min_idx) < 2) break h <- residual SD <- 1 # sifting停止条件用 iter <- 0 max_iter <- 100 while (SD > 0.3 && iter < max_iter) { iter <- iter + 1 extrema <- find_extrema(h) if (length(extrema$max_idx) < 4 || length(extrema$min_idx) < 4) break # --- 端点処理 --- i.max_orig <- extrema$max_idx i.min_orig <- extrema$min_idx if (i.min_orig[1] > i.max_orig[1]) { i.min <- c(1, i.min_orig) x.min <- c(h[1], h[i.min_orig]) i.max <- c(1, i.max_orig) x.max <- c(h[i.max_orig[1]], h[i.max_orig]) } else { i.min <- c(1, i.min_orig) x.min <- c(h[i.min_orig[1]], h[i.min_orig]) i.max <- c(1, i.max_orig) x.max <- c(h[1], h[i.max_orig]) } if (i.min_orig[length(i.min_orig)] < i.max_orig[length(i.max_orig)]) { i.min <- c(i.min, length(h)) x.min <- c(x.min, h[length(h)]) i.max <- c(i.max, length(h)) x.max <- c(x.max, x.max[length(i.max_orig) + 1]) } else { i.max <- c(i.max, length(h)) x.max <- c(x.max, h[length(h)]) i.min <- c(i.min, length(h)) x.min <- c(x.min, x.min[length(i.min_orig) + 1]) } # --- 端点処理ここまで --- indices <- 1:length(h) upper_env <- spline(i.max, x.max, xout = indices, method = "natural")$y lower_env <- spline(i.min, x.min, xout = indices, method = "natural")$y m <- (upper_env + lower_env) / 2 h1 <- h - m SD <- sum((h - h1)^2) / sum(h^2) h <- h1 } imfs[[imf_index]] <- h imf_index <- imf_index + 1 residual <- residual - h } # 結果のプロット n.imfs <- length(imfs) num_plots <- n.imfs + 1 # 元信号、各IMF、残差 par(mfrow = c(num_plots,1), mar = c(2, 4, 2, 2)) plot(x0, type = "l", main = "Original Signal", xlab = "Index", ylab = "",col=2,lwd=2) for (j in 1:(n.imfs-1)){ plot(imfs[[j]], type = "l", main = paste("IMF", j), xlab = "Index", ylab = "",col=4,lwd=2) } residual <- residual + imfs[[n.imfs]] plot(residual, type = "l", main = "Residual", xlab = "Index", ylab = "",col=4,lwd=2)
端点の処理について
端点の処理については、ヒルベルト変換の結果の安定性を優先するか、元信号の振動を検出することを優先するかで、異なるアプローチが考えられます。今回は，後者を想定しました。 上のスクリプトの端点処理は、私が適当に考えたものです。みなさんは、パッケージを使うか、生成AI先生に聞くかして、適当に端点の処理をしてください。
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。
それと、下の はてなブログランキング はあまり信用できる指標ではなさそうですが（笑）、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
