【C++】【Eigen】NumPyのconvolve()関数(っぽいもの)を実装する
先日C++でプログラムを作成していたところ、PythonのNumPyのconvolve()のようなものが必要になったので、Eigenを使って同様の関数を実装してみました。
NumPyのconvolve()関数とは
NumPyのconvolve()関数は、畳み込み(コンボリューション)を行う関数です。
畳み込みは2つの異なる関数から新しい関数を作り出す操作です。数学的には以下のように定義されます。
$$
(f*g)(t) = \int^\infty_{-\infty}f(\tau)g(t-\tau)\rm{d}\tau
$$
積分を用いて定義されていますが、離散値においては以下のように置き換えられます。
$$
(f*g)[n] = \sum^\infty_{m=-\infty}f[m]g[n-m]
$$
ここでは畳み込みについて詳細は説明しません。検索すると多数の説明ページが見つかります。NumPyのconvolve()関数については以下をご覧ください。
numpy.convolve — NumPy v2.4 Manual
このリファレンスによると、境界データの計算の仕方によって3種類のモードがあります。
| モード | 説明 |
|---|---|
"full" | デフォルト。1点でも重なる部分は計算を行う。 2つの入力信号の長さをそれぞれ$N, M$としたとき、長さ$N+M-1$の出力信号を返す |
"same" | 入力信号と同じ長さとなるように計算を行う。 長さ$\max(N, M)$の出力信号を返す |
"valid" | 2つの信号が完全に重なりあう部分のみ計算を行う。 長さ$\max(N, M) – \min(N, M) + 1$の出力信号を返す。 |
convolve()関数の各モードが実際にどのような計算をしているかは、以下のページが詳しいです。
Eigenによる実装
Eigenを使い、C++でこれと似たような関数を実装します。"full", "same", "valid"の3つのモードにも対応させます。3つのモードの違いは、両端のデータをどのように処理するかの違いなのですが、今回は"full"相当の処理を行い、後からデータを抜き出すという方法で対応させました。
なお、畳み込み
$$
(f*g)[n] = \sum_{i}f[i]g[n-i]
$$
の計算は、「$n$番目の出力信号は、2つの入力信号の添字の和が$n$となるものを足し合わせたもの」という特徴があります。これに気付くと、処理が簡単に書けます。
また、計算モードの指定方法には、enum classを使いました。enum classの使い方については、こちらのエントリをご覧ください。
ソースコードはGitHubにも置いてあります。Eigenと同じようにヘッダファイルのインクルードのみで動作します。
#include <algorithm>
#include <Eigen/Eigen>
// 計算モード
enum class ConvolveMode
{
Full, Same, Valid
};
// NumPyのconvolve()関数っぽいもの
//
// v1 : 1次元ベクトル(VectorX<Scalar>またはRowVectorX<Scalar>)
// v2 : 1次元ベクトル(VectorX<Scalar>またはRowVectorX<Scalar>)
// mode: 計算モード(Full, Same, Valid)
template <typename Derived1, typename Derived2>
const typename Derived1::PlainObject
Convolve(
const Eigen::MatrixBase<Derived1> & v1,
const Eigen::MatrixBase<Derived2> & v2,
const ConvolveMode mode = ConvolveMode::Full
)
{
// MatrixBase<Derived>で受けているが、有効なクラスはベクトルのみとする
// それ以外はコンパイル時に弾く
static_assert(Derived1::IsVectorAtCompileTime, "Error: v1 is not vector.");
static_assert(Derived2::IsVectorAtCompileTime, "Error: v2 is not vector.");
int length, start_index, max_size, min_size;
typename Derived1::PlainObject result;
max_size = std::max(v1.size(), v2.size());
min_size = std::min(v1.size(), v2.size());
length = v1.size() + v2.size() - 1;
start_index = 0;
result = Derived1::PlainObject::Zero(length);
for(int j = 0; j < v2.size(); j++) {
for(int i = 0; i < v1.size(); i++) {
result(i + j) += v1(i) * v2(j);
}
}
if (mode == ConvolveMode::Valid) {
length = max_size - min_size + 1;
start_index = min_size - 1;
}
else if(mode == ConvolveMode::Same) {
length = max_size;
start_index = (int)(min_size / 2.0 + 0.5) - 1;
}
if(result.rows() == 1) {
return result.block(0, start_index, 1, length);
}
else {
return result.block(start_index, 0, length, 1);
}
}使い方
#include <iostream>
#include "convolve.h"
int main()
{
Eigen::RowVectorXd a(11), b(5);
a << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10;
b << 0.2, 0.3, 0.5, 0.8, 1;
std::cout << Convolve(a, b) << std::endl;
std::cout << Convolve(a, b, ConvolveMode::Same) << std::endl;
std::cout << Convolve(a, b, ConvolveMode::Valid) << std::endl;
return 0;
}出力結果
0 0.2 0.7 1.7 3.5 6.3 9.1 11.9 14.7 17.5 20.3 20.9 20.2 17 10
0.7 1.7 3.5 6.3 9.1 11.9 14.7 17.5 20.3 20.9 20.2
3.5 6.3 9.1 11.9 14.7 17.5 20.3- Pythonの出力結果
>>> import numpy as np
>>> np.convolve([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [0.2, 0.3, 0.5, 0.8, 1])
array([ 0. , 0.2, 0.7, 1.7, 3.5, 6.3, 9.1, 11.9, 14.7, 17.5, 20.3,
20.9, 20.2, 17. , 10. ])
>>> np.convolve([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [0.2, 0.3, 0.5, 0.8, 1], "same")
array([ 0.7, 1.7, 3.5, 6.3, 9.1, 11.9, 14.7, 17.5, 20.3, 20.9, 20.2])
>>> np.convolve([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [0.2, 0.3, 0.5, 0.8, 1], "valid")
array([ 3.5, 6.3, 9.1, 11.9, 14.7, 17.5, 20.3])Pythonの結果と出力が一致していることがわかります。
動作を確認した環境
- OS
- Windows 11
- Microsoft C/C++ Optimizing Compiler 19.44.35213 for x64
- MinGW-W64 GCC version 13.2.0
- Clang version 20.1.0
- Linux
- GCC version 12.2.0
- Windows 11
- Eigen version 3.4.0
まとめ
Eigenを使用してNumPyのconvolve()関数のようなものを自作してみました。2つのベクトルオブジェクト(VectorX<Scalar>, RowVectorX<Scalar>)に対して、畳み込みを行った結果を1番目のベクトルと同じ型で返します。



ディスカッション
コメント一覧
まだ、コメントがありません