优化贝塞尔滤波器实现

Question

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

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

代码如下，以及对所发生情况的细分和解释：

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 个条目，我可能会在一次运行中处理数千个这样的数组），所以我猜很多运行时只是将数据加载到缓存中，这可能会有所收获。贝塞尔滤波的工作原理如下：我们有一些系数数组（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 是一个与paddedsignal 具有相同维度的临时数组，这是必要的，因为不可能就地进行前向和后向过滤。两者的填充部分都填充了前几个样本的平均值，从而减少了错误。

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

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

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

Answer 1

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