实现由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推力元。

差分方程操作 - 差分方程

通过简单的代数显式的形式，你可以找到以下内容：

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 THRUST IMPLEMENTATION

可以实现通过下列步骤在上述明确的形式：

1. 初始化的输入序列d_input到除了d_input[0] = 1.alpha;
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, ...。

EDIT

上述方法可推荐用于线性时变（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上也是如此;特别是如果你能以单一（浮点）精度完成所有的计算。