ベッセルフィルタの実装を最適化する

Question

私は博士論文の研究に使用する信号処理コードを開発しています。コードのプロファイリングと最適化を開始する段階にあります。私の最後の大きな最適化（Optimizing Disk IO）以来、コードのランタイムの50％を取っている新しいボトルネックは、私のデジタルローパスベッセルフィルタの実装です。

私が探しているのは、この計算を高速化する方法の提案です。

コードは何が起こっているの内訳と説明とともに、以下の通りです：すべての

void filter_signal(double *signal, bessel *lpfilter, int64_t length) 
{ 
    int64_t i; 
    int64_t p; 
    int64_t end; 
    int64_t order = lpfilter->order; 
    int64_t padding = lpfilter->padding; 
    double *paddedsignal = lpfilter->paddedsignal; 
    double *temp = lpfilter->temp; 
    double *dcof = lpfilter->dcof; 
    double *ccof = lpfilter->ccof; 
    end = length+2*(order+padding); 
    int64_t imax = order+padding; 
    double padval = signal_average(signal,padding); 

    memcpy(&paddedsignal[imax],signal,length*sizeof(double)); 

    for (i=0; i<imax; i++) 
    { 
     temp[i] = padval; 
     paddedsignal[i] = padval; 
     paddedsignal[end-1-i] = padval; 
     temp[end-1-i] = padval; 
    } 
    for (i=order; i<end; i++) 
    { 
     temp[i] = ccof[0]*paddedsignal[i]; 
     for (p=1; p<=order; p++) 
     { 
      temp[i] += ccof[p]*paddedsignal[i-p] - dcof[p]*temp[i-p]; 
     } 
    } 
    padval = signal_average(&temp[order],padding); 
    for (i=0; i<imax; i++) 
    { 
     paddedsignal[end-1-i] = padval; 
     paddedsignal[i] = padval; 
    } 
    for (i=order; i<end; i++) 
    { 
     paddedsignal[end-1-i] = ccof[0]*temp[end-1-i]; 
     for (p=1; p<=order; p++) 
     { 
      paddedsignal[end-1-i] += ccof[p]*temp[end-1-i+p] - dcof[p]*paddedsignal[end-1-i+p]; 
     } 
    } 
    memcpy(signal,&paddedsignal[order+padding],length*sizeof(double)); 
} 

まず：signal配列は最大length = 1E7エントリを言う（非常に大きく、私は多くを処理するかもしれません何千ものこれらの配列が1回の実行で）、私はランタイムの多くがデータをキャッシュにロードするだけで、何らかの利益が得られるかもしれないと推測しています。ベッセルフィルタリングは次のように機能します：係数の配列（dcofおよびccof）の長さがそれぞれorderで、2〜10の間になります。各点のフィルタ処理された信号は、配列内の前の点と前の点の加重和ですすでにフィルタリングされた配列では、ccofとdcofが重み付けされています。

2つの複雑さがあります：1つは有限長配列のため、このようなフィルタリングはエッジ効果を導入します。これを回避する方法は、配列を平均でパディングし、フィルタリング後にパディングを破棄して、実際のデータが開始されるまでにエッジ効果がなくなるようにすることです。第2の複雑さは、フィルタリングがデータに位相シフトを導入することである（フィルタリングされたアレイは、元のアレイからのいくつかのサンプル数によってオフセットされる）。これを回避する方法は、ノイズの高周波成分を除去してデータを位相シフトした後、周波数成分にはほとんど影響を与えないが、位相シフトを逆転させるフォワード（forward）を2回フィルタリングする方法です。これらの修正は両方とも以下に実装されています。もう少し詳細にコードをステップ

：

memcpy(&paddedsignal[imax],signal,length*sizeof(double)); 

for (i=0; i<imax; i++) 
{ 
    temp[i] = padval; 
    paddedsignal[i] = padval; 
    paddedsignal[end-1-i] = padval; 
    temp[end-1-i] = padval; 
} 

paddedsignalは、中央部にsignalを保持しているが、エッジ効果を回避するために、order+padding試料と両端に埋め込まれた一時的なアレイです。 tempは、paddedsignalと同じディメンションを持つ一時的な配列です。これは、前方および後方のフィルタリングを行うことができないため必要です。両方のパッディング部分は、最初の数サンプルの平均で満たされ、エラーを低減します。

for (i=order; i<end; i++) 
{ 
    temp[i] = ccof[0]*paddedsignal[i]; 
    for (p=1; p<=order; p++) 
    { 
     temp[i] += ccof[p]*paddedsignal[i-p] - dcof[p]*temp[i-p]; 
    } 
} 

これはフォワードフィルタリングループです。それが終わると、tempには、パディングされた順方向フィルタリングされ位相シフトされた信号が含まれます。

for (i=order; i<end; i++) 
{ 
    paddedsignal[end-1-i] = ccof[0]*temp[end-1-i]; 
    for (p=1; p<=order; p++) 
    { 
     paddedsignal[end-1-i] += ccof[p]*temp[end-1-i+p] - dcof[p]*paddedsignal[end-1-i+p]; 
    } 
} 

位相シフトを元に戻す逆方向フィルタリングループです。それが完了すると、パディング信号には、パディングされた、フィルタされた、および位相シフトされていないデータが含まれます。次に、中央部分を信号配列に戻して格納し、パディングを削除します。

特に、キャッシュミスを避けるためにこれを最適化するクリーンな方法があるのだろうかと思います。実際のフィルタリングループでは、サンプルiでのフィルタリングされた配列の値は、両方ともに依存しています。フィルタリングされていない配列の値はi-pになります。これがキャッシュ非友好的になるのではないかと思いますか？また、2番目のループは配列を後方に移動します。それは問題なのですか？

Answer 1

最初と最後にmemcpyコールを削除することをお勧めします。これらはおそらく高価な操作であり、元のアレイ上で作業の大部分を「インプレース」で実行できるようです。私はパディングに関するあなたのコメントを見ました - それで、メソッドを呼び出す前に元の配列にパディングを作成するか、この関数のローカル配列にのみパッドの値を格納するロジックを追加することで、それを処理することを検討します。できるだけコピーを避けてください。