实现由 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 方法。我在对FIR filter in CUDA 的回答中使用 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();
}

Answer 2

对于另一种方法，您可以截断指数移动平均窗口，然后通过在信号和加窗指数之间进行卷积来计算过滤后的信号。可以使用免费的 CUDA FFT 库 (cuFFT) 计算卷积，因为您可能知道，卷积可以表示为傅里叶域中两个信号的逐点乘法（这就是卷积定理的恰当名称，它以 O(n log(n)) 的复杂度运行。这种方法将最小化您的 CUDA 内核代码并运行得非常非常快，即使在 GeForce 570 上也是如此；如果您可以以单（浮点）精度进行所有计算，则尤其如此。