2

現在、CUDA を使用して、指数移動平均フィルターを表す差分方程式を評価しようとしています。フィルタは次の差分方程式で表されます。

y[n] = y[n-1] * beta + alpha * x[n]

ここでalpha、 とbetaは次のように定義された定数です。

alpha = (2.0 / (1 + Period))
beta = 1 - alpha

上記の差分方程式を操作して、このフィルターのシステム応答を得るにはどうすればよいでしょうか? このフィルターを GPU に実装する効率的な方法は何でしょうか?

私はGTX 570に取り組んでいます。

4

2 に答える 2

5

上記の差分方程式を以下に示すように操作し、CUDA Thrust プリミティブを使用することを提案します。

差分方程式の操作 - 差分方程式の明示的な形式

簡単な代数によって、次のことがわかります。

y[1] = beta * y[0] + alpha * x[1]

y[2] = beta^2 * y[0] + alpha * beta * x[1] + alpha * x[2]

y[3] = beta^3 * y[0] + alpha * beta^2 * x[1] + alpha * beta * x[2] + alpha * x[3]

したがって、明示的な形式は次のとおりです。

y[n] / beta^n = y[0] + alpha * x[1] / beta + alpha * x[2] / beta^2 + ...

CUDA スラストの実装

上記の明示的な形式は、次の手順で実装できます。

  1. 入力シーケンスd_inputをalphaexcept forに初期化しd_input[0] = 1.ます。
  2. d_1_over_beta_to_the_nに等しいベクトルを定義し1, 1/beta, 1/beta^2, 1/beta^3, ...ます。
  3. d_input要素ごとにd_1_over_beta_to_the_n;を掛けます。
  4. inclusive_scanのシーケンスを取得するために を実行しy[n] / beta^nます。
  5. 上記のシーケンスを で割り1, 1/beta, 1/beta^2, 1/beta^3, ...ます。

編集

上記のアプローチは、線形時変 (LTV) システムに推奨できます。線形時不変 (LTI) システムの場合、Paul が言及した FFT アプローチをお勧めします。CUDA のFIR フィルターへの回答で CUDA Thrust と cuFFT を使用して、そのアプローチの例を提供しています。

完全なコード

#include <thrust/sequence.h>

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>

int main(void)
{
    int N = 20;

    // --- Filter parameters
    double alpha    = 2.7;
    double beta     = -0.3;

    // --- Defining and initializing the input vector on the device
    thrust::device_vector<double> d_input(N,alpha * 1.);
    d_input[0] = d_input[0]/alpha;

    // --- Defining the output vector on the device
    thrust::device_vector<double> d_output(d_input);

    // --- Defining the {1/beta^n} sequence
    thrust::device_vector<double> d_1_over_beta(N,1./beta);
    thrust::device_vector<double> d_1_over_beta_to_the_n(N,1./beta);
    thrust::device_vector<double> d_n(N);
    thrust::sequence(d_n.begin(), d_n.end());
    thrust::inclusive_scan(d_1_over_beta.begin(), d_1_over_beta.end(), d_1_over_beta_to_the_n.begin(), thrust::multiplies<double>());
    thrust::transform(d_1_over_beta_to_the_n.begin(), d_1_over_beta_to_the_n.end(), d_input.begin(), d_input.begin(), thrust::multiplies<double>());    
    thrust::inclusive_scan(d_input.begin(), d_input.end(), d_output.begin(), thrust::plus<double>());
    thrust::transform(d_output.begin(), d_output.end(), d_1_over_beta_to_the_n.begin(), d_output.begin(), thrust::divides<double>());

    for (int i=0; i<N; i++) {
        double val = d_output[i];
        printf("Device vector element number %i equal to %f\n",i,val);
    }

    // --- Defining and initializing the input vector on the host
    thrust::host_vector<double> h_input(N,1.);

    // --- Defining the output vector on the host
    thrust::host_vector<double> h_output(h_input);

    h_output[0] = h_input[0];
    for(int i=1; i<N; i++)
    {
        h_output[i] = h_input[i] * alpha + beta * h_output[i-1];
    }

    for (int i=0; i<N; i++) {
        double val = h_output[i];
        printf("Host vector element number %i equal to %f\n",i,val);
    }

    for (int i=0; i<N; i++) {
        double val = h_output[i] - d_output[i];
        printf("Difference between host and device vector element number %i equal to %f\n",i,val);
    }

    getchar();
}
于 2014-04-30T16:45:22.757 に答える
4

別のアプローチとして、指数移動平均ウィンドウを切り捨ててから、信号とウィンドウ処理された指数との間で畳み込みを行うことにより、フィルター処理された信号を計算できます。たたみ込みは、無料の CUDA FFT ライブラリ (cuFFT) を使用して計算できます。これは、ご存知かもしれませんが、たたみ込みは、フーリエ ドメインでの 2 つの信号の点ごとの乗算として表現できるためです (これは、適切な名前の Convolution Theorem です。 O(n log(n)) の複雑さで実行されます)。このタイプのアプローチにより、CUDA カーネル コードが最小限に抑えられ、GeForce 570 でも非常に高速に実行されます。すべての計算を単精度 (float) で実行できる場合は特にそうです。

于 2014-04-30T17:04:19.127 に答える