c - 优化贝塞尔滤波器的实现

Question

我正在开发一段用于我的博士研究的信号处理代码，并且我正处于开始分析和优化我的代码的阶段。自从我上次进行重大优化 ( Optimizing Disk IO ) 以来，新的瓶颈是我对数字低通贝塞尔滤波器的实现，占用了代码运行时间的 50%。

我正在寻找有关加快计算速度的方法的建议。

代码如下，以及所发生情况的 segmentation 和解释:

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个条目，我可能会在一次运行中处理数千个这样的数组) ，所以我猜测很多运行时只是花费在将数据加载到缓存中，这可能会带来一些 yield 。贝塞尔滤波的工作原理如下:我们有一些系数数组(dcof 和 ccof)，每个系数的长度为 order，介于 2 和10. 每个点处的滤波信号是数组中先前点和已滤波数组中先前点的加权和，其中 ccof 和 dcof 是权重。

有两个复杂问题:一是对于有限长度数组，这样的过滤会引入边缘效应。解决这个问题的方法是用平均值填充数组，并在过滤后丢弃填充，以便在实际数据开始时边缘效应消失。第二个复杂之处是滤波会在数据中引入相移(滤波后的阵列将偏移原始阵列中的一些样本)。解决这个问题的方法是进行两次滤波:一次向前，消除噪声的高频分量并对数据进行相移，然后再次向后，对频率分量几乎没有影响，但反转相移。这两个修复均在下面实现。

更详细地单步执行代码:

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 是一个临时数组，其尺寸与 patedsignal 相同，这是必需的，因为不可能就地进行前向和后向过滤。两者的填充部分都填充了前几个样本的平均值，从而减少了误差。

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];
    }
}

这是后向滤波循环，可消除相移。完成后，patedsignal 将包含填充的、过滤的和未相移的数据。然后，我们将中心部分存储回信号数组并删除填充。

特别是，我想知道是否有任何干净的方法来优化它以避免缓存丢失。有几个地方可能会感兴趣，但我希望在我花太多时间找出错误的树之前得到更有经验的程序员的意见:在实际的过滤循环中，示例 i< 处过滤数组的值 取决于它自己的值和 i-p 处未过滤数组的值。我想知道这是否会对缓存不友好？此外，第二个循环向后遍历数组。就速度而言，这是一个问题吗？

Answer 1

我建议尝试消除开头和结尾的 memcpy 调用。这些可能是昂贵的操作，似乎您可以在原始阵列上“就地”完成大部分工作。我看到了您关于填充的评论 - 因此我会考虑通过在调用该方法之前将填充构建到原始数组中或添加一些逻辑以仅将填充值存储在该函数本地的数组中来处理该问题。一般来说，尽可能避免复制。