【Rで瞬時周波数推定】Savitzky–Golay微分フィルタを用いたTeager–Kaiser Energy Methodの改善策

瞬時周波数や振幅を同時に推定できる TKEO（Teager–Kaiser Energy Operator） は、振幅変調・周波数変調（AM–FM）信号解析の古典的な手法です。しかし、実データへの適用を考えると、従来の TKEO 法には明確な問題があります。ごくわずかな加法ノイズが含まれるだけで，推定結果が急激に不安定になるのです。
　下図は、振幅変調・周波数変調信号に小さな加法ノイズを加えた例に対して、従来の TKEO を用いて瞬時周波数と振幅を推定した結果を示しています。

上段：観測された信号 $x(t)$
中段：真の瞬時周波数（青）と TKEO 推定値（赤）
下段：真の振幅（青）と TKEO 推定値（赤）

一般的な差分演算を用いたTKEO法の適用例。本記事の最後に掲載したRスクリプトの実行結果。解析対象の信号は0を中心に上下にほぼ対象に振動する必要があることに注意。

まず中段の瞬時周波数を見ると、本来は時間とともに滑らかに変化しているはずの周波数に対して、TKEO による推定値は全面的に激しくばらついており、真の変動構造をほとんど捉えられていません。ノイズは非常に小さいにもかかわらず、推定された瞬時周波数はほぼランダムに散乱してしまっています。
　同様の問題は、下段の振幅推定にも現れています。真の振幅変調は滑らかな正弦波状であるにもかかわらず、TKEO による推定結果はスパイク状の大きな変動を示し、振幅の時間変化として解釈することが困難な状態になっています。
　このように、

瞬時周波数が激しくばらつく
振幅推定がスパイク状になる

という現象は、従来の TKEO 法をノイズを含む実データに適用した際に典型的に見られる挙動です。その結果、理論的には有用なはずの TKEO であっても、実用的な瞬時周波数推定には使いにくい、という評価につながってきました。
　そこで今回は、Savitzky–Golay微分フィルタを用いて，従来の TKEO の欠点を改善してみます。

Savitzky–Golay 微分フィルタを使ったTKEO法（ケィオスのオリジナル）
■ Rスクリプト

Savitzky–Golay 微分フィルタを使ったTKEO法（ケィオスのオリジナル）

　Savitzky–Golay 微分フィルタ を用いたTKEO法のポイントは、TKEO の考え方そのものは一切変えず、「微分の推定方法」だけを置き換えるという点です。
　従来の TKEO では、信号の微分を単純な差分で近似していました。差分は計算が簡単である一方、高周波成分を強調する性質があり、わずかなノイズであっても推定結果に大きな影響を与えてしまいます。このことが、前節で見たような瞬時周波数推定の破綻や、振幅推定のスパイク状変動の主因となっています。
　そこで今回は、差分の代わりに Savitzky–Golay 微分フィルタ を用いて、信号の一次および二次微分を推定します。Savitzky–Golay 微分フィルタは、局所的に多項式近似を行った上で微分を評価する手法であり、平滑化と微分を同時に行えるという特徴を持っています。そのため、微分演算に伴うノイズの増幅を効果的に抑えることができます。
　この Savitzky–Golay 微分フィルタで得られた一次微分と二次微分を用いて、連続時間型の TKEO

$\displaystyle{ \Psi[x(t)] = \left(\frac{dx}{dt}\right)^2 - x(t)\frac{d^2x}{dt^2} }$

をそのまま離散化して計算します。重要なのは、TKEO のエネルギー演算や、比率を用いた瞬時周波数・振幅の推定手順自体は、従来法と全く同じであるという点です。変更しているのは、「微分をどう求めるか」だけです。
　下図は、先ほどと同じ振幅変調・周波数変調信号に対して、Savitzky–Golay 微分フィルタを用いた TKEO を適用した結果を示しています。従来の TKEO と比べると、瞬時周波数の推定値は真の変動を滑らかに追従しており、ノイズによるばらつきが大幅に抑えられていることが分かります。また、振幅推定についても、スパイク状の不安定な挙動はほとんど見られず、真の振幅変調と良好に一致しています。

この結果から分かるように、TKEO の実用性を損なっていた原因は、TKEO 自体ではなく、微分を差分で近似していた点にありました。Savitzky–Golay 微分フィルタを用いることで、TKEO の本来の能力を引き出し、ノイズを含む実データに対しても安定した瞬時周波数・振幅推定が可能になります。

■ Rスクリプト

Savitzky–Golay 微分フィルタを用いたTKEO法

　Savitzky–Golay 微分フィルタを用いたTKEO法をRで実装したものを掲載しておきます。このスクリプトでは、

w.SG   <- 11   # odd window length

の値を小さくすれば（w.SG <- 9 など）、高周波数成分の検出精度が向上し、大きくすれば（w.SG <- 15 など）、低い周波数成分の検出精度が向上します。この特性を利用し、目的の周波数領域の瞬時周波数・瞬時振幅の推定精度を向上することができます。　詳しい説明は、そのうち論文を書きます。

############################################################
# TKEO demo in R (Savitzky–Golay derivative-based TKEO)
#
# Continuous-time TKEO:
#   Psi{x(t)} = (dx/dt)^2 - x(t) * (d^2x/dt^2)
#
# Discrete implementation (SG-derivative w.r.t. sample index n):
#   x1[n] = d/dn x[n]  (estimated by SG deriv=1)
#   x2[n] = d^2/dn^2 x[n] (estimated by SG deriv=2)
#   Psi_SG[n] = (x1[n])^2 - x[n]*x2[n]
#
############################################################

# -----------------------------
# 0) Savitzky–Golay coefficients (general)
# -----------------------------
savgol_coeffs <- function(w, poly = 4L, deriv = 0L) {
  stopifnot(is.numeric(w), length(w) == 1L, w %% 2L == 1L, w >= 3L)
  stopifnot(is.numeric(poly), length(poly) == 1L, poly >= 0L)
  stopifnot(is.numeric(deriv), length(deriv) == 1L, deriv >= 0L, deriv <= poly)

  w <- as.integer(w); poly <- as.integer(poly); deriv <- as.integer(deriv)

  k <- (w - 1L) %/% 2L
  j <- (-k):k

  # Vandermonde design matrix: j^0..j^poly
  A <- outer(j, 0:poly, `^`)

  # derivative at 0: only a_deriv contributes with factor deriv!
  e <- rep(0, poly + 1L)
  e[deriv + 1L] <- factorial(deriv)

  ATA <- crossprod(A)
  sol <- qr.solve(ATA, e)
  as.numeric(A %*% sol)
}

# Centered SG convolution
sg_filter <- function(x, w, poly = 4L, deriv = 0L) {
  stopifnot(is.numeric(x))
  c <- savgol_coeffs(w = w, poly = poly, deriv = deriv)
  as.numeric(stats::filter(x, filter = c, sides = 2, method = "convolution"))
}

# -----------------------------
# 1) TKEO using SG derivatives (continuous-form discretization)
#   Psi_SG[n] = (x'[n])^2 - x[n]*x''[n]
# -----------------------------
tkeo_sg <- function(x, w, poly = 4L, return_derivs = FALSE) {
  stopifnot(is.numeric(x), length(x) >= 3L)
  stopifnot(w %% 2L == 1L, w >= 5L)   # recommend >=5 for stable 2nd derivative
  stopifnot(poly >= 2L)              # need at least quadratic for 2nd derivative

  x1 <- sg_filter(x, w = w, poly = poly, deriv = 1L)  # d/dn
  x2 <- sg_filter(x, w = w, poly = poly, deriv = 2L)  # d2/dn2
  psi <- x1^2 - x * x2

  if (return_derivs) return(list(psi = psi, x1 = x1, x2 = x2))
  psi
}

# -----------------------------
# 2) Synthetic AM–FM signal (noise-free-ish)
# -----------------------------
fs <- 200
dt <- 1 / fs
Tsec <- 10
n <- Tsec * fs
t <- (0:(n - 1)) * dt

# True amplitude (AM)
a_true <- 1.0 + 0.4 * sin(2 * pi * 0.2 * t)

# True instantaneous frequency (FM) [Hz]
f_true <- 8 + 3 * sin(2 * pi * 0.1 * t)

# Phase: phi(t) = 2*pi * integral f(t) dt
phi <- 2 * pi * cumsum(f_true) * dt

# Signal
x <- a_true * cos(phi) + rnorm(n, 0, 0.005)

# -----------------------------
# 3) SG-derivative TKEO estimation
# -----------------------------
w.SG   <- 11   # odd window length
poly.SG <- 4L # SG polynomial order

# Psi[x] via SG derivatives (continuous-form)
psi_x <- tkeo_sg(x, w = w.SG, poly = poly.SG)

# Also compute SG first derivative x' (w.r.t. sample index n) for the ratio method
x1 <- sg_filter(x, w = w.SG, poly = poly.SG, deriv = 1L)

# Psi[x'] via SG derivatives (continuous-form)
psi_x1 <- tkeo_sg(x1, w = w.SG, poly = poly.SG)

# Safety
eps <- 1e-12
psi_x_pos  <- pmax(psi_x, eps)
psi_x1_pos <- pmax(psi_x1, 0)

# Ratio for SG-derivative version:
# For x[n] = A cos(nΩ), if x1 approximates x' (d/dn),
#   Psi[x'] / Psi[x] ≈ Ω^2
r <- psi_x1_pos / psi_x_pos
r <- pmin(pmax(r, 0), pi^2)

Omega_hat <- sqrt(r)                 # rad/sample
f_hat <- (fs / (2 * pi)) * Omega_hat # Hz

# Amplitude estimate from Psi[x] = A^2 sin^2(Ω)
a_hat <- sqrt(psi_x_pos) / pmax(abs(sin(Omega_hat)), eps)

# -----------------------------
# 4) Plot (requested colors)
# -----------------------------
op <- par(no.readonly = TRUE)
par(mfrow = c(3, 1), mar = c(5, 5, 3, 1), las = 1)

# (Top) signal col=4 (show original only; no extra smoothing step requested)
plot(t, x, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "x(t)",
     main = sprintf("AM–FM signal (SG-derivative TKEO, w=%d)", w.SG))

# (Middle) frequency overlay: true col=4, estimated col=2
plot(t, f_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Frequency [Hz]",
     main = "Instantaneous frequency: true vs TKEO (SG-derivative, continuous-form)")
lines(t, f_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True f(t)", "TKEO (SG-derivative)"),
       col = c(4, 2), lwd = 2, bty = "n")

# (Bottom) amplitude overlay: true col=4, estimated col=2
plot(t, a_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Amplitude",
     main = "Amplitude: true vs TKEO (SG-derivative, continuous-form)")
lines(t, a_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True a(t)", "TKEO (SG-derivative)"),
       col = c(4, 2), lwd = 2, bty = "n")

par(op)

従来のTKEO法

############################################################
# TKEO demo in R
#  - Estimate instantaneous frequency f[n] and amplitude a[n]
#  - Discrete-time formulas consistent with Psi[x[n]]=A^2 sin^2(Omega)
############################################################

# -----------------------------
# 0) Helper: discrete TKEO
# -----------------------------
tkeo <- function(x) {
  n <- length(x)
  psi <- rep(NA_real_, n)
  if (n >= 3L) {
    psi[2:(n - 1L)] <- x[2:(n - 1L)]^2 - x[1:(n - 2L)] * x[3:n]
  }
  psi
}

# -----------------------------
# 1) Synthetic AM–FM signal (noise-free)
# -----------------------------
fs <- 200
dt <- 1 / fs
Tsec <- 10
n <- Tsec * fs
t <- (0:(n - 1)) * dt

# True amplitude (AM)
a_true <- 1.0 + 0.4 * sin(2 * pi * 0.2 * t)

# True instantaneous frequency (FM) [Hz]
f_true <- 8 + 3 * sin(2 * pi * 0.1 * t)

# Phase: phi(t) = 2*pi * integral f(t) dt
phi <- 2 * pi * cumsum(f_true) * dt

# Signal
x <- a_true * cos(phi) + rnorm(n,0,0.005)

# -----------------------------
# 2) TKEO-based estimation (discrete, correct formulas)
# -----------------------------
# Difference signal
y <- c(NA_real_, diff(x))   # y[n] = x[n] - x[n-1]

psi_x <- tkeo(x)
psi_y <- tkeo(y)

# Safety
eps <- 1e-12
psi_x_pos <- pmax(psi_x, eps)
psi_y_pos <- pmax(psi_y, 0)

# Ratio r = Psi[y]/Psi[x] = 4 sin^2(Omega/2)  (ideal sinusoid)
r <- psi_y_pos / psi_x_pos
r <- pmin(pmax(r, 0), 4)  # theoretical range [0,4]

# Omega_hat (rad/sample), internal only (not shown in the blog text if you prefer)
Omega_hat <- 2 * asin(0.5 * sqrt(r))

# Frequency estimate in Hz
f_hat <- (fs / (2 * pi)) * Omega_hat

# Amplitude estimate (consistent with Psi[x]=a^2 sin^2(Omega))
a_hat <- sqrt(psi_x_pos) / pmax(abs(sin(Omega_hat)), eps)

# -----------------------------
# 3) Plot (requested colors)
# -----------------------------
op <- par(no.readonly = TRUE)
par(mfrow = c(3, 1), mar = c(5, 5, 3, 1), las = 1)

# (Top) signal col=4
plot(t, x, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "x(t)",
     main = "AM–FM signal with small additive noise")

# (Middle) frequency overlay: true col=4, estimated col=2
plot(t, f_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Frequency [Hz]",
     main = "Instantaneous frequency: true vs TKEO")
lines(t, f_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True f(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

# (Bottom) amplitude overlay: true col=4, estimated col=2
plot(t, a_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Amplitude",
     main = "Amplitude: true vs TKEO estimate")
lines(t, a_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True a(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

par(op)

■ 従来法の改善のためのthree-pass Savitzky–Golay filter

　従来の離散 TKEO でも、前処理で平滑化処理すれば、結果は劇的に改善します。そこで、今回は私の必殺技 "three-pass Savitzky–Golay filter" で前処理します。これは、Savitzky–Golay filterを繰り返し３回かけるという単純な方法ですが、Savitzky–Golay filterの欠点を補うことができます。詳細は、以下の記事を読んでください。

chaosmemo.com

　詳しい説明は面倒なので、デモ用のRスクリプトのみ掲載しておきます。

############################################################
# TKEO demo in R
#  - Estimate instantaneous frequency f[n] and amplitude a[n]
#  - Discrete-time formulas consistent with Psi[x[n]]=A^2 sin^2(Omega)
#
# Preprocessing:
#  - three-pass Savitzky–Golay filter (w.SG = 9, poly = 3)
############################################################

# -----------------------------
# 0) Helper: discrete TKEO
# -----------------------------
tkeo <- function(x) {
  n <- length(x)
  psi <- rep(NA_real_, n)
  if (n >= 3L) {
    psi[2:(n - 1L)] <- x[2:(n - 1L)]^2 - x[1:(n - 2L)] * x[3:n]
  }
  psi
}

# -----------------------------
# 0.5) Preprocess: three-pass Savitzky–Golay filter (w=9, poly=3)
# -----------------------------
if (!requireNamespace("signal", quietly = TRUE)) {
  stop("Package 'signal' is required. Install it by install.packages('signal').")
}

w.SG   <- 9L   # window length (odd)
p.SG   <- 3L   # polynomial order
n_pass <- 3L   # three-pass

sg_smooth_3pass <- function(x, w = 9L, p = 2L) {
  stopifnot(is.numeric(x))
  stopifnot(w %% 2L == 1L, w >= 3L)
  stopifnot(p >= 0L, p < w)
  y <- x
  for (k in 1:3) {
    y <- signal::sgolayfilt(y, p = p, n = w)  # smoothing (deriv=0)
  }
  y
}

# -----------------------------
# 1) Synthetic AM–FM signal (noise-free)
# -----------------------------
fs <- 200
dt <- 1 / fs
Tsec <- 10
n <- Tsec * fs
t <- (0:(n - 1)) * dt

# True amplitude (AM)
a_true <- 1.0 + 0.4 * sin(2 * pi * 0.2 * t)

# True instantaneous frequency (FM) [Hz]
f_true <- 8 + 3 * sin(2 * pi * 0.1 * t)

# Phase: phi(t) = 2*pi * integral f(t) dt
phi <- 2 * pi * cumsum(f_true) * dt

# Signal
x <- a_true * cos(phi) + rnorm(n, 0, 0.005)

# --- Preprocessing: three-pass Savitzky–Golay filter (w=9, poly=3) ---
x <- sg_smooth_3pass(x, w = w.SG, p = p.SG)

# -----------------------------
# 2) TKEO-based estimation (discrete, correct formulas)
# -----------------------------
# Difference signal
y <- c(NA_real_, diff(x))   # y[n] = x[n] - x[n-1]

psi_x <- tkeo(x)
psi_y <- tkeo(y)

# Safety
eps <- 1e-12
psi_x_pos <- pmax(psi_x, eps)
psi_y_pos <- pmax(psi_y, 0)

# Ratio r = Psi[y]/Psi[x] = 4 sin^2(Omega/2)  (ideal sinusoid)
r <- psi_y_pos / psi_x_pos
r <- pmin(pmax(r, 0), 4)  # theoretical range [0,4]

# Omega_hat (rad/sample)
Omega_hat <- 2 * asin(0.5 * sqrt(r))

# Frequency estimate in Hz
f_hat <- (fs / (2 * pi)) * Omega_hat

# Amplitude estimate (consistent with Psi[x]=a^2 sin^2(Omega))
a_hat <- sqrt(psi_x_pos) / pmax(abs(sin(Omega_hat)), eps)

# -----------------------------
# 3) Plot (requested colors)
# -----------------------------
op <- par(no.readonly = TRUE)
par(mfrow = c(3, 1), mar = c(5, 5, 3, 1), las = 1,xaxs="i")

# (Top) signal col=4
plot(t, x, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "x(t)",
     main = "AM–FM signal with small additive noise (three-pass SG prefilter)")

# (Middle) frequency overlay: true col=4, estimated col=2
plot(t, f_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Frequency [Hz]",
     main = "Instantaneous frequency: true vs TKEO")
lines(t, f_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True f(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

# (Bottom) amplitude overlay: true col=4, estimated col=2
plot(t, a_true, type = "l", col = 4, lwd = 2,
     xlab = "Time [s]", ylab = "Amplitude",
     main = "Amplitude: true vs TKEO estimate")
lines(t, a_hat, col = 2, lwd = 3, lty = 2)
legend("topright",
       legend = c("True a(t)", "TKEO estimate"),
       col = c(4, 2), lwd = 2, bty = "n")

par(op)

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。