4

I am trying to implement an IIR filter.

I tried to implement the following filter but

FIR : y(n) = b0(x[n]) + ... +bM-1(x[n-M+1])

IIR : y(n) = {b0(x[n]) + ... +bM-1(x[n-M+1])} - {a1(y[n-1]) + ... +aN(y[n-N}

I am confused on how to implement y[n-1]......

Here is my code.

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *coeffa, *coeffb;
    double *inputp;
    double *outputp;
    int n,k;

    //filter length =7
    for (n = 0; n < length; n++) {
        coeffa = coeffs_A;
        coeffb = coeffs_B;
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) {
            bcc += (*coeffb++) * (*inputp--); //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 1; k < filterLength; k++) {
            acc += (*coeffa++) * (*output--); //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

I will not copy here the code for memove function for seek of brevity.

1

4 Answers 4

3

An IIR filter is recursive which means it takes the past values of the output (equivalently it could be expressed as infinite sequence of the input).

Suppose you have the filter

y[n] = b*x[n]-a*y[n-1]

Generally the first output is initialized to a given value, for example 0.

Suppose you want to filter the signal x of length N then you should do something like:

double out[N]; //output vector
double a = 0.3, b=0.5; //assign the coefficients a value
out[0] = 0; //initialize the first element

for (int i=1; i<N; i++)
{
    out[i] = b*x[i] -a*[i-1];  
}

In the case of your code I cannot know what you are doing in the line inputp = &insamp[filterLength - 1 + n]; and that might be a problem.

I will assume inputp is the signal you want to filter.

Another issue: you used filter filterLength to indicate the lengths of both the input elements and the output elements: that is not generally in a IIR filter.

The output elements from 0 to filterlength should be initialized somehow, suppose to 0. Then in the second loop you made the loop index starting from 1 but the coefficient array should start from 0.

(DOUBT: If the oldest element is y[5] how can the length be 7?)

Using index instead of dereferencing the array your code should be something like this:

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *inputp;
    int n,k;


    for (int ii=0; ii<filterLength; ii++)
    {
        output[ii] = 0;
    }

    //filter length =7
    for (n = 0; n < length; n++) {
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) 
        {
            bcc += coeffb[k] * inputp[filterLength-k-1]; //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 0; k < filterLength; k++) 
        {
            acc += coeffa[k] * output[filterLength-k-1]; //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

EDIT I think the two for being the same length can be merged together:

for (k = 0; k < filterLength; k++) 
{
    output[n] += (coeffb[k] * inputp[filterLength-k-1] - coeffa[k] * output[filterLength-k-1]); 
}
1

If you really want to do it with pointers:

void filter1(const double *b, const double *a, size_t filterLength, const double *in, double *out, size_t length) {
    const double a0 = a[0];
    const double *a_end = &a[filterLength-1];
    const double *out_start = out;
    a++;
    out--;
    size_t m;
    for (m = 0; m < length; m++) {
        const double *b_macc = b;
        const double *in_macc = in;
        const double *a_macc = a;
        const double *out_macc = out;
        double b_acc = (*in_macc--) * (*b_macc++);
        double a_acc = 0;
        while (a_macc <= a_end && out_macc >= out_start) {
            b_acc += (*in_macc--) * (*b_macc++);
            a_acc += (*out_macc--) * (*a_macc++);
        }
        *++out = (b_acc - a_acc) / a0;
        in++;
    }
}

I compared the result of this algorithm to MATLAB's filter function.

Note: You can get a large performance gain if you normalize your coefficients (i.e. a0 == 1). To do that, you just divide your a and b vectors by a0, and then you don't have to divide b_acc - a_acc by a0 on each iteration.

1
  • in_macc-- will result in access violations when m < filter_length
    – Zhu Mengxi
    Commented May 26, 2023 at 19:56
0

y[n-1] is just the result from the previous time step ex: output[n-1] (y[n-N] is for the Nth time step from the current one). So you need to index into the output array by the amount thats appropriate output[n-k] (?) (be careful to avoid indexing beyond the beginning of the array however, you'll need code to prevent this).

Also you are not indexing into your coeffa/coeffb coefficients properly (I think). To get the result you want you might need to do this *(coeffa++) to make sure you increment the pointer before you dereference it. I would urge you to do this by indexing into the array directly ex: coeff_A[k], however, as it is much easier to see what is going on.

Last side note: the operation (*output--) is moving the output pointer in a way that I think is unintentional. Index into the array instead.

0

Biggest thing is that afaik u never want to specify A and B separately because B should always equal (1-A). So, that being said an implementation of a basic IIR filter is:

class IIR_Filter
{
    IIR_Filter(double _alpha, double initialValue = 0) 
    { 
        alpha = _alpha; 
        lastVal = initialValue;
    }
    
    double processNextStep(double newVal)
    {
        lastVal = alpha * newVal + (1.0 - alpha) * lastVal;
        return lastVal;
    }
private:
    double alpha;
    double lastVal;
};

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.