次のコードを使用して、インプレースの前方後方 FIR フィルタリングを作成します。
lena = len(a)
lenb = len(b)
convol = zeros(a.shape)
code = """
// Forward convolution
int pad = (lenb-1)/2;
int i, j;
for (i=pad; i<pad+lena; i++)
{
int kmin, kmax, k;
// Reverse indexing for the next pass
j = lena-1-i+pad;
convol(j) = 0;
kmin = (i >= lenb - 1) ? i - (lenb - 1) : 0;
kmax = (i < lena - 1) ? i : lena - 1;
for (k = kmin; k <= kmax; k++)
{
convol(j) += a(k)*b(i - k);
}
}
// Backward convolution (the signal in convol has been
// reversed using reversed indexes)
for (i=pad; i<pad+lena; i++)
{
int kmin, kmax, k;
// Reverse indexing for reordering the output vector
j = lena-1-i+pad;
a(j) = 0;
kmin = (i >= lenb - 1) ? i - (lenb - 1) : 0;
kmax = (i < lena - 1) ? i : lena - 1;
for (k = kmin; k <= kmax; k++)
{
a(j) += convol(k)*b(i - k);
}
}
return_val = 1;
"""
weave.inline(code, [ 'a', 'b', 'lena', 'lenb', 'convol'],
type_converters=converters.blitz, compiler = 'g++')
もちろん、「convol」変数は C/C++ スコープの外にある必要はなく、スペースと処理の最適化の両方が必要です。したがって、このコードを次のように置き換えることは (少なくとも私にとっては) 理にかなっています。
lena = len(a)
lenb = len(b)
code = """
// Forward convolution
int pad = (lenb-1)/2;
int i, j;
float* convol = (float*) calloc(lena, sizeof(float));
for (i=pad; i<pad+lena; i++)
{
int kmin, kmax, k;
// Reverse indexing for the next pass
j = lena-1-i+pad;
convol[j] = 0;
kmin = (i >= lenb - 1) ? i - (lenb - 1) : 0;
kmax = (i < lena - 1) ? i : lena - 1;
for (k = kmin; k <= kmax; k++)
{
convol[j] += a(k)*b(i - k);
}
}
// Backward convolution (the signal in convol has been
// reversed using reversed indexes)
for (i=pad; i<pad+lena; i++)
{
int kmin, kmax, k;
// Reverse indexing for reordering the output vector
j = lena-1-i+pad;
a(j) = 0;
kmin = (i >= lenb - 1) ? i - (lenb - 1) : 0;
kmax = (i < lena - 1) ? i : lena - 1;
for (k = kmin; k <= kmax; k++)
{
a(j) += convol[k]*b(i - k);
}
}
free(convol);
return_val = 1;
"""
weave.inline(code, [ 'a', 'b', 'lena', 'lenb'],
type_converters=converters.blitz, compiler = 'g++')
唯一の違いは、numpy 配列を使用する代わりに、C float 配列を直接使用したことです。問題は、2 番目のコードの処理に 1 番目のコードよりも約 2 倍の時間がかかることです...なぜでしょうか? 2番目のバージョンに何か問題がありますか?早くなるはず…!