实现由CUDA中的差分方程描述的指数移动平均滤波器

Question

我目前正在尝试使用CUDA来评估代表指数移动平均滤波器的差分方程.滤波器由以下差分方程描述

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

其中alpha和beta常量定义为

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

如何操纵上述差分方程以获得该滤波器的系统响应？在GPU上实现此过滤器的有效方法是什么？

我正在研发GTX 570.

Answer 1

我建议操纵如下所示的上述差分方程,然后使用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,以alpha除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();
}