问Matlab滤波器的实现
EN

Stack Overflow用户

提问于 2011-12-12 13:11:37

回答 2查看 2.6K关注 0票数 2

我读到过Matlab filter命令是用来解差分方程的。filter()是在内部使用z变换，还是简单地使用递归，也就是说，在初始值x(0)，y(0)的情况下，它只是在时间上向前运行差分方程？如果这个问题没有意义，很抱歉，我是这个领域的新手。

谢谢,

matlab

signal-processing

回答 2

Stack Overflow用户

发布于 2011-12-12 14:46:49

iir ()函数实现MatLab滤波器(递归滤波器)，即它求解差分方程。在频域中的实现将具有更高的成本。时域为O(N)，频域理想为log(N)，如果使用FFT，则为O(N²)。

票数 2

Stack Overflow用户

发布于 2014-12-03 11:32:47

这是我前段时间如何用C++实现MATLAB的内置filter函数，如果有人需要它的话。在Zi中，您可以传递初始条件，并在需要时收集最终条件。

#include <vector>
#include <exception>
#include <algorithm>

typedef vector<double> vectord;

void filter(vectord B, vectord A, const vectord &X, vectord &Y, vectord &Zi)
{
    if (A.empty())
        throw std::domain_error("The feedback filter coefficients are empty.");
    if (std::all_of(A.begin(), A.end(), [](double coef){ return coef == 0; }))
        throw std::domain_error("At least one of the feedback filter coefficients has to be non-zero.");
    if (A[0] == 0)
        throw std::domain_error("First feedback coefficient has to be non-zero.");

    // Normalize feedback coefficients if a[0] != 1;
    auto a0 = A[0];
    if (a0 != 1.0)
    {       
        std::transform(A.begin(), A.end(), A.begin(), [a0](double v) { return v / a0; });
        std::transform(B.begin(), B.end(), B.begin(), [a0](double v) { return v / a0; });
    }

    size_t input_size = X.size();
    size_t filter_order = std::max(A.size(), B.size());
    B.resize(filter_order, 0);
    A.resize(filter_order, 0);  
    Zi.resize(filter_order, 0);
    Y.resize(input_size);

    for (size_t i = 0; i < input_size; ++i)
    {
        size_t order = filter_order - 1;
        while (order)
        {
            if (i >= order)
                Zi[order - 1] = B[order] * X[i - order] - A[order] * Y[i - order] + Zi[order];
            --order;
        }
        Y[i] = B[0] * X[i] + Zi[0];
    }
    Zi.resize(filter_order - 1);
}

您可以使用以下代码对其进行测试：

TEST_METHOD(TestFilter)
{
    vectord b_coeff = { /* Initialise your feedforward coefficients here */ };
    vectord a_coeff = { /* Initialise your feedback coefficients here */ };

    vectord input_signal = { /* Some input data to be filtered */ };
    vectord y_filter_ori = { /* MATLAB output from calling y_filter_ori = filter(b_coeff, a_coeff, input_signal); */ };

    vectord y_filter_out; vectord zi = { 0 };  // Set empty initial conditions
    filter(b_coeff, a_coeff, input_signal, y_filter_out, zi);
    Assert::IsTrue(compare(y_filter_out, y_filter_ori, 0.0001));
}

bool compare(const vectord &original, const vectord &expected, double tolerance = DBL_EPSILON)
{
    if (original.size() != expected.size())
        return false;
    size_t size = original.size();

    for (size_t i = 0; i < size; ++i)
    {
        auto diff = abs(original[i] - expected[i]);
        if (diff >= tolerance)
            return false;
    }
    return true;
}