ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【デジタルフィルタの基礎３】積分はローパス、微分はハイパス—— フラクタル時系列解析の基礎として学べ

デジタルフィルタ 数学 物理 確率過程 1/fゆらぎ R言語

これまで FIR フィルタの基本的な考え方を概観してきましたが、今回はその中でも特に基礎的でありながら、信号解析のあらゆる分野に通じる重要なテーマを扱います。
　それは、
「積分（累積和）」と「微分（差分）」という基本操作が、なぜローパス、ハイパスとして働くのか？
という問題です。
　離散信号における積分と微分は、単なる数式上の操作ではありません。波形がどのように変化し、どの成分が強調され、どの成分が抑えられるかという、波形の形そのものを支配する基本原理であり、これを理解しておくことはその後のフィルタ設計の直観を劇的に強めます。
　さらに重要なのは、積分と微分の性質は、長時間相関のスケーリング解析にも直結しているという点です。

これらは、累積和や差分によってスケーリング特性を変換したり、スケール依存の変動を抽出したりしています。
　したがって、今回扱う「積分と微分の周波数特性」は、単にフィルタの話にとどまらず、1/f ゆらぎ、長時間相関、スケーリング解析の理解にも直結する非常に重要な基礎となります。
　以下では、離散信号に限定し、無限長の累積和と有限長の累積和の違い、そして 一次差分がなぜ典型的な FIR ハイパスになるのかを、z 変換と周波数特性を用いて解説します。

無限長の累積和 (左)、有限長の累積和 (中)、一次差分 (右)のフィルタとしての周波数応答特性。

1. 「積分＝累積和」のローパス特性

　離散時間において積分に相当する操作は 累積和（running sum） です。累積和には、理想的に無限長で考えるか、現実的に有限長を扱うかで、違いがあります。デジタルフィルタの区別で言えば、無限は IIR (無限インパルス応答)、有限は FIR (有限インパルス応答)という違いがあります。

1.1 無限長の累積和は理想だが実用上は架空の話

　無限の過去から離散積分する場合は、次のように無限に積算する形です。 \displaystyle{ y[n] = y[n-1] + x[n] } これは漸化式の形ですが、これを和として展開すると、

\displaystyle{ y[n] = x[n] + x[n-1] + x[n-2] + x[n-3] + \cdots }

すなわち、

\displaystyle{ y[n] = \sum_{k=0}^{\infty} x[n-k] }

と書けます。この式は、過去すべての入力を無限に足し合わせていることを意味します。
　したがってインパルス応答

\displaystyle{ h[n] = 1 \quad (n = 0, 1, 2, \dots) }

となり、無限長のインパルス応答を持つため、典型的な IIR（Infinite Impulse Response）です。

[補足] z 変換の基礎

　ここからの解説で迷子になる学生を少しでも減らすため、念のため、「z 変換」を使って周波数応答を調べる方法について軽く説明しておきます。
　離散時間信号 \displaystyle{ x[n] } の z 変換は、次のように定義されます。

\displaystyle{ X(z) = \sum_{n=-\infty}^{\infty} x[n]\, z^{-n} }


　線形時不変系 (差分方程式の係数が定数の場合)に入力 \displaystyle{ x[n] } を与えたとき、出力 \displaystyle{ y[n] } の z 変換 \displaystyle{ Y(z) } との比

\displaystyle{ H(z) = \frac{Y(z)}{X(z)} }

を 伝達関数（transfer function） と呼びます。
　周波数特性は、複素平面上の単位円

\displaystyle{ z = e^{j\omega} = e^{j 2\pi f} }

を代入することで求まります。したがって、

\displaystyle{ H(e^{j\omega}) = H(e^{j 2\pi f}) }

が、周波数 f [Hz] に対する振幅・位相応答になります。
　z 変換の基本性質で良く使う公式が、時間シフト（time shift）の公式です。
入力 \displaystyle{ x[n] } の 1 サンプル遅延信号 \displaystyle{ x[n-1] } の z 変換は：

\displaystyle{ x[n-1] \;\;\Longleftrightarrow\;\; z^{-1} X(z) }

　一般に \displaystyle{ k } サンプルシフトで、

\displaystyle{ x[n-k] \;\;\Longleftrightarrow\;\; z^{-k} X(z) }

が成り立ちます。

1.2 無限長の累積和の周波数特性

　話は戻り、無限長の累積和は

\displaystyle{ y[n] = y[n-1] + x[n] }

と書けることまで説明しました。この差分方程式に z 変換を適用すると、

\displaystyle{ Y(z) = z^{-1} Y(z) + X(z) }

となるので、伝達関数 \displaystyle{ H(z) } は

\displaystyle{ H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1 - z^{-1}} }

となります。
　このとき、複素平面上の単位円 \displaystyle{ z = e^ {j 2 \pi f} } 上での値をとることで周波数特性が得られます。すなわち、

\displaystyle{ H(e^{j2\pi f}) = \frac{1}{1 - e^{-j 2\pi f}} }

となります。大きさは

\displaystyle{ \left|H(e^{j2\pi f})\right| = \frac{1}{2 |\sin(\pi f)|} }

となります。
　ここでは、サンプリング周波数を１で考えていますので、意味のある周波数 \displaystyle{ f } は、0 から 1/2（ナイキスト周波数）の範囲です。したがって、無限長の累積和には以下と性質があることがわかります。

  • \displaystyle{f \to 0}（ゼロ周波数＝平均値）では
\displaystyle{ \sin(\pi f) \approx \pi f \;\Rightarrow\; |H| \to \infty }

→ 平均成分を非常に強く通過させる。

\displaystyle{ \sin(\pi f) \to \sin\left(\frac{\pi}{2}\right)=1 }

→ 高周波成分のゲイン (振幅の倍率)は最小。

■ DFAなどのフラクタル時系列解析とな関係

　DFA（detrended fluctuation analysis）などのスケーリング解析では、最初に時系列を「積分する」手続きがあります。では、なぜわざわざ積分する必要があるのでしょうか。
　もし DＦA を使わず、パワースペクトルから直接スケーリング指数を推定するのであれば、時系列を積分などする必要はありません。周波数 \displaystyle{ f } が小さい領域でパワースペクトルが

\displaystyle{ S(f) \sim f^{-\beta} }

のように振る舞う場合、両対数プロットの直線領域の傾き \displaystyle{ \hat{\beta} } を求めれば、スケーリング指数 \displaystyle{ \beta } を推定できます。
　一方、DFAでは推定可能なスケーリング指数に 下限 が存在します。DFA で得られるスケーリング指数 \displaystyle{ \alpha } はパワースペクトルの指数 \displaystyle{ \beta } と

\displaystyle{ \alpha = \frac{\beta + 1}{2} }

という関係をもつと覚えている方がいると思います。ただし注意すべきなのは、この式が \displaystyle{ \alpha \gt 0 } の場合にのみ成立することです。すなわち、DFAでは \displaystyle{ \alpha \le 0 } を推定できません。
　長時間相関を示す定常過程を前提とすると、パワースペクトルのスケーリング指数については、

\displaystyle{ -1 < \beta < 1 }

の範囲を考えることになります。したがって、上の関係式から対応する

\displaystyle{ 0 < \alpha < 1 }

を推定できそうに見えます。しかし、それが可能になるのは、DFA が最初に時系列を積分しているからです。
　もしDFAで時系列を積分しなければ、対応する指数は

\displaystyle{ -1 < \alpha < 0 }

となり、DFAの推定可能範囲の外に出てしまいます。したがって、DFAで時系列を積分する目的は、スケーリング指数 \displaystyle{ \alpha } を 1 だけ押し上げて、\displaystyle{ \alpha \gt 0 } の領域に移すためなのです。
　この理由を教えてくれるのが、積分フィルタの周波数特性です。積分フィルタの振幅特性は、周波数 \displaystyle{ f } が小さい領域で

\displaystyle{ \left|H(e^{j2\pi f})\right| = \frac{1}{2|\sin(\pi f)|} \approx \frac{1}{2\pi} |f|^{-1} }

となります。パワースペクトルは振幅の二乗に比例する量なので、ある時系列のパワースペクトルが \displaystyle{ S(f) } であるとき、これを積分すると新しいパワースペクトルは

\displaystyle{ S_{\text{int}}(f) = \left|H(e^{j2\pi f})\right|^{2} S(f) \approx \frac{1}{(2\pi)^{2}} |f|^{-2}\, S(f) }

となります。
　つまり、元の時系列のスペクトルが

\displaystyle{ S(f) \sim f^{-\beta} }

であれば、積分後のスペクトルは

\displaystyle{ S_{\text{int}}(f) \sim f^{-(\beta + 2)}. }

となり、スペクトルの傾きが 2 だけ急になる（スケーリング指数 \betaが 2 増える）ことがわかります。スペクトル指数が 2 増えるということは、スケーリング指数 \displaystyle{ \alpha } が 1 だけ増えることと同値です。
　したがって、DFAにおける「積分」は、単なる前処理ではなく、本来のスケーリング指数を DFA が扱える領域に移すための変換操作なのです。

1.3. 有限長の累積和は実用的な FIR ローパス

　実際の信号処理では、無限の過去の和をとるのではなく、最近の M サンプルだけを足し合わせるという、有限長の積分を用います。つまり、

\displaystyle{ y[n] = \sum_{k=0}^{M-1} x[n-k] }

です。これを、\displaystyle{ 1/M } 倍したものは、単純な 移動平均フィルタ、

\displaystyle{ y[n] = \frac{1}{M} \sum_{k=0}^{M-1} x[n-k] }

に一致します。これは、以前の記事で説明した通り、インパルス応答が有限長（長さ M）なので、完全な FIR です。そして、平均の計算は、低周波を通すローパスフィルタになります。

1.4 有限長の累積和の周波数特性

　有限長の累積和

\displaystyle{ y[n] = \sum_{k=0}^{M-1} x[n-k] }

の z 変換は、

\displaystyle{ H(z) = (1 + z^{-1} + z^{-2} + \cdots + z^{-(M-1)}) }

です。したがって、周波数応答は

\displaystyle{ H(e^{j2\pi f}) = \sum_{k=0}^{M-1} e^{-j 2\pi f\, k}= \frac{1 - e^{-j2\pi f M}}{1 - e^{-j2\pi f}} }

となります。
　ゲインについては、低周波で大きく、高周波で小さくなるローパス特性を持ちます。最初に示した図を見てください。さらに、おまけのRスクリプトを実行して、周波数応答のグラフを確認しておいてください。

■ 無限長と有限長の積分フィルタの特性の違い

　無限長の累積和と有限長の累積和の周波数特性は、分母は同じですが、分子の性質がまったく異なります。
無限長の場合、周波数応答は

\displaystyle{ H_{\infty}(f) = \frac{1}{\,1 - e^{-j 2\pi f}\,} }

となり、分子が常に 1 のため、低周波では

\displaystyle{ |H_{\infty}(f)| \approx (2\pi |f|)^{-1} }

となって無限に発散する「理想積分器」になります。
　一方、有限長 M 点の累積和では、

\displaystyle{ H_{M}(f) = \frac{\,1 - e^{-j 2\pi f M}\,}{\,1 - e^{-j 2\pi f}\,} }

となります。ここで分子

\displaystyle{ 1 - e^{-j 2\pi f M} }

の絶対値は

\displaystyle{ \left|1 - e^{-j 2\pi f M}\right| = 2\left|\sin(\pi f M)\right| }

のように、0 から 2 の間で周期的に振動します。
　そのため、低周波でも（ f \to 0 の極限で）

\displaystyle{ |H_{M}(f)| \approx M }

という有限の値に収まり、発散しません。
　要するに、有限長の積分フィルタでは、分子が振動するおかげで、無限長の場合に生じる低周波での無限大のゲインが自然に打ち消され、現実的なローパスとしての特性を示すようになるのです。

2. 「微分＝差分」のハイパス特性

　積分がローパスなら、微分に相当する操作はハイパスになります。

2.1 微分はとても単純：差分の計算を周波数の世界で見る

　離散微分は 一次差分で与えられます：

\displaystyle{ y[n] = x[n] - x[n-1]. }

　これはインパルス応答が、\displaystyle{ 1, -1 } の 2 点だけなので、有限長のFIRでとても単純な形です。

2.2 微分の周波数特性

　差分の計算

\displaystyle{ y[n] = x[n] - x[n-1]. }

を z 変換すると、

\displaystyle{ H(z) = 1 - z^{-1} }

となります。
　周波数特性は、

\displaystyle{ |H(e^{j2\pi f})| = 2,\left|\sin\left(\frac{2\pi f}{2}\right)\right| }

となります。 　 　したがって、\displaystyle{ f=0 } で 0、\displaystyle{ f=1/2 }（ナイキスト周波数）で最大となるハイパス特性を示します。最初に示した図を見てください。さらに、おまけのRスクリプトを実行して、周波数応答のグラフを確認しておいてください。

無限長の累積和 (左)、有限長の累積和 (中)、一次差分 (右)のフィルタとしての周波数応答特性（両対数プロット）。最後に示したRスクリプトで描いたグラフ。

5. まとめ

　今回は、離散信号における「積分＝累積和」「微分＝差分」が、なぜローパス・ハイパスとして働くのかを、z 変換と周波数特性を使って説明しました。無限長（IIR）と有限長（FIR）の累積和の違い、累積和フィルタのローパス性、そして一次差分がもっとも単純なハイパス FIR になることを確認してきました。

　信号処理に馴染みがない読者にとっては、「微分・積分を周波数特性で語られても正直ピンとこない」という感覚があると思います。実際、私自身も理論物理のバックグラウンドで生体信号研究を始めた頃、微積分といえば「時間領域の連続極限操作」というイメージが強く、フィルタとしての動作（ローパス・ハイパス）という視点で見たことはほとんどありませんでした。

　ところが、実際の時系列解析では、微分や積分が、どの周波数成分を強調し、どの成分を抑制するかというフィルタとしての特性を持ちます。この視点を得ると、DFA や DMA、Higuchi 法、さらにはパワースペクトルによる 1/f ゆらぎ解析など、多くのスケーリング解析手法がなぜあのような操作を行うのか、その背後にある「波形のどこを取り出そうとしているか」が、直感的に理解できるようになります。

　信号処理は、難しそうに見えて、実は「波形の形をどう変えるか」という非常に素朴な操作の積み重ねです。微分は“速い変化に敏感”、積分は“ゆっくりした成分を大事にする”。この古典的な性質が、周波数特性としてきれいに表現されるからこそ、解析手法の設計や結果の解釈も格段にやりやすくなります。

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

ランキング参加中
数学・科学・工学
ランキング参加中
知識
ランキング参加中
テクノロジー

[おまけ] Rで周波数特性を描画

　今回説明した積分・微分の周波数特性のグラフを描くRスクリプトを掲載しておきます。両対数のプロってになっていることに注意してください。

############################################################
# Frequency responses of:
#  1) Infinite-length cumulative sum (IIR "integrator")
#  2) Finite-length cumulative sum (FIR LPF, unnormalized)
#  3) First difference (FIR HPF)
############################################################

# 周波数軸の生成（点数を増やして解像度アップ）
Nf <- 1000
f  <- seq(1e-4, 0.5, length.out = Nf)

############################################################
# 1. 無限長の累積和 (IIR)
############################################################
H_inf_mag <- 1 / (2 * abs(sin(pi * f)))

############################################################
# 2. 有限長の累積和 (FIR, 非正規化)
############################################################
M <- 1001
Nf <- 10000
ff  <- seq(1e-4, 0.5, length.out = Nf)
H_cum_mag <- abs(sin(pi * M * ff) / sin(pi * ff))

############################################################
# 3. 一次差分 (FIR ハイパス)
############################################################
H_diff_mag <- 2 * abs(sin(pi * f))

############################################################
# log(0) を避ける
############################################################
eps <- 1e-12
H_inf_plot  <- pmax(H_inf_mag,  eps)
H_cum_plot  <- pmax(H_cum_mag,  eps)
H_diff_plot <- pmax(H_diff_mag, eps)

############################################################
# 対数軸補助関数
############################################################
log_axis <- function(side, lim, cex.axis = 1.4,
                     tck.major = -0.02, tck.minor = -0.01) {

  lim <- lim[lim > 0]
  if (length(lim) != 2) return(NULL)

  kmin <- floor(log10(min(lim)))
  kmax <- ceiling(log10(max(lim)))

  major <- 10^(kmin:kmax)
  minor <- as.vector(outer(major, 1:9, "*"))
  minor <- minor[minor >= min(lim) & minor <= max(lim)]

  axis(side, at = minor, labels = FALSE, tck = tck.minor)

  labs <- parse(text = paste0("10^", kmin:kmax))
  axis(side, at = major, labels = labs,
       las = ifelse(side %in% c(2, 4), 1, 0),
       cex.axis = cex.axis, tck = tck.major)
}

############################################################
# 両対数プロット
############################################################
old_par <- par(no.readonly = TRUE)
par(mfrow = c(1, 3), mar = c(6, 7, 4, 1), xaxs = "i", yaxs = "i", cex.lab=1.5)

xlim <- range(f)

############################################################
# (1) 無限長累積和
############################################################
ylim1 <- range(H_inf_plot)

plot(f, H_inf_plot,
     type = "l", col = 2, lwd = 2,
     xaxt = "n", yaxt = "n", log = "xy",
     xlim = xlim, ylim = ylim1,
     xlab = "Frequency f (normalized, 0 ... 0.5)",
     ylab = "")

log_axis(1, xlim)
log_axis(2, ylim1)

# ★ y軸ラベルを mtext で記述（数式表記）
mtext(side = 2, line = 4.5, cex = 1.3,
      expression( "|" * H[inf](e^{j*2*pi*f}) * "|" ) )

title("Infinite-length cumulative sum\n(IIR-like integrator)")

############################################################
# (2) 有限長累積和（FIR）
############################################################
ylim2 <- c(1, M)

plot(ff, H_cum_plot,
     type = "l", col = 2, lwd = 2,
     xaxt = "n", yaxt = "n", log = "xy",
     xlim = xlim, ylim = ylim2,
     xlab = "Frequency f (normalized, 0 ... 0.5)",
     ylab = "")

log_axis(1, xlim)
log_axis(2, ylim2)

# ★ y軸 mtext
mtext(side = 2, line = 4.5, cex = 1.3,
      expression( "|" * H[cum](e^{j*2*pi*f}) * "|" ) )

title(paste0("Finite-length cumulative sum\n(FIR, M = ", M, ")"))

############################################################
# (3) 一次差分（FIR HPF）
############################################################
ylim3 <- range(H_diff_plot)

plot(f, H_diff_plot,
     type = "l", col = 2, lwd = 2,
     xaxt = "n", yaxt = "n", log = "xy",
     xlim = xlim, ylim = ylim3,
     xlab = "Frequency f (normalized, 0 ... 0.5)",
     ylab = "")

log_axis(1, xlim)
log_axis(2, ylim3)

# ★ y軸 mtext
mtext(side = 2, line = 4.5, cex = 1.3,
      expression( "|" * H[diff](e^{j*2*pi*f}) * "|" ) )

title("First difference\nFIR high-pass")

par(old_par)