CODEWITHSUNDEEP
matlabsignal-processing
BBSysDyn
BBSysDyn

Reputation: 4601

Matlab Filter Implementation

I read that Matlab filter command is used to solve difference equations. Does filter() internally use z-Transform or does it simply use recursion, that is, with a starting value x(0),y(0), it simply runs the difference equations forward in time? Sorry if this question does not make sense, I am new in this field.

Thanks,

Upvotes: 2

Views: 2844

Answers (3)

Darien Pardinas
Darien Pardinas

Reputation: 6186

Here is how I implemented MATLAB's built-in filter funtion in C++ some time ago, if someone needs it. In Zi you pass the initial conditions and collect the final conditions if needed. 

#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);
}

Feel free to test it with this code:

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

Upvotes: 1

Andr&#233; Bergner
Andr&#233; Bergner

Reputation: 1439

The MatLab filter() function implements iir filters (recursive filters), i.e. it solves difference equations. An implementation in frequency domain would have a much higher cost. Time-domain is O(N), frequency domain ideally log(N) if one would use FFT, and O(N²).

Upvotes: 2

Phonon
Phonon

Reputation: 12727

Filter implementation is likely some clever usage of frequency-domain techniques, however according to MATLAB documentation, it is mathematically equivalent to solving a difference equation:

 Y = FILTER(B,A,X) filters the data in vector X with the
    filter described by vectors A and B to create the filtered
    data Y.  The filter is a "Direct Form II Transposed"
    implementation of the standard difference equation:

    a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

The initial conditions are picked to be all zero, since we assume simply absence of signal (all zeros) before we start filtering. If you would like to specify those initial conditions, the filter command allows you to specify vectors of initial (Zi) and final (Zf) conditions:[Y,Zf] = FILTER(B,A,X,Zi)

Upvotes: 3

Related Questions