Rの関数:convolve {stats}

Rの関数から convolve {stats} を確認します。

関数 `convolve` とは

convolve は、高速フーリエ変換（FFT: Fast Fourier Transform）を利用して、2つの数値シーケンス（ベクトル）の「畳み込み（Convolution）」を計算する関数です。

時系列分析や信号処理において、畳み込みは「移動平均（平滑化フィルタ）の適用」や「2つの信号の類似度の計算（相互相関）」、「多項式の積の計算」など、多岐にわたる用途で使用されます。

FFTを利用することで、計算の計算量を通常の \(O(N^2)\) から \(O(N \log N)\) へと削減できるため、長い時系列データを扱う際に有用となります。

関数 `convolve` の引数

args(convolve)

function (x, y, conj = TRUE, type = c("circular", "open", "filter")) 
NULL

x
- 1つ目の数値シーケンス（ベクトル）。通常は「信号」や「元のデータ」を入力します。
- 数値型（実数または複素数）である必要があります。
y
- 2つ目の数値シーケンス（ベクトル）。「フィルタ（重み）」や「比較したいパターン」を入力します。
conj
- 計算過程で y のフーリエ変換に対して「複素共役（Complex Conjugate）」(複素数の虚数部分の符号（±）を逆にしたもの)をとるかどうかの論理値です。
- デフォルト: TRUE
type
- 計算方法および出力する結果の長さを指定する文字列です。
- "circular":
  - 巡回畳み込み（Circular Convolution）。x と y は完全に同じ長さである必要があり、端から終端へループするような計算を行います。
- "open":
  - 線形畳み込み（Linear Convolution）。計算前にゼロ埋め（ゼロパディング）を行い、データがループしないようにします。出力の長さは length(x) + length(y) - 1 になります。
- "filter":
  - open と同じ計算を行いますが、両端の「データとフィルタが完全に重なっていない部分（エッジ効果）」を切り落とし、有効な計算結果のみを返します。出力の長さは length(x) - length(y) + 1 となります（ただし length(x) >= length(y) の場合）。

なお、関数のヘルプには以下の一文が記されています。

Note that the usual definition of convolution of two sequences x and y is given by convolve(x, rev(y), type = “o”).（注：一般的な定義でのxとyの畳み込みは、convolve(x, rev(y), type = “o”) で求められます）
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/convolve.html

シミュレーションコード

以下に、convolve 関数の機能を確認するためのシミュレーションコードを示します。

引数（conj および type）の違いを確認するシミュレーション:
- 短い非対称なデータを用いて、計算結果の「長さ」と「向き」がどう変わるかを確認します。
信号処理（ノイズ除去）:
- ノイズを含む波形データに対し、ガウシアンフィルタを畳み込んで平滑化する実践例を示します。

引数 conj と type の違いを確認

# 基礎データの作成
# x_sig: 1箇所だけ値を持つインパルス信号
x_sig <- c(0, 0, 1, 0, 0, 0)
# y_fil: 非対称なフィルタ（重み）
y_fil <- c(10, 20, 30)

cat("入力信号 x_sig: ", paste(x_sig, collapse = ", "), "\n")
cat("フィルタ y_fil: ", paste(y_fil, collapse = ", "), "\n\n")

# --- 1-1. 引数 conj の違い (type = "open" で統一) ---
cat("--- conj = FALSE (畳み込み積分 / 多項式の積) ---\n")
res_conj_F <- convolve(x_sig, y_fil, conj = FALSE, type = "open")
cat(paste(round(res_conj_F, 1), collapse = ", "), "\n")

cat("\n--- conj = TRUE (デフォルト: 相互相関 / 移動平均的挙動) ---\n")
res_conj_T <- convolve(x_sig, y_fil, conj = TRUE, type = "open")
cat(paste(round(res_conj_T, 1), collapse = ", "), "\n")

入力信号 x_sig:  0, 0, 1, 0, 0, 0 
フィルタ y_fil:  10, 20, 30 

--- conj = FALSE (畳み込み積分 / 多項式の積) ---
0, 0, 0, 0, 10, 20, 30, 0 

--- conj = TRUE (デフォルト: 相互相関 / 移動平均的挙動) ---
0, 0, 30, 20, 10, 0, 0, 0

conj=FALSE ではフィルタがそのままの向き (10, 20, 30) で適用されています。

conj=TRUE ではフィルタが逆向き (30, 20, 10) で適用されています。

時系列の加重移動平均としてフィルタを設計した通りに機能させたい場合、conj=TRUE (デフォルト) にしておくか、入力するフィルタをあらかじめ rev() で反転させておく必要があります。

なお、本サンプルの場合であれば、以下のコードで convolve と同じ結果が得られます(引数 type を “open” とした場合)

`conj = FALSE` の場合

convolve_conj_false <- function(x, y) {
  nx <- length(x)
  ny <- length(y)
  out_len <- nx + ny - 1

  # 関数 convolve では計算前に、xの「左側」にフィルタの長さ-1個のゼロを追加しています。
  n1 <- ny - 1
  x_pad <- c(rep(0, n1), x)

  # yの「右側」にもゼロを追加し、長さを合わせます。
  y_pad <- c(y, rep(0, nx - 1))

  out_seq <- numeric(out_len)

  # FFTの掛け算は、時間領域での「循環畳み込み (Circular Convolution)」に相当します。
  # これをループで記述します。
  for (k in 0:(out_len - 1)) {
    sum_val <- 0
    for (j in 0:(out_len - 1)) {
      # 循環インデックス: はみ出した分は反対側からループさせる
      idx_y <- (k - j) %% out_len

      # Rは1番目から始まるため、インデックスに+1します
      sum_val <- sum_val + x_pad[j + 1] * y_pad[idx_y + 1]
    }
    out_seq[k + 1] <- sum_val
  }

  return(out_seq)
}

res <- convolve_conj_false(x_sig, y_fil)

paste(round(res, 1), collapse = ", ")

[1] "0, 0, 0, 0, 10, 20, 30, 0"

`conj = TRUE` の場合

convolve_conj_true <- function(x, y) {
  nx <- length(x)
  ny <- length(y)

  # 出力されるシーケンスの長さ (type = "open" の仕様)
  out_len <- nx + ny - 1
  out_seq <- numeric(out_len)

  # フィルタが端から入り込み、通り抜けるまでの全過程を計算するため、
  # 入力信号の前後に「余白（ゼロパディング）」を追加します。
  # 前後に (ny - 1) 個のゼロを追加します。
  x_pad <- c(rep(0, ny - 1), x, rep(0, ny - 1))

  # フィルタを1ステップずつスライドさせながら積和（内積）を計算します
  for (k in 1:out_len) {
    # スライド位置 k における、フィルタと同じ長さの入力部分（窓）を切り出します
    window_x <- x_pad[k:(k + ny - 1)]

    # conj=TRUE の場合、フィルタ y を「反転させずにそのままの向き」で掛け合わせます。
    # これにより、移動平均や相互相関（Cross-correlation）と同じ挙動になります。
    sum_val <- sum(window_x * y)

    out_seq[k] <- sum_val
  }

  return(out_seq)
}

res <- convolve_conj_true(x_sig, y_fil)

paste(round(res, 1), collapse = ", ")

[1] "0, 0, 30, 20, 10, 0, 0, 0"

引数 type の違いを確認 (conj = TRUE で統一)

cat("--- type = \"open\" (全ての重なりを計算) ---\n")
res_open <- convolve(x_sig, y_fil, type = "open")
cat("長さ: ", length(res_open), " -> ", paste(round(res_open, 1), collapse = ", "), "\n")

cat("\n--- type = \"filter\" (不完全なエッジを切り落とす) ---\n")
res_filter <- convolve(x_sig, y_fil, type = "filter")
cat("長さ: ", length(res_filter), " -> ", paste(round(res_filter, 1), collapse = ", "), "\n")

cat("\n--- type = \"circular\" (巡回・長さを合わせる必要あり) ---\n")
y_circ <- c(y_fil, 0, 0, 0) # 長さを x_sig (長さ6) と合わせる
res_circ <- convolve(x_sig, y_circ, type = "circular")
cat("長さ: ", length(res_circ), " -> ", paste(round(res_circ, 1), collapse = ", "), "\n")

--- type = "open" (全ての重なりを計算) ---
長さ:  8  ->  0, 0, 30, 20, 10, 0, 0, 0 

--- type = "filter" (不完全なエッジを切り落とす) ---
長さ:  4  ->  30, 20, 10, 0 

--- type = "circular" (巡回・長さを合わせる必要あり) ---
長さ:  6  ->  30, 20, 10, 0, 0, 0

open は、長さが (6 + 3 - 1 = 8) に拡張されます。
filter は、有効な重なりのみを取り出し、長さが (6 - 3 + 1 = 4) になります。
circular は、はみ出した部分が反対側からループして戻ってきます。

信号処理 (ノイズ除去フィルタ)

# 乱数シードの固定
seed <- 20260305
set.seed(seed)

# パッケージの読み込み
library(ggplot2)
library(gridExtra)

# 時系列データの作成 (サイン波 + ランダムノイズ)
n_time <- 200
time_idx <- 1:n_time
true_signal <- sin(2 * pi * time_idx / 50)
noise <- rnorm(n_time, mean = 0, sd = 0.5)
observed_signal <- true_signal + noise

# ガウシアンフィルタ（正規分布の形をした平滑化の重み）の作成
filter_length <- 15
x_seq <- seq(-3, 3, length.out = filter_length)
gaussian_filter <- dnorm(x_seq)
# フィルタの合計が1になるように正規化（信号の大きさを変えないため）
gaussian_filter <- gaussian_filter / sum(gaussian_filter)

# 畳み込みの実行 (type="filter", conj=TRUE)
# ノイズを含む信号にガウシアンフィルタをスライドさせながら適用します
smoothed_signal <- convolve(observed_signal, gaussian_filter, type = "filter")

# プロット用のデータ整形
# type="filter" は両端が切り落とされるため、元の時間軸と合わせるためのオフセットを計算
offset <- floor(filter_length / 2)
time_smoothed <- (1 + offset):(n_time - offset)

df_signal <- data.frame(
  Time = time_idx,
  Observed = observed_signal,
  TrueSignal = true_signal
)

df_smoothed <- data.frame(
  Time = time_smoothed,
  Smoothed = smoothed_signal
)

# 可視化
p1 <- ggplot() +
  # 観測データ (ノイズあり)
  geom_point(
    data = df_signal, aes(x = Time, y = Observed),
    color = "gray60", alpha = 0.6, size = 1.5
  ) +
  # 真の信号 (正解ライン)
  geom_line(
    data = df_signal, aes(x = Time, y = TrueSignal, color = "真の信号"),
    linewidth = 1, linetype = "dashed"
  ) +
  # 畳み込みによる平滑化結果
  geom_line(
    data = df_smoothed, aes(x = Time, y = Smoothed, color = "平滑化結果 (convolve)"),
    linewidth = 1.5
  ) +
  scale_color_manual(
    values = c("真の信号" = "blue", "平滑化結果 (convolve)" = "red")
  ) +
  labs(
    title = "convolve関数を用いた時系列データの平滑化",
    x = "時間 (Time)",
    y = "振幅 (Amplitude)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom", legend.title = element_blank())

print(p1)

灰色の点

プロット全体に散らばっている灰色の点は、センサーや計測器から得られた「生データ」を模したものです。

本来の滑らかな波形に対して標準偏差0.5のランダムノイズが付加されているため、上下に振動しており、元の規則正しい波の形を直接読み取ることは困難な状態です。

真の信号（青い破線）

青の破線は、ノイズが一切ない状態のサイン波（正解データ）を示しています。

現実のデータ分析においてはこの線は未知であり、観測データ（灰色の点）のみからこの形状を復元することが分析の目的となります。

平滑化結果（赤い実線）

この赤い実線が、convolve 関数によって観測データに「ガウシアンフィルタ（正規分布の形をした重み付き移動平均）」を適用した結果です。

ノイズの除去効果:
- 灰色の点が持つノイズが平滑化（スムージング）され、赤い線は青い破線（真の信号）に近い軌跡を描いています。畳み込み処理が、データの局所的な外れ値を吸収し、基調となるうねりだけを残したことが確認できます。
端の処理（エッジ効果）:
- 赤い線が時間軸の0付近および200付近で途切れています。これは、関数実行時に引数 type = "filter" を指定した結果です。フィルタ（長さ15）がデータに完全に重ならない両端の部分は「信頼性の低い計算結果」となるため、convolve 関数が自動的にその部分を切り落として出力していることを示しています。

以上です。

関数 convolve とは

関数 convolve の引数