信号処理の要!Python NumPyで相互相関を極める完全ガイド
numpy.correlate()
は、NumPy ライブラリが提供する関数で、2つの1次元シーケンス(配列)間の相関を計算します。ここでいう「相関」は、統計学的な意味での相関係数とは少し異なり、信号処理の分野で使われる**相互相関(cross-correlation)**を指します。
相互相関とは?
相互相関は、ある信号を別の信号に対してずらしながら、その類似度を測るためのものです。直感的に言うと、あるパターンが別のパターンの中にどれだけ隠れているか、あるいは、2つの波形が時間的にどれくらい似ているかを定量化するのに使われます。
numpy.correlate()
の基本的な使い方
numpy.correlate(a, v, mode='valid')
-
mode
: 相関の計算方法を指定します。以下の3つのモードがあります。'valid'
(デフォルト): 出力の長さはmax(M, N) - min(M, N) + 1
となります。ここで M はa
の長さ、N はv
の長さです。これは、v
がa
の完全に重なる部分のみで計算され、結果の配列の長さが最も短くなるモードです。このモードでは、境界効果が最も少なくなります。'same'
: 出力の長さはmax(M, N)
となります。結果の中心が出力配列の中心になるようにパディングが行われます。出力は入力のうち長い方の配列と同じ長さになります。'full'
: 出力の長さはM + N - 1
となります。これは、v
がa
とわずかにでも重なるすべての位置で計算が行われるため、結果の配列の長さが最も長くなるモードです。
-
v
: 2番目の入力配列(1次元) -
a
: 最初の入力配列(1次元)
計算の仕組み(イメージ)
numpy.correlate(a, v)
は、基本的に v
を反転(フリップ)させてから、a
に対して畳み込み(convolution)を実行するのと同等です。
もし a=[a0​,a1​,...,aM−1​] と v=[v0​,v1​,...,vN−1​] の2つの配列があったとします。
numpy.correlate(a, v)
の出力 ckは、大まかには以下のように計算されます。
ck​=i=0∑N−1​ai+k​⋅vi​
ただし、これは v
が反転していない場合の畳み込みの定義と似ていますが、相互相関では v
が反転していることに注意が必要です。より厳密には、numpy.correlate(a, v)
は np.convolve(a, v[::-1])
と同じ結果を返します。(ここで [::-1]
は配列の反転を意味します。)
numpy.correlate()
の応用例
- 音響処理: エコーの検出など。
- 時系列分析: 2つの時系列データ間の関連性の分析。
- 画像処理: 画像内の特定パターンの検出(テンプレートマッチング)。
- 信号処理: 信号の遅延時間の推定、パターン検出、ノイズ除去など。
import numpy as np
# 例1: シンプルな相関
a = np.array([1, 2, 3, 4, 5])
v = np.array([1, 2, 3])
# mode='valid' (デフォルト)
# v が a の完全に重なる部分のみ計算
# (1*3 + 2*2 + 3*1) = 3+4+3 = 10
# (2*3 + 3*2 + 4*1) = 6+6+4 = 16
# (3*3 + 4*2 + 5*1) = 9+8+5 = 22
corr_valid = np.correlate(a, v, mode='valid')
print(f"mode='valid': {corr_valid}") # 出力: [10 16 22]
# mode='same'
# 出力の長さは max(len(a), len(v)) = 5
corr_same = np.correlate(a, v, mode='same')
print(f"mode='same': {corr_same}") # 出力例: [ 4 10 16 22 20] (パディングにより両端の値が小さくなる)
# mode='full'
# 出力の長さは len(a) + len(v) - 1 = 5 + 3 - 1 = 7
corr_full = np.correlate(a, v, mode='full')
print(f"mode='full': {corr_full}") # 出力例: [ 3 10 16 22 20 15 5]
# 例2: パターン検出
signal = np.array([0, 0, 1, 2, 3, 2, 1, 0, 0])
pattern = np.array([1, 2, 1])
# パターンがどこにあるかを探す
# 相関値が高いところがパターンと類似している場所
corr_pattern = np.correlate(signal, pattern, mode='valid')
print(f"\nパターン検出の相関: {corr_pattern}")
# 出力例: [ 8 10 8] -> 中央の値 (10) が最も高いため、[1,2,3,2,1] の中の [1,2,1] と一番相関が高いことを示唆
# 実際のインデックスは相関配列の最大値のインデックスから計算する必要がある
# 最大値のインデックス
max_idx = np.argmax(corr_pattern)
print(f"最大値のインデックス: {max_idx}")
# 'valid' モードの場合、元の信号の開始インデックスは max_idx
# 信号のどこにパターンがあるか(おおよその位置)
# 注意: 'valid'モードの場合、結果のインデックスと元の信号のインデックスの関係を正しく理解する必要があります。
# 上記の例では、[1,2,1] のパターンが信号の [1,2,3,2,1] の部分と最も高い相関を持ちます。
numpy.correlate()
は信号処理などで非常に便利な関数ですが、いくつか誤解しやすい点や、予期せぬ結果につながる落とし穴があります。
相関係数と混同する (Common Misconception: Confusing with Correlation Coefficient)
エラーの原因
numpy.correlate()
は、統計学で一般的に使われる**相関係数(correlation coefficient)とは異なる概念である相互相関(cross-correlation)**を計算します。相関係数は −1 から 1 の間の正規化された値ですが、numpy.correlate()
の出力は正規化されていません。そのため、値のスケールが大きく異なることがあります。
トラブルシューティング
-
別の関数を使う
統計的な相関係数が欲しい場合は、numpy.corrcoef()
やscipy.stats.pearsonr()
など、目的に合った別の関数を使用してください。 -
正規化
もし相関係数のような−1から1の範囲の値が欲しい場合は、numpy.correlate()
の結果を別途正規化する必要があります。一般的な正規化方法としては、それぞれの入力配列のノルム(大きさ)で割る方法があります。例:
import numpy as np a = np.array([1, 2, 3, 4, 5]) v = np.array([1, 2, 3]) corr = np.correlate(a, v, mode='full') # 正規化の一例 (相互相関係数の定義に近い) # ここでは、scipy.signal.correlate の normalized=True に近い動作を目指す # ただし、厳密な正規化は文脈によって異なる normalized_corr = corr / (np.linalg.norm(a) * np.linalg.norm(v)) print(f"非正規化された相関: {corr}") print(f"正規化された相関 (一例): {normalized_corr}")
-
定義の理解
numpy.correlate()
はa
とv
の要素の積の和を計算していることを理解してください。これは、v
を反転させてa
と畳み込み (convolution) を行うのと同等です。
出力配列の長さとインデックスの解釈 (Output Array Length and Index Interpretation)
エラーの原因
mode
引数('valid'
、'same'
、'full'
)によって出力配列の長さが異なるため、結果のインデックスがどの「ラグ」(時間的なずれ)に対応するのかを誤解することがあります。特に、最大相関値のインデックスが、そのまま望む「ずれ量」を示すとは限りません。
トラブルシューティング
-
np.argmax() の使用時の注意
np.argmax()
は相関値が最大となるインデックスを返しますが、そのインデックスが意味する実際のラグを正しく計算する必要があります。import numpy as np a = np.array([0, 0, 1, 2, 3, 2, 1, 0, 0]) v = np.array([1, 2, 1]) # フルモードで相関を計算 corr_full = np.correlate(a, v, mode='full') print(f"フルモードの相関: {corr_full}") # 最大値のインデックス max_idx = np.argmax(corr_full) print(f"最大値のインデックス: {max_idx}") # ラグの計算 (mode='full' の場合) # v を a の左端からスライドさせると考えると、 # 0番目の出力は v が a の len(v)-1 個分左にずれている状態 (a[0]*v[len(v)-1]...) # したがって、ラグは max_idx - (len(v) - 1) lag = max_idx - (len(v) - 1) print(f"推定されるラグ: {lag}") # この例では 2 となるはず (a[2]からパターンが始まるため)
-
図で考える
2つの配列がどのようにスライドして積和が計算されるかを視覚的に描いてみると、インデックスとラグの関係を理解しやすくなります。 -
mode='valid' の理解
valid
モードは、v
がa
に完全に含まれる範囲でのみ計算されるため、出力は最も短くなります。このモードでの出力インデックスは、元の配列の開始インデックスと密接に関連します。 -
mode='full' の理解
mode='full'
は、v
がa
とわずかにでも重なるすべての可能な位置で相関を計算するため、最も長い出力配列を返します。この場合、出力配列の中心付近がゼロラグ(ずれなし)に対応し、両端に行くほど大きなラグに対応します。- 出力配列のインデックス k が示すラグは、
full
モードの場合、通常$k - (len(a) - 1)$
もしくは$k - (len(v) - 1)$
のように計算されます(どちらの配列が長いかによる)。正確なラグの計算式はドキュメントを参照するか、簡単な例で試して確認するのが確実です。 - 例えば、
numpy.correlate(a, v, mode='full')
の場合、結果のインデックスi
に対応するラグ(v
がa
に対してどれだけずれているか)は$i - (len(v) - 1)$
となります。
- 出力配列のインデックス k が示すラグは、
入力配列が1次元ではない (Input Arrays are Not 1-Dimensional)
エラーの原因
numpy.correlate()
は1次元配列専用の関数です。2次元以上の配列を渡すとエラーが発生します。
トラブルシューティング
- 画像処理など2次元相関の場合
画像処理などで2次元の相互相関が必要な場合は、scipy.signal.correlate2d
のような関数を使用してください。 - リシェイプまたはスライス
2次元以上の配列の一部を相関させたい場合は、適切な行や列をスライスして1次元配列として渡すか、array.flatten()
などで1次元に変換してから渡します。 - 次元の確認
array.ndim
やarray.shape
を使って入力配列の次元を確認してください。
実行速度の問題 (Performance Issues)
エラーの原因
numpy.correlate()
は、内部的にFFT(高速フーリエ変換)を使用せず、直接的な積和演算で相関を計算します。このため、非常に長い配列(数万以上の要素を持つ配列など)に対しては計算が遅くなることがあります。
トラブルシューティング
-
scipy.signal.correlate の利用
大規模な配列の相互相関を計算する場合は、FFT を使用して高速に計算するscipy.signal.correlate
の使用を検討してください。これは、numpy.correlate
と同じようなインターフェースを持っています。import numpy as np from scipy.signal import correlate a = np.random.rand(100000) v = np.random.rand(1000) # NumPyのcorrelate # %timeit np.correlate(a, v, mode='full') # SciPyのcorrelate (FFT使用) # %timeit correlate(a, v, mode='full') # 一般的にSciPyの方が高速
データ型に関する問題 (Data Type Issues)
エラーの原因
NumPy はデータ型に厳密なため、意図しないデータ型(例えば整数型)で計算を行うと、期待する浮動小数点数の結果が得られないことがあります。特に、正規化を行う際に問題になることがあります。
トラブルシューティング
- 浮動小数点数への変換
入力配列が整数型の場合、明示的に浮動小数点数に変換してからcorrelate()
を実行してください。a_int = np.array([1, 2, 3]) v_int = np.array([1, 2, 1]) # 整数型で計算すると結果も整数になる可能性がある (NumPyのバージョンや演算による) corr_int = np.correlate(a_int, v_int) print(f"整数型の相関: {corr_int}") # 浮動小数点数に変換してから計算 a_float = a_int.astype(float) v_float = v_int.astype(float) corr_float = np.correlate(a_float, v_float) print(f"浮動小数点数の相関: {corr_float}")
注意点
numpy.correlate()
は複素数配列も扱うことができますが、その定義は v
の共役複素数 (v.conj()
) との畳み込みになります。信号処理の文脈ではこれが一般的ですが、数学的な定義によっては異なる場合があるため、注意が必要です。
トラブルシューティング
- ドキュメントの確認
複素数配列を扱う際は、NumPy の公式ドキュメントで定義を確認し、意図した計算が行われているかを確認してください。
numpy.correlate()
は、主に信号処理やパターンマッチングで2つの1次元シーケンス間の類似度(相互相関)を計算するために使用されます。
例1: 基本的な相互相関の計算 (mode='valid'
)
最も基本的な使い方です。v
が a
に完全に重なる部分のみで計算を行います。
import numpy as np
# 入力配列 a
a = np.array([1, 2, 3, 4, 5])
# 入力配列 v
v = np.array([1, 2, 3])
# numpy.correlate() を使用して相関を計算
# mode='valid' はデフォルト値なので省略可能ですが、明示的に記述しています
correlation_valid = np.correlate(a, v, mode='valid')
print(f"配列 a: {a}")
print(f"配列 v: {v}")
print(f"mode='valid' での相関結果: {correlation_valid}")
# 出力:
# 配列 a: [1 2 3 4 5]
# 配列 v: [1 2 3]
# mode='valid' での相関結果: [10 16 22]
解説
mode='valid'
は、v
がa
に完全に収まる範囲でのみ計算されるため、出力の長さはmax(len(a), len(v)) - min(len(a), len(v)) + 1
となります(この例では5 - 3 + 1 = 3
)。correlation_valid
の結果[10, 16, 22]
は以下のように計算されます。v
を反転させると[3, 2, 1]
になります。a
の最初の3つの要素[1, 2, 3]
と[3, 2, 1]
の積和:(1*3) + (2*2) + (3*1) = 3 + 4 + 3 = 10
v
を1つ右にずらし、a
の[2, 3, 4]
と[3, 2, 1]
の積和:(2*3) + (3*2) + (4*1) = 6 + 6 + 4 = 16
v
をさらに1つ右にずらし、a
の[3, 4, 5]
と[3, 2, 1]
の積和:(3*3) + (4*2) + (5*1) = 9 + 8 + 5 = 22
例2: mode='full'
を使用した相関計算とラグの解釈
mode='full'
は、v
が a
のわずかにでも重なるすべての位置で計算を行います。結果の配列が最も長くなります。これは、信号の遅延を検出する際などに便利です。
import numpy as np
# 入力配列 a (信号)
a = np.array([0, 0, 1, 2, 3, 2, 1, 0, 0])
# 入力配列 v (検出したいパターン)
v = np.array([1, 2, 1])
# mode='full' で相関を計算
correlation_full = np.correlate(a, v, mode='full')
print(f"信号 a: {a}")
print(f"パターン v: {v}")
print(f"mode='full' での相関結果: {correlation_full}")
# 相関値が最大となるインデックスを見つける
max_correlation_index = np.argmax(correlation_full)
print(f"最大相関値のインデックス: {max_correlation_index}")
# ラグの計算(重要!)
# fullモードの場合、出力配列のインデックスから実際のラグを計算する必要があります。
# vがaの左端からスライドし始めると考えると、
# 結果のインデックス i は、vがaに対して (len(v) - 1 - i) だけ左にずれている状態を示します。
# したがって、ラグ(vがaに対してどれだけ右にずれているか)は i - (len(v) - 1) となります。
lag_at_max_correlation = max_correlation_index - (len(v) - 1)
print(f"最大相関となるラグ: {lag_at_max_correlation}")
# 出力例:
# 信号 a: [0 0 1 2 3 2 1 0 0]
# パターン v: [1 2 1]
# mode='full' での相関結果: [ 0 0 1 4 10 8 4 0 0 0 0]
# 最大相関値のインデックス: 4
# 最大相関となるラグ: 2
解説
- 計算されたラグ
2
は、パターンv
の開始が信号a
のインデックス2
に位置することを示しています。これはa
の[1, 2, 3]
の部分に対応します。 - このインデックス
4
は、パターンv
が信号a
の中でa[2], a[3], a[4]
(つまり[1, 2, 3]
) の位置と最も相関が高いことを示しています。 - 結果
[0, 0, 1, 4, 10, 8, 4, 0, 0, 0, 0]
の中で、最大値10
はインデックス4
にあります。 mode='full'
の出力長さはlen(a) + len(v) - 1
(9 + 3 - 1 = 11
) となります。
例3: mode='same'
を使用した相関計算
mode='same'
は、出力の長さが入力のうち長い方の配列と同じになります。結果の中心が相関の中心に対応します。
import numpy as np
# 入力配列 a
a = np.array([1, 2, 3, 4, 5])
# 入力配列 v
v = np.array([1, 2, 3])
# mode='same' で相関を計算
correlation_same = np.correlate(a, v, mode='same')
print(f"配列 a: {a}")
print(f"配列 v: {v}")
print(f"mode='same' での相関結果: {correlation_same}")
# 出力:
# 配列 a: [1 2 3 4 5]
# 配列 v: [1 2 3]
# mode='same' での相関結果: [ 4 10 16 22 20]
解説
- このモードでは、出力の中心がゼロラグに対応するようにパディングが行われます。ただし、結果の解釈は
full
モードに比べると直感的ではない場合があります。通常、ラグ検出にはfull
モードが推奨されます。 - 出力の長さは
len(a)
と同じ5
です。
周期的な信号とその遅延版の相関を計算し、遅延時間を推定します。
import numpy as np
import matplotlib.pyplot as plt # 結果を可視化するためにMatplotlibを使用
# サンプルレート
fs = 100 # Hz
t = np.arange(0, 1, 1/fs) # 1秒間の時間軸
# 基準信号(サイン波)
signal_a = np.sin(2 * np.pi * 5 * t) # 5Hzのサイン波
# 遅延させた信号(信号aを0.1秒遅延させる)
delay_samples = int(0.1 * fs) # 100Hzで0.1秒は10サンプル
signal_b = np.roll(signal_a, delay_samples) # 配列をシフトして遅延をシミュレート
# signal_b の最初の方は0で埋まるので、よりリアルなシミュレーションのためにはゼロパディングされた信号を生成するのが良い場合もありますが、
# ここでは単純化のために np.roll を使用
# 相互相関を計算 (fullモードで)
# np.correlate(signal_a, signal_b) と np.correlate(signal_b, signal_a) は結果が反転する点に注意
correlation = np.correlate(signal_a, signal_b, mode='full')
# 相関結果の可視化
plt.figure(figsize=(10, 6))
plt.subplot(3, 1, 1)
plt.plot(t, signal_a)
plt.title('Signal A (Reference)')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t, signal_b)
plt.title(f'Signal B (Delayed by {delay_samples} samples)')
plt.grid(True)
plt.subplot(3, 1, 3)
# 相関結果のX軸(ラグ)を計算
# fullモードの出力長さ
n = len(signal_a)
m = len(signal_b)
lags = np.arange(-(n - 1), m) # ラグの範囲は -(N-1) から M-1 まで
plt.plot(lags / fs, correlation) # ラグを秒単位で表示
plt.title('Cross-Correlation between A and B')
plt.xlabel('Lag (seconds)')
plt.ylabel('Correlation Value')
plt.grid(True)
# 最大相関値を見つける
max_corr_idx = np.argmax(correlation)
estimated_lag_samples = lags[max_corr_idx]
estimated_lag_time = estimated_lag_samples / fs
plt.axvline(estimated_lag_time, color='r', linestyle='--', label=f'Estimated Lag: {estimated_lag_time:.3f} s')
plt.legend()
plt.tight_layout()
plt.show()
print(f"\n実際の遅延サンプル数: {delay_samples}")
print(f"推定された遅延サンプル数: {estimated_lag_samples}")
print(f"推定された遅延時間: {estimated_lag_time:.3f} 秒")
# 出力例:
# 実際の遅延サンプル数: 10
# 推定された遅延サンプル数: -10 # ここが注意点: np.correlate(a, b) は b が a に対してどれだけ進んでいるかを示す
# 推定された遅延時間: -0.100 秒
numpy.correlate()
は相互相関を計算するのに便利ですが、パフォーマンス、機能、あるいは計算の概念的な違いから、他の方法を使用したい場合があります。ここでは、いくつかの主要な代替方法とその使い分けについて解説します。
scipy.signal.correlate() (推奨される代替手段)
科学技術計算ライブラリSciPyに含まれる scipy.signal.correlate()
は、numpy.correlate()
の直接的な代替として最もよく使われます。主な利点は、大規模なデータセットに対して高速なFFT(高速フーリエ変換)ベースのアルゴリズムを使用できる点です。
- 用途
- ほとんどの場合、
numpy.correlate()
の代わりにこちらを使うことを推奨します。特に、信号処理、画像処理、時系列分析などで長い信号や大きな画像を取り扱う場合。
- ほとんどの場合、
- 欠点
- NumPy だけでなく SciPy もインストールする必要がある。
- 利点
- 大規模な配列に対してはるかに高速(特に
mode='full'
の場合)。 method
引数で'auto'
,'direct'
,'fft'
を選択できる。boundary
とfillvalue
を指定して境界処理を細かく制御できる。normalize
引数で結果を自動的に正規化できる(ただし、相関係数のような厳密な正規化とは異なる場合がある)。
- 大規模な配列に対してはるかに高速(特に
コード例
import numpy as np
from scipy.signal import correlate
a = np.array([1, 2, 3, 4, 5])
v = np.array([1, 2, 3])
# numpy.correlate と同様の計算
scipy_corr_full = correlate(a, v, mode='full', method='auto')
print(f"scipy.signal.correlate (full): {scipy_corr_full}")
# 正規化された相関 (注意: 相関係数とは異なる)
scipy_corr_normalized = correlate(a, v, mode='valid', normalize=True)
print(f"scipy.signal.correlate (valid, normalized): {scipy_corr_normalized}")
numpy.convolve() を利用する (相互相関の定義から)
相互相関は、一方の信号を反転させてから畳み込みを行うことと同等です。この性質を利用して、numpy.convolve()
を相互相関の計算に使用できます。
- 用途
numpy.correlate()
の内部動作を理解したい場合。- SciPy をインストールできない環境で、かつデータサイズが小さい場合。
- 欠点
- 手動で一方の配列を反転させる必要がある(忘れやすい)。
numpy.correlate()
と同じく、大規模な配列ではFFTベースのアルゴリズムに比べて遅い。
- 利点
numpy
だけで完結する。- 相互相関の数学的な定義を理解するのに役立つ。
コード例
import numpy as np
a = np.array([1, 2, 3, 4, 5])
v = np.array([1, 2, 3])
# v を反転させる
v_reversed = v[::-1]
# 畳み込みを実行
# mode='full', 'same', 'valid' は numpy.correlate() と同じように指定可能
conv_corr_full = np.convolve(a, v_reversed, mode='full')
conv_corr_valid = np.convolve(a, v_reversed, mode='valid')
print(f"numpy.convolve (full) を使った相関: {conv_corr_full}")
print(f"numpy.convolve (valid) を使った相関: {conv_corr_valid}")
# numpy.correlate と比較
# print(f"numpy.correlate (full): {np.correlate(a, v, mode='full')}")
# print(f"numpy.correlate (valid): {np.correlate(a, v, mode='valid')}")
フーリエ変換 (FFT) を直接利用する
相互相関定理(Cross-correlation theorem)によれば、2つの信号の相互相関は、一方の信号のフーリエ変換と、もう一方の信号の共役複素フーリエ変換の積の逆フーリエ変換に等しいとされます。これにより、大規模なデータで非常に高速な計算が可能です。
- 用途
- 最高速の性能が求められる場合。
- 周波数領域での信号処理に慣れている場合。
scipy.signal.correlate
の内部動作を理解したい場合。
- 欠点
- パディング(ゼロパディング)や正規化の処理を自分で実装する必要がある。
- 概念的にやや複雑。
- 利点
- 非常に長い信号に対して究極的に高速(
scipy.signal.correlate
のFFT実装の基礎)。 - 相互相関の周波数領域での理解が深まる。
- 非常に長い信号に対して究極的に高速(
コード例
import numpy as np
a = np.array([1, 2, 3, 4, 5], dtype=float)
v = np.array([1, 2, 3], dtype=float)
# FFTに最適な長さを決定 (通常は2のべき乗にゼロパディング)
n = len(a) + len(v) - 1
n_fft = 2**int(np.ceil(np.log2(n)))
# 各配列をゼロパディングしてFFT
A_fft = np.fft.fft(a, n_fft)
V_fft = np.fft.fft(v, n_fft)
# 相互相関定理: XCORR(f) = A*(f) * CONJ(V(f))
# ただし、numpy.correlate の定義は V を反転させた畳み込みなので、
# V の共役複素数をとる代わりに V のフーリエ変換をそのまま使う(Vの共役の逆フーリエ変換に相当)
# または Aの共役 * V で計算することもできる (相互相関の別の定義)
# ここでは、np.correlate の結果に近づけるため、np.fft.fft(v[::-1], n_fft) と同じ効果を得る
# つまり、A_fft * np.conj(V_fft) で計算すると、v を反転させない相関になるため、
# np.correlate の動作に合わせるには A_fft * np.fft.fft(v[::-1], n_fft) と同じ効果を出す必要がある。
# 実際には、A_fft * np.conj(V_fft) は、a の相関をとったものを返すため、順序が逆になる。
# 正しいのは、np.fft.ifft(np.fft.fft(a, n_fft) * np.conj(np.fft.fft(v, n_fft))) の実部
# と、結果をシフトする。
# より簡潔に、しかし定義に従う方法
# 相互相関定理: R_xy(tau) = IFFT( FFT(x) * CONJ(FFT(y)) )
corr_fft = np.fft.ifft(A_fft * np.conj(V_fft))
# 結果は複素数なので実部を取り出す
corr_fft_real = np.real(corr_fft)
# 結果のシフト(ゼロラグが中央に来るように)
# np.correlate の 'full' モードと結果を合わせるには、特定の範囲を切り出す
result_full_mode = corr_fft_real[:n] # 必要な長さだけ切り出す
print(f"FFT を使った相関 (full mode に相当): {result_full_mode}")
# np.correlate と比較
# print(f"numpy.correlate (full): {np.correlate(a, v, mode='full')}")
解説
np.fft.ifft
の結果は通常、ゼロラグが配列の先頭に来るため、numpy.correlate
の'full'
モードの出力と直接比較するには、結果を循環的にシフトする必要があります。上記のコードでは、必要な長さだけを切り出すことで'full'
モードに似た結果を得ています。np.conj(V_fft)
がポイントで、これによりv
が反転された状態での畳み込み(つまり相互相関)が実現されます。- FFT を使った相互相関の計算は、
numpy.correlate
やscipy.signal.correlate
が内部でどう動いているかを理解する上で非常に重要です。
最も原始的な方法で、Python のループを使って相関を計算します。これは非常に遅いため、実際のアプリケーションではほとんど使いません。しかし、相互相関の計算原理を理解するには役立ちます。
- 用途
- 相互相関の概念を学習する初期段階。
- 欠点
- 非常に低速。特に長い配列では実用的ではない。
- 利点
- 計算の仕組みを最も明確に理解できる。
コード例
import numpy as np
def manual_correlate(a, v, mode='full'):
"""
手動で相互相関を計算する関数 (教育目的)
numpy.correlate と同じ結果を返すように調整
"""
na = len(a)
nv = len(v)
# v を反転
v_rev = v[::-1]
# 出力配列の長さを決定
if mode == 'full':
output_len = na + nv - 1
elif mode == 'valid':
output_len = max(0, na - nv + 1)
elif mode == 'same':
output_len = na # numpy.correlate の same モードは、a の長さに合わせる
else:
raise ValueError("Invalid mode")
result = np.zeros(output_len)
# 畳み込み(相関)の計算
for k in range(output_len):
sum_val = 0
# convolution と同じ計算だが、v_rev を使う
# output[k] = sum(a[i] * v_rev[k-i])
# これは、a[i+lag] * v_rev[i] に相当
for i in range(nv):
# 'full' モードのインデックス計算
# 畳み込み結果のk番目の要素に対応する a の開始インデックス
a_idx = k - (nv - 1) + i
if 0 <= a_idx < na:
sum_val += a[a_idx] * v[i] # ここでは v は反転せずに使っている
# np.correlate(a,v) は np.convolve(a, v[::-1]) と同じ
# よって sum(a[idx+k] * v[idx]) となるようにループを組む
# => 結果的に v を反転させないで sum(a_idx * v_idx) の形にする
# np.correlate(a, v) のk番目の要素 = sum_{j=0}^{len(v)-1} a[j+k] * v[j]
# という定義に基づくと、ループは下記になる
# np.correlate(a,v) の定義
# out[k] = sum_i a[i+k] * v[i] (kはラグ)
# v を反転させる必要はない。これは畳み込みの定義と異なるため注意が必要
# np.correlate は v を反転させずに定義されることもある
# 実際の np.correlate の結果は np.convolve(a, v[::-1]) と一致するため、
# v[::-1] を使って畳み込みを行う方が混乱が少ない
# 手動で numpy.correlate の動作を再現
# numpy.correlate(a,v) = np.convolve(a, v[::-1])
# なので、畳み込みとして計算し、v を反転させる
# v_rev を使った畳み込みの計算
if mode == 'full':
a_idx = k - (nv - 1) + i
if 0 <= a_idx < na:
sum_val += a[a_idx] * v_rev[i]
elif mode == 'valid':
a_idx = k + i
if 0 <= a_idx < na and 0 <= i < nv:
sum_val += a[a_idx] * v_rev[i]
elif mode == 'same':
# 'same'モードのインデックス調整はより複雑になる
# 一般的には 'full' で計算後、中央部分を切り出す
pass # この例では複雑になるため省略
result[k] = sum_val
# same モードの簡易実装 (full から切り出し)
if mode == 'same':
full_result = manual_correlate(a, v, mode='full')
start = (len(full_result) - output_len) // 2
end = start + output_len
return full_result[start:end]
return result
a = np.array([1, 2, 3, 4, 5])
v = np.array([1, 2, 3])
manual_corr_full = manual_correlate(a, v, mode='full')
manual_corr_valid = manual_correlate(a, v, mode='valid')
manual_corr_same = manual_correlate(a, v, mode='same')
print(f"手動 (full): {manual_corr_full}")
print(f"手動 (valid): {manual_corr_valid}")
print(f"手動 (same): {manual_corr_same}")
print(f"numpy.correlate (full): {np.correlate(a, v, mode='full')}")
print(f"numpy.correlate (valid): {np.correlate(a, v, mode='valid')}")
print(f"numpy.correlate (same): {np.correlate(a, v, mode='same')}")
注意点
上記の手動実装は教育目的であり、numpy.correlate
の正確な実装を完全に再現するのは複雑です(特に mode='same'
のパディングやオフセット)。実際には、ネイティブな NumPy 関数や SciPy 関数を使うべきです。
- 最高のパフォーマンスを追求し、FFTの原理を理解している場合
NumPy のfft
モジュールを直接使用して、ゼロパディングと結果のシフトを自分で管理します。 - 相互相関の定義を理解したい場合
numpy.convolve()
を使って手動で一方の配列を反転させる方法を試すと良いでしょう。 - SciPy を使いたくない、またはFFTのオーバーヘッドを避けたい小規模な配列の場合
numpy.correlate()
をそのまま使用します。 - ほとんどの場合
scipy.signal.correlate()
を使用するのが最も良い選択です。パフォーマンスと機能のバランスが優れています。