平滑算法

Question

我编写了这段代码用于曲线平滑处理。 它会取一个点旁边的5个点，将它们相加并求平均值。

/* Smoothing */
void smoothing(vector<Point2D> &a)
{
    //How many neighbours to smooth
    int NO_OF_NEIGHBOURS=10;
    vector<Point2D> tmp=a;
    for(int i=0;i<a.size();i++)
    {

        if(i+NO_OF_NEIGHBOURS+1<a.size())
        {
            for(int j=1;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x+=a.at(i+j).x;
                a.at(i).y+=a.at(i+j).y;
            }
            a.at(i).x/=NO_OF_NEIGHBOURS;
            a.at(i).y/=NO_OF_NEIGHBOURS;

        }
        else
        {
            for(int j=1;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x+=tmp.at(i-j).x;
                a.at(i).y+=tmp.at(i-j).y;
            }
            a.at(i).x/=NO_OF_NEIGHBOURS;
            a.at(i).y/=NO_OF_NEIGHBOURS;
        }

    }

}

但是我得到了每个点的非常高的值，而不是类似于上一个点的相似值。形状被大大地最大化了，这个算法出了什么问题？

- rajat

---

使用 += 和 /= 可以使您的代码更易读（并且减少对 at(i) 的调用）。 - Nobody moving away from SE

如果您为它提供可预测结果的简单数据，并对其进行调试，会发生什么？ - Martin James

j是从0开始计数的，因此在您的内部循环中，您不是使用]i，i+ NO_OF_NEIGHBOURS]来平滑，而是使用[i，i+ NO_OF_NEIGHBOURS]。 - Tom Knapen

请注意，在else子句中，您正在访问已被平滑算法修改的元素！还要注意，将一个点与其邻居平均时，必须将总和除以“1 + NO_OF_NEIGHBORS”，因为该点本身包含在此总和中。 - Nobody moving away from SE

@Nobody 我知道我正在访问已经更改的点，有没有更好的替代方法来包含边界条件？ - rajat

Answer 1

这里看起来是一个反向实现有限脉冲响应（FIR）滤波器的例子，它实现了一个矩形窗口函数。从DSP的角度思考问题，您需要使用NO_OF_NEIGHBOURS个相等的FIR系数过滤输入的vector，每个系数的值为1/NO_OF_NEIGHBOURS。通常最好使用已有的算法而不是重新发明轮子。

这是我快速编写的用于过滤双精度浮点数的简陋实现。您可以轻松修改它以过滤您的数据类型。演示仅显示了几个上升锯齿函数（0，.25，.5，1）的周期，仅供演示目的。它可以编译，所以您可以随意尝试。

#include <iostream>
#include <vector>

using namespace std;

class boxFIR
{
    int numCoeffs; //MUST be > 0
    vector<double> b; //Filter coefficients
    vector<double> m; //Filter memories

public:
    boxFIR(int _numCoeffs) :
    numCoeffs(_numCoeffs)
    {
        if (numCoeffs<1)
            numCoeffs = 1; //Must be > 0 or bad stuff happens

        double val = 1./numCoeffs;
        for (int ii=0; ii<numCoeffs; ++ii) {
            b.push_back(val);
            m.push_back(0.);
        }
    }    

    void filter(vector<double> &a)
    {
        double output;

        for (int nn=0; nn<a.size(); ++nn)
        {
            //Apply smoothing filter to signal
            output = 0;
            m[0] = a[nn];
            for (int ii=0; ii<numCoeffs; ++ii) {
                output+=b[ii]*m[ii];
            }

            //Reshuffle memories
            for (int ii = numCoeffs-1; ii!=0; --ii) {
                m[ii] = m[ii-1];
            }                        
            a[nn] = output;
        }
    }


};

int main(int argc, const char * argv[])
{
    boxFIR box(1); //If this is 1, then no filtering happens, use bigger ints for more smoothing

    //Make a rising saw function for demo
    vector<double> a;
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);

    box.filter(a);

    for (int nn=0; nn<a.size(); ++nn)
    {
        cout << a[nn] << endl;
    }
}

使用此行增加滤波器系数，可以看到逐渐平滑的输出。仅使用1个滤波器系数时，没有平滑效果。

boxFIR box(1);

这段代码足够灵活，如果您想要改变窗口形状，可以通过修改构造函数中定义的系数来实现。

注意：由于这是一个因果滤波器（仅依赖于当前样本和之前的样本），因此这将与您的实现略有不同。您的实现不是因果的，因为它向未来的样本望去以进行平均，并且这就是为什么您需要条件语句来处理接近向量末尾的情况。如果您想使用此算法获得与您正在尝试使用滤波器获得的输出类似的输出，请将您的向量反向通过此算法（只要窗口函数对称即可）。这样，您就可以在没有算法中不良条件部分的情况下获得类似的输出。

- learnvst

---

我对FIR滤波器一无所知，我只是想平滑这些点。时域实现的FIR滤波器比这个算法更好吗？ - rajat

@rajat 绝对没错，过滤函数中间的 if 语句很讨厌，因为你必须为边界编写单独的算法。请查看我的编辑。 - learnvst

YEY。有人给了这个+1。这是最佳解决方案。咕哝咕哝咕哝。 - learnvst

反向传递可能有助于防止相位失真（在这种情况下它是零相位失真滤波器）。 - trig-ger

Answer 2

过滤对于“内存”平滑很有用。这是learnvst答案的反向传递，以防止相位失真：

for (int i = a.size(); i > 0; --i)
{
    // Apply smoothing filter to signal
    output = 0;
    m[m.size() - 1] = a[i - 1];

    for (int j = numCoeffs; j > 0; --j) 
        output += b[j - 1] * m[j - 1];

    // Reshuffle memories
    for (int j = 0; j != numCoeffs; ++j) 
        m[j] = m[j + 1];

    a[i - 1] = output;
}

关于MATLAB中的零相位失真FIR滤波器，更多信息请参考：http://www.mathworks.com/help/signal/ref/filtfilt.html

- trig-ger

Answer 3

在以下代码块中：

            for(int j=0;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x=a.at(i).x+a.at(i+j).x;
                a.at(i).y=a.at(i).y+a.at(i+j).y;
            }

对于每个邻居，你需要分别将a.at(i)的x和y添加到邻居的值中。

如果我理解正确，应该是这样的。

            for(int j=0;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x += a.at(i+j+1).x
                a.at(i).y += a.at(i+j+1).y
            }

- sardok

---

使用此更改时，“NO_OF_NEIGHBOURS”将比实际使用的邻居数量多一个。 - Nobody moving away from SE

正确的写法是，条件应该是j<NO_OF_NEIGHBOURS+1，为了更清晰，j应该从0开始，a.at(i+j) 应该改为 a.at(i+j+1)。 - sardok

变量名为NO_OF_NEIGHBOURS，表示要平滑的邻居数量。根据您更正的循环，j将有NO_OF_NEIGHBOURS - 1个不同的值，这意味着将平滑NO_OF_NEIGHBOURS - 1个邻居，这与名称有点相反直觉。 - Nobody moving away from SE

Answer 4

当您需要获取相邻点时，可以使用点本身进行加法运算 - 只需将索引偏移1：

for(int j=0;j<NO_OF_NEIGHBOURS;j++)
 {
    a.at(i).x += a.at(i+j+1).x
    a.at(i).y += a.at(i+j+1).y
 }

- Denis Ermolin

Answer 5

当前点的值被使用了两次：一次是因为你使用了+=，一次是如果y==0。所以你正在构建6个点的总和，但只除以5。这个问题在IF和ELSE情况下都存在。另外：你应该检查向量是否足够长，否则你的ELSE情况将读取负索引。

以下本身不是问题，只是一个想法：你考虑过使用一个算法，只触摸每个点两次吗？：你可以存储一个临时的x-y值（初始化为与第一个点相同），然后当你访问每个点时，你只需添加新点并减去最旧的点，如果它比你的NEIGHBOURS后面更远。你为每个点更新这个“运行总和”，并将这个值除以NEIGHBOURS数存储到新点中。

- Bernd Elkemann

---

谢谢你提供检查向量长度的提示。我很难理解备选算法。你能用伪代码解释一下吗？ - rajat

好的，明白了，谢谢 :)。时间复杂度在我的问题中并不需要担心，因为我有一个只有128个点的向量。 - rajat

Answer 6

这对我来说很好用：

for (i = 0; i < lenInput; i++)
{
    float x = 0;
    for (int j = -neighbours; j <= neighbours; j++)
    {
        x += input[(i + j <= 0) || (i + j >= lenInput) ? i : i + j];
    }
    output[i] = x / (neighbours * 2 + 1);
}

- SinOfficial