1

次のコードを使用して、インプレースの前方後方 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番目のバージョンに何か問題がありますか?早くなるはず…!

4

0 に答える 0