【デジタルフィルタの基礎２】FIRのローパスフィルタの仲間たち—— Rで体験しよう

前回の記事では、FIR フィルタの基本について概観しました（以下のリンク）。

【デジタルフィルタの基礎１】まずは、FIRフィルタの世界を一望してみよう - ケィオスの時系列解析メモランダム

　しかし、フィルタを「使いこなせる」ようになるためには、抽象的な知識だけでは不十分です。

　フィルタは本来、教科書の数式に閉じ込められている概念ではなく、私たちが日常的に触れる数多くのデータを「見やすく」「意味のある形」に整えるための実用的な道具です。

雑音でガタガタした波形をなめらかにしたい
ピークや立ち上がりの形を保ったままノイズだけ除きたい
ある周波数成分だけを取り出したい
データに含まれるトレンドを見たい

こうした目的のために使うのが、ローパス FIR フィルタです。ローパスフィルタ (low-pass filter) とは、ゆっくり変化する成分（低周波成分: low-frequency components）だけを通し、速く変動する成分（高周波成分: high-frequency components）を抑えるフィルタのことです。

　そこで今回の記事では、FIR の基礎をふまえたうえで、代表的なローパス FIR フィルタとして、以下の４種とおまけの１種を紹介します。

Moving Average FIR
Savitzky–Golay FIR 平滑フィルタ
Windowed-Sinc FIR Low-Pass（Hamming 窓）
Frequency-Sampling FIR Low-Pass
(おまけ) Parks–McClellan の equiripple FIR

説明だけではブログとして書く意義がありませんので、ここでは、R を使ってインパルス応答と周波数応答を可視化したり、実データへ適用したいできるように、Rスクリプトの例も示しておきます。

　「簡単に使える道具として身近に感じ、体験を通じて使いこなせるようになる」ことを目指します。

　簡単さや便利さという点では、R のパッケージを使う方法もとても魅力的です。しかし本稿では、以下に示す R スクリプトの処理手順そのものが理解の助けになると考え、あえてパッケージは一切用いないことにします。スクリプトを部分ごとに実行しながら、「どのような計算をしているのか」「どの式がどの処理に対応しているのか」を追いかけてみてください。

0.【Rで必ず最初に実行】Rスクリプトを使うために準備
- ■ 描画のための関数
- ■ サンプル時系列の生成
1. Moving Average FIR（単純移動平均ローパス）
2. Savitzky–Golay FIR 平滑フィルタ
3. Windowed-Sinc FIR Low-Pass（Hamming 窓）
4. Frequency-Sampling FIR Low-Pass
5. おまけ：Parks–McClellan による equiripple FIR
6. おわりに

0.【Rで必ず最初に実行】Rスクリプトを使うために準備

　この記事で示すRスクリプトを実行する前に、以下の共通の２つのRスクリプトを一度だけ実行してください。ここで定義した関数を繰り返し使います。

■ 描画のための関数

　ますは、インパルス応答と周波数応答（振幅）を描画する関数です。

###############################################
# FIR のインパルス応答と周波数応答（振幅）を描画する関数
#  - h : FIR 係数ベクトル（0中心の非因果でも、0〜N-1 でも可）
#  - Fs: サンプリング周波数 [任意]（軸ラベル上は正規化周波数）
###############################################

plot_fir <- function(h, Fs = 1, main_prefix = "") {

  N <- length(h)

  # インデックス n を決める
  #  奇数長なら -M..M で 0 中心、偶数長なら 0..N-1
  if (N %% 2 == 1) {
    M <- (N - 1) / 2
    n <- -M:M
  } else {
    n <- 0:(N - 1)
  }

  # タイトル
  if (nzchar(main_prefix)) {
    main.time  <- paste0(main_prefix)
    main.freq  <- "Amplitude Response (0 to Nyquist)"
  } else {
    main.time  <- "FIR impulse response"
    main.freq  <- "Amplitude Response (0 to Nyquist)"
  }

  ################################################
  # 周波数応答 H(ω) を DTFT で計算
  #   H(ω) = Σ_n h[n] exp(-j ω n)
  ################################################
  n_freq <- 1024
  omega  <- seq(0, pi, length.out = n_freq)   # 0〜π
  freq   <- omega / (2 * pi)                  # 正規化周波数 (0〜0.5)

  H_omega <- sapply(omega, function(w) {
    sum(h * exp(-1i * w * n))
  })

  mag_raw <- Mod(H_omega)
  max_mag <- max(mag_raw, na.rm = TRUE)
  if (!is.finite(max_mag) || max_mag == 0) {
    mag_norm <- rep(0, length(mag_raw))
  } else {
    mag_norm <- mag_raw / max_mag
  }

  mag_safe <- pmax(mag_norm, 1e-12)          # log(0) 回避
  mag_db   <- 20 * log10(mag_safe)
  mag_db[mag_db < -80] <- -80                # -80 dB で打ち切り

  ################################################
  # プロット（左：インパルス応答 / 右：周波数応答）
  ################################################
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  par(mfrow = c(1, 2), cex.lab = 1.3, cex.axis = 1.1)

  ## 左：インパルス応答
  plot(
    n, h,
    type = "h", lwd = 3, col = "lightblue",
    xlab = "n", ylab = "h[n]",
    main = main.time,
    xaxs = "r"
  )
  points(n, h, pch = 16, col = "blue")

  ## 右：周波数応答（0〜0.5, -80〜0 dB）
  plot(
    freq, mag_db,
    type = "l", lwd = 2, col = "blue",
    xlab = "Normalized frequency (cycles/sample)",
    ylab = "Magnitude [dB]",
    main = main.freq,
    ylim = c(-80, 0),
    xaxs = "i"
  )
  abline(h = -3,  lty = 3, col = "gray")
  abline(h = -40, lty = 3, col = "gray")
}

■ サンプル時系列の生成

　フィルタリングを試すために、サンプル時系列を生成するＲスクリプトです。

###############################################
# Generate Input Signal
###############################################

n.peaks <- 8
sig     <- exp(seq(log(50), log(1), length = n.peaks))
mu      <- max(sig) * 6 * (1:n.peaks)
dt      <- 1

t.range <- seq(0, max(mu) + max(sig) * 3, by = dt)
n.t     <- length(t.range)

# Clean signal: normalized Gaussian peaks
x.clean <- numeric(n.t)
for (i in 1:n.peaks) {
  x.clean <- x.clean +
    dnorm(t.range, mean = mu[i], sd = sig[i]) /
    dnorm(mu[i],    mean = mu[i], sd = sig[i])
}

# Noisy signal
x    <- x.clean + rnorm(n.t, sd = 0.05)
time <- seq(0, (n.t - 1) * dt, by = dt)

1. Moving Average FIR（単純移動平均ローパス）

　Moving Average（移動平均）フィルタは、もっとも素朴な FIR ローパスです。移動平均は、もともと統計学・時系列解析（経済時系列の平滑化など）で古くから使われてきた手法であり、デジタル信号処理においても最初に登場する平滑化フィルタの一つです。

　窓幅 $\displaystyle{ M }$ の単純移動平均は

$\displaystyle{ h[n] = \frac{1}{M},\quad n=0,1,\dots,M-1 }$

それ以外は 0 というインパルス応答になります。これは「一定時間幅の単純平均」を取る操作そのものであり、

時間領域では：
- 因果型フィルタの場合は、「直前の $\displaystyle{M}$ 個のサンプルを均等に平均する」操作になります。
- 非因果型フィルタの場合は、「現在点を中心とした左右対称の区間にある $\displaystyle{M}$ 個のサンプルを均等に平均する」ことに相当します。
周波数領域では： 矩形窓のフーリエ変換である sinc 型の特徴をもつローパス特性を示します。

■ 特徴

線形位相 係数を左右対称に配置すれば、純粋な遅延だけを与える線形位相フィルタになります。
ノイズの平滑化に有効 白色ノイズであれば、分散が $\displaystyle{1/M}$ に縮小するので、ざっくりとしたノイズ低減には有効です。
周波数特性は良くない 矩形窓に対応するため、
- サイドローブが大きい（阻止域での減衰が弱い）
- 過渡帯域（通過域と阻止域の境界）が広いという意味で、「鋭いローパス」としては性能がよくありません。

■ Rで描くインパルス応答と周波数応答特性

　Rを使って、インパルス応答を計算し、周波数特性を描いてみましょう。これは、非因果型の移動平均フィルタについて、フィルタ長を M = 7 とした例です。

M <- 7
h_ma <- rep(1 / M, M)

plot_fir(h_ma, main_prefix = "Moving-average")

■ Rでフィルタリング

　時系列 x に対し、フィルタをかけてみましょう。上のRスクリプトを実行して、インパルス応答を計算した後に、このスクリプトを実行してください。

###############################################
# Moving Average FIR を x に適用して可視化
###############################################

# 非因果（左右対称）移動平均として適用（sides = 2）
y_ma <- as.numeric(stats::filter(x, filter = h_ma, sides = 2))

# 端の NA をとりあえず元のデータで埋める（簡便な処理）
half <- floor(M / 2)
y_ma[1:half] <- x[1:half]
y_ma[(length(y_ma) - half + 1):length(y_ma)] <- x[(length(x) - half + 1):length(x)]

# プロット
oldpar <- par(no.readonly = TRUE); on.exit(par(oldpar))
par(mfrow = c(2, 1), cex.lab = 1.2, cex.axis = 1.0,xaxs="i")

# 上段：元データとクリーン波形
plot(time, x, type = "l",
     main = "Input signal (noisy) and clean reference",
     xlab = "time", ylab = "amplitude", col = "gray70")
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Clean x.clean"),
       col = c("gray70", "red"), lwd = c(1, 2), bty = "n")

# 下段：移動平均フィルタ出力
plot(time, x, type = "l",
     main = "Moving-average FIR",
     xlab = "time", ylab = "amplitude", col = "gray80")
lines(time, y_ma, col = "blue", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Filtered"),
       col = c("gray80", "blue"),
       lwd = c(1, 2), lty = c(1, 1), bty = "n")

2. Savitzky–Golay FIR 平滑フィルタ

　Savitzky–Golay（S-G）フィルタは、「窓の中のデータに多項式をフィットして、その多項式の値を平滑化値とする」タイプの FIR フィルタです。これは、1964 年に Abraham Savitzky と Marcel J. E. Golay が「Smoothing and Differentiation of Data by Simplified Least Squares Procedures」と題する論文で提案した手法です。

当初は赤外／可視スペクトルなど「連続測定されたスペクトルデータ」の平滑化と微分のために開発され、
その後、Numerical Recipes などの文献を通じて広く普及し、現在では信号処理・化学計測・生体信号解析など、多くの分野で用いられています。

　長さ $\displaystyle{ 2m+1 }$ の窓をとり、相対インデックス $\displaystyle{ k=-m,\dots,0,\dots,m }$ に対して

$\displaystyle{ f(k) = c_0 + c_1 k + c_2 k^2 + \dots + c_p k^p }$

という $\displaystyle{ p }$ 次多項式を最小二乗法でフィットし、中心（ $\displaystyle{ k=0 }$ ）の値 ( $\displaystyle{ f(0)=c _ 0 }$ ) を「平滑化された値」とみなします。この「線形回帰の解」を行列表現で書き下すと、結局

$\displaystyle{ \tilde{x}[n] = \sum_{k=-m}^{m} h[k]\, x[n-k] }$

という、FIR畳み込みになます。

■ 特徴：波形の形を保つローパス

S-G フィルタは、周波数領域から設計されたものではなく、時間領域での「局所近似」の精度を最大化するよう設計されている点が特徴的です。

ピークの位置や幅を保ちやすい 多項式近似により、ピーク周辺のカーブを再現するため、単純移動平均に比べて「ピークがつぶれにくい」「形が保たれる」傾向があります。
線形位相 FIR 係数が対称になるため、位相特性は線形で、波形の時間方向の歪みを生じません。
ローパスとしては「ほどほど」 高周波ノイズはある程度抑えますが、
- 阻止域減衰は windowed-sinc や equiripple に比べると弱い
- 特定のカットオフ周波数を厳密に指定する設計ではないといった性質があります。

その代わり、スペクトルの細かな形よりも「波形の見た目」を重視した解析（生体信号、スペクトルのピーク解析など）に向いています。

■ Rで描くインパルス応答と周波数応答特性

　Rを使って、インパルス応答を計算し、周波数特性を描いてみましょう。以下は、m = 7、p = 2 の例です。

## Savitzky-Golay フィルタ係数を計算する関数 -------------------

sgolay_coeff <- function(m, p) {
  # m : 窓の片側の点数（窓幅は 2m+1）
  # p : 多項式次数
  # 例: m = 7, p = 2 なら窓幅 15 の 2 次多項式フィット

  # 相対インデックス k = -m, ..., m
  k <- -m:m
  s <- length(k)   # = 2m+1

  # デザイン行列 A: A_{ij} = k_i^{j-1}, j = 1..(p+1)
  A <- matrix(0, nrow = s, ncol = p + 1)
  for (i in 1:s) {
    for (j in 0:p) {
      A[i, j + 1] <- k[i]^j
    }
  }

  # 最小二乗解: 中心点 (k=0) の近似値に対する係数を求める
  # 中心点の値 = 多項式の 0 次係数 c0 → 単位ベクトル e0 = (1, 0, 0, ..., 0)
  e0 <- c(1, rep(0, p))

  # 通常の LS 理論から、FIR 係数ベクトル h は
  # h^T = e0^T (A^T A)^{-1} A^T  となる。
  # よって h = A (A^T A)^{-1} e0
  ATA_inv <- solve(t(A) %*% A)
  h <- A %*% (ATA_inv %*% e0)

  as.vector(h)   # 長さ 2m+1 の対称 FIR 係数
}

## Savitzky-Golay LPF (例: m=5, p=3) -----------------------------

m <- 7   # 片側 7 → 全体で 15 点
p <- 2   # 2 次多項式

h_sg <- sgolay_coeff(m = m, p = p)

plot_fir(h_sg, main_prefix = "Savitzky-Golay (p=3, width=11)")

Savitzky–Golay FIR 平滑フィルタのインパルス応答と周波数応答特性

■ Rでフィルタリング

###############################################
# Savitzky–Golay FIR を x に適用して可視化
###############################################

L_sg <- length(h_sg)

# 対称 FIR として適用（sides = 2）
y_sg <- as.numeric(stats::filter(x, filter = h_sg, sides = 2))

# 端の NA をとりあえず元のデータで埋める（簡便な処理）
half <- floor(L_sg / 2)
y_sg[1:half] <- x[1:half]
y_sg[(length(y_sg) - half + 1):length(y_sg)] <- x[(length(x) - half + 1):length(x)]

# プロット
oldpar <- par(no.readonly = TRUE); on.exit(par(oldpar))
par(mfrow = c(2, 1), cex.lab = 1.2, cex.axis = 1.0, xaxs = "i")

# 上段：元データとクリーン波形
plot(time, x, type = "l",
     main = "Input signal (noisy) and clean reference",
     xlab = "time", ylab = "amplitude", col = "gray70")
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Clean x.clean"),
       col = c("gray70", "red"), lwd = c(1, 2), bty = "n")

# 下段：Savitzky–Golay フィルタ出力
plot(time, x, type = "l",
     main = "Savitzky-Golay FIR",
     xlab = "time", ylab = "amplitude", col = "gray80")
lines(time, y_sg, col = "blue", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Filtered"),
       col = c("gray80", "blue"),
       lwd = c(1, 2), lty = c(1, 1), bty = "n")

3. Windowed-Sinc FIR Low-Pass（Hamming 窓）

　理想的なローパスフィルタの周波数応答は

通過域： $\displaystyle{|\omega| \le \omega _ c}$ で振幅 1
阻止域：( $\displaystyle{|\omega| \gt \omega _ c}$ ) で振幅 0

という「理想ブリックウォール形」を持ちます。そのインパルス応答を逆フーリエ変換すると、無限長の sinc 型の応答になります：

$\displaystyle{ h_{\text{ideal}}[n] = \begin{cases} \displaystyle \frac{\sin(\omega_c n)}{\pi n}, & n \ne 0,\\ \displaystyle \frac{\omega_c}{\pi}, & n = 0. \end{cases} }$

　しかし、これは

無限長
非因果（対称にとると $\displaystyle{n\in(-\infty,+\infty}$ )）

であり、そのままでは実装できません。

　そこで、

sinc のまわりを有限長 $\displaystyle{-m,\dots,m}$ で切り出し、
さらに「窓関数」 $\displaystyle{w[n]}$ を掛けて端をなめらかにする

という設計を行います：

$\displaystyle{ h[n] = h_{\text{ideal}}[n]\cdot w[n],\quad n=-m,\dots, m. }$

これが、windowed-sinc FIR ローパスであり、窓関数として Hamming, Hanning, Blackman, Kaiser などを選ぶことで

サイドローブ高さ（阻止域リップル）
主ローブ幅（過渡帯域の広さ）

のトレードオフを設計者の好みに合わせて調整できます。

　ここでは、Hamming 窓、

$\displaystyle{ w[n] = 0.54 - 0.46\cos\left(\frac{2\pi n}{N-1}\right) }$

を使いました。これは、windowed-sinc FIR の標準的な窓だそうです。

■ 特徴

線形位相 FIR 係数を対称に取るため、入力波形の形を時間方向に平行移動させるだけで保つことができます。
設計が簡単で、特性も良い カットオフ周波数とフィルタ長（窓幅）を決めるだけで、「それなりに良い」ローパスが得られるため、実用上のデフォルト選択肢として使われることが多いです。
等リップル最適化（Parks–McClellan）に比べるとやや鈍い 同じフィルタ長で比べると、Parks–McClellan による equiripple FIR の方が過渡帯域を狭くできますが、 windowed-sinc は計算が単純で理解しやすく、教育・実務面でバランスの良い方法です。

■ Rで描くインパルス応答と周波数応答特性

　Rを使って、インパルス応答を計算し、周波数特性を描いてみましょう。以下は、m=30、fc = 0.1の例です。

## Hamming 窓付き windowed-sinc ローパス -----------------------

# フィルタ長（奇数が扱いやすい）
m <- 30
M <- 2*m+1     # フィルタ長：M tap
fc <- 0.1            # 正規化カットオフ (0〜0.5 推奨)

n <- 0:(M - 1)

# 理想 LPF のインパルス応答 (非因果形, 中心が 0)
h_ideal <- numeric(M)
for (i in 1:M) {
  k <- n[i] - m
  if (k == 0) {
    h_ideal[i] <- 2 * fc
  } else {
    h_ideal[i] <- sin(2 * pi * fc * k) / (pi * k)
  }
}

# Hamming 窓
w <- 0.54 - 0.46 * cos(2 * pi * n / (M - 1))

# 窓付き FIR 係数
h_win <- h_ideal * w

# プロット
plot_fir(h_win, main_prefix = "Windowed-sinc (Hamming)")

Windowed-Sinc FIR Low-Pass（Hamming 窓）のインパルス応答と周波数応答特性

■ Rでフィルタリング

###############################################
# Windowed-sinc (Hamming) FIR を x に適用して可視化
###############################################

L_win <- length(h_win)

# 対称 FIR として適用（sides = 2）
y_win <- as.numeric(stats::filter(x, filter = h_win, sides = 2))

# 端の NA をとりあえず元のデータで埋める（簡便な処理）
half <- floor(L_win / 2)
y_win[1:half] <- x[1:half]
y_win[(length(y_win) - half + 1):length(y_win)] <- x[(length(x) - half + 1):length(x)]

# プロット
oldpar <- par(no.readonly = TRUE); on.exit(par(oldpar))
par(mfrow = c(2, 1), cex.lab = 1.2, cex.axis = 1.0, xaxs = "i")

# 上段：元データとクリーン波形
plot(time, x, type = "l",
     main = "Input signal (noisy) and clean reference",
     xlab = "time", ylab = "amplitude", col = "gray70")
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Clean x.clean"),
       col = c("gray70", "red"), lwd = c(1, 2), bty = "n")

# 下段：windowed-sinc フィルタ出力
plot(time, x, type = "l",
     main = "Windowed-sinc LPF (Hamming)",
     xlab = "time", ylab = "amplitude", col = "gray80")
lines(time, y_win, col = "blue", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Filtered"),
       col = c("gray80", "blue"),
       lwd = c(1, 2), lty = c(1, 1), bty = "n")

Windowed-Sinc FIR Low-Pass（Hamming 窓）を使ったフィルタリングの例

4. Frequency-Sampling FIR Low-Pass

　Frequency-Sampling 法は、周波数領域から直接 FIR を設計する古典的な方法です。この方法は、デジタルフィルタ設計の初期（1960年代〜）から知られている手法で、Rader & Gold らによる線形位相 FIR フィルタ設計の文献の中で登場します。

まだ計算機資源が限られていた時代には、「DFT グリッド上で少数のサンプルを与え、その IFFT を取る」という構成が実装しやすいというメリットがありました。
現在では、Parks–McClellan や least-squares などの最適設計法が普及していますが、任意形状フィルタの prototyping や教育用デモとして、今も有用な設計手法です。
まず、設計したいローパスの「理想周波数応答」 $\displaystyle{ H(\omega) }$ を周波数グリッド（たとえば DFT の格子）上のサンプル $\displaystyle{ \{H[k]\} }$ として与えます。
それに対して逆離散フーリエ変換（IDFT）を行い、インパルス応答 $\displaystyle{ h[n] }$ を得ます：

$\displaystyle{ h[n] = \frac{1}{N}\sum_{k=0}^{N-1} H[k] \, e^{j 2\pi kn/N}. }$

必要であれば、線形位相になるように $\displaystyle{ H[k] }$ を複素共役対称にしたり、位相を調整したりします。

この方法は、

周波数応答の形を直感的に指定できる
高域だけでなく、中間帯域に細かいピーク・谷を持つような「特殊な」フィルタも簡単に作れる

という意味で大きな自由度を持ちます。

■ 特徴

任意形状の周波数応答を実現しやすい 「この周波数ではゲイン 1、ここでは 0.5、ここでは 0」といった指定をそのまま DFT グリッド上で与えられます。
補間によりリップルが生じる 周波数サンプル点の間は補間（sinc の重ね合わせ）の形で埋まるため、理想応答と一致しない部分にリップルが現れます。
最適性という意味では Parks–McClellan に及ばない 通過域・阻止域リップルの最大値を理論的に最小化する設計（equiripple）と比べると、Frequency-Sampling 法は「指定したサンプルを通る補間」であって、誤差最小化という意味での最適性は保証されません。

■ Rで描くインパルス応答と周波数応答特性

　Rを使って、インパルス応答を計算し、周波数特性を描いてみましょう。以下は、N=56 (偶数)、'fc = 0.1' の例です。

## Frequency-sampling 法による LPF -----------------------------

N <- 56            # FFT / FIR の長さ (偶数)
fc <- 0.1           # 正規化カットオフ (0〜0.5)

# 周波数グリッド
k <- 0:(N - 1)
f <- (k / N)    # 0～1

# 理想 LPF の振幅特性：0〜fc を 1, fc〜Fs/2 を 0
H_desired <- ifelse(f <= fc, 1, 0)

# 複素共役対称にしたければもう少し細工するが、
# ここでは簡便にそのまま IFFT を取る
h_tmp <- fft(H_desired, inverse = TRUE) / N

# 実部を取り，循環シフトして「ほぼ線形位相」に近い形にする
h_fs <- Re(h_tmp)
# 中心を 0 に揃えるために N/2 だけシフト
h_fs <- c(h_fs[(N/2 + 1):N], h_fs[1:(N/2)])

plot_fir(h_fs, main_prefix = "Frequency-sampling LPF")

Frequency-Sampling FIR Low-Passのインパルス応答と周波数応答特性

■ Rでフィルタリング

###############################################
# Frequency-sampling FIR を x に適用して可視化
###############################################

L_fs <- length(h_fs)

# 対称 FIR として適用（sides = 2）
y_fs <- as.numeric(stats::filter(x, filter = h_fs, sides = 2))

# 端の NA をとりあえず元のデータで埋める（簡便な処理）
half <- floor(L_fs / 2)
y_fs[1:half] <- x[1:half]
y_fs[(length(y_fs) - half + 1):length(y_fs)] <- x[(length(x) - half + 1):length(x)]

# プロット
oldpar <- par(no.readonly = TRUE); on.exit(par(oldpar))
par(mfrow = c(2, 1), cex.lab = 1.2, cex.axis = 1.0, xaxs = "i")

# 上段：元データとクリーン波形
plot(time, x, type = "l",
     main = "Input signal (noisy) and clean reference",
     xlab = "time", ylab = "amplitude", col = "gray70")
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Clean x.clean"),
       col = c("gray70", "red"), lwd = c(1, 2), bty = "n")

# 下段：frequency-sampling フィルタ出力
plot(time, x, type = "l",
     main = "Frequency-sampling FIR LPF",
     xlab = "time", ylab = "amplitude", col = "gray80")
lines(time, y_fs, col = "blue", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Filtered"),
       col = c("gray80", "blue"),
       lwd = c(1, 2), lty = c(1, 1), bty = "n")

Frequency-Sampling FIR Low-Passを使ったフィルタリングの例

5. おまけ：Parks–McClellan による equiripple FIR

ここまでの記事では、「自分で式から係数を作る」タイプの FIR ローパスを４種類見てきました。最後におまけとして、実務ではよく使われる「等リップル最適設計（equiripple FIR）」 を、Parks–McClellan 法で一気に設計してしまう方法も紹介しておきます。

■ Parks–McClellan 法とは？

　Parks–McClellan 法は、Remez 交換アルゴリズムに基づく FIR フィルタ設計法で、

「通過域・阻止域での振幅誤差の最大値（Chebyshev ノルム）」を最小化
その結果、通過域・阻止域のリップル高さが、ほぼ等しく並ぶ（equiripple）特性になる
同じフィルタ長で比較すると、windowed-sinc や frequency-sampling よりも「過渡帯域を狭く・阻止域減衰を深く」できることが多い

という意味で、「与えた条件の範囲で、かなり頑張って最適化された FIR」を与えてくれる設計法です。ただし、アルゴリズムはかなり複雑なので、自前で書くよりも、パッケージに任せるのが現実的です。

■ Rで描くインパルス応答と周波数応答特性

　Rを使って、インパルス応答を計算し、周波数特性を描いてみましょう。ここでは、以下の条件

正規化カットオフ周波数：fc = 0.1（Nyquist = 0.5 のとき）
短めの遷移帯域：0.1 → 0.13
フィルタ長：61 tap（次数 order = 60）

で、equiripple ローパス FIR を設計してみます。

　R では、signal パッケージの remez() 関数が Parks–McClellan 法を実装しています。初回だけ install.packages("signal") を実行してください。

###############################################
# Parks–McClellan equiripple FIR
###############################################

# 初回のみインストール
# install.packages("signal")

library(signal)

# 条件設定 ------------------------------------
order <- 60   # フィルタ次数（→ tap 数 = order + 1 = 61）
fc    <- 0.1  # 正規化カットオフ (cycles/sample, Nyquist = 0.5)
tr    <- 0.03 # 遷移帯域の幅（0.1 → 0.13）

# remez() は 0〜1 を 0〜Nyquist に正規化した周波数を使う
# ここでは Nyquist = 0.5 → 割り算して [0,1] に変換
f_edge <- c(0, fc, fc + tr, 0.5) / 0.5   # = c(0, 0.4, 0.5, 1.0)

# 通過域・阻止域の理想振幅
a_edge <- c(1, 1, 0, 0)

# 通過域と阻止域の重み（阻止域を強めに）
w_band <- c(1, 10)

# Parks–McClellan で equiripple FIR を設計
h_pm <- as.numeric(
  remez(n = order, f = f_edge, a = a_edge, w = w_band)
)

# すでに定義済みの plot_fir() でインパルス応答と周波数特性を確認
plot_fir(h_pm, main_prefix = "Equiripple LPF (Parks–McClellan)")

f_edge では、f を 0〜1 の範囲で指定します。ここでの 1 は Nyquist 周波数（0.5 cycles/sample）に対応します。
a_edge はその周波数境界での理想振幅（1 or 0）。
w_band は各バンドの誤差にかける重みで、ここでは「阻止域のリップルを通過域より厳しめ」にしています。

plot_fir() を使うと、windowed-sinc や frequency-sampling と比べて、

通過域と阻止域のリップルがほぼ等しい高さでそろっている（equiripple）
同じくらいの tap 数でも、過渡帯域がやや鋭くなっている

といった違いが見て取れるはずです。

■ Rでフィルタリング

###############################################
# Equiripple FIR (Parks–McClellan) を x に適用して可視化
###############################################

L_pm <- length(h_pm)

# 対称 FIR として適用（sides = 2）
y_pm <- as.numeric(stats::filter(x, filter = h_pm, sides = 2))

# 端の NA をとりあえず元のデータで埋める（簡便な処理）
half <- floor(L_pm / 2)
y_pm[1:half] <- x[1:half]
y_pm[(length(y_pm) - half + 1):length(y_pm)] <- x[(length(x) - half + 1):length(x)]

# プロット
oldpar <- par(no.readonly = TRUE); on.exit(par(oldpar))
par(mfrow = c(2, 1), cex.lab = 1.2, cex.axis = 1.0, xaxs = "i")

# 上段：元データとクリーン波形
plot(time, x, type = "l",
     main = "Input signal (noisy) and clean reference",
     xlab = "time", ylab = "amplitude", col = "gray70")
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Clean x.clean"),
       col = c("gray70", "red"), lwd = c(1, 2), bty = "n")

# 下段：equiripple フィルタ出力
plot(time, x, type = "l",
     main = "Equiripple FIR LPF (Parks–McClellan)",
     xlab = "time", ylab = "amplitude", col = "gray80")
lines(time, y_pm, col = "blue", lwd = 2)
legend("topright",
       legend = c("Noisy x", "Filtered"),
       col = c("gray80", "blue"),
       lwd = c(1, 2), lty = c(1, 1), bty = "n")

6. おわりに

　FIR フィルタを設計するとき、「とりあえずカットオフ周波数を決めて、そのまま使う」だけで終わってしまうことが少なくありません。もちろん、カットオフ周波数やフィルタ長を決めることは大事ですが、それ自体が目的ではありません。

　本当に意識したいのは、「このフィルタで何を守り、何を捨てたいのか」というフィルタリングの目的です。ピークの形を崩さずにノイズだけを落としたいのか、ゆっくりしたトレンドだけを取り出したいのか、あるいは特定の帯域だけを強調したいのか──目的によって、選ぶべきフィルタは大きく変わります。

　同じカットオフ周波数でも、Moving Average、Savitzky–Golay、windowed-sinc、frequency-sampling、equiripple では「守られるもの」と「犠牲になるもの」が違います。ぜひ、R スクリプトを動かしながら、スペクトルだけでなく波形の変化も眺めてみてください。「このフィルタは自分のデータのどんな情報を残し、何を捨てているのか」を意識できるようになると、フィルタは単なる周波数の加工ではなく、「データからどんな物語を読み取るか」をデザインするための道具として、ぐっと頼もしく感じられるはずです。