Implementation of versatile IIR digital filter in C++

Question

I have realized that all the digital filters of the IIR type have the same structure. They are described by difference equation in following form:

$$ y(k) = b_0\cdot x(k) + b_1\cdot x(k-1) + \ldots + b_M\cdot x(k-M) -a_1\cdot y(k-1) - a_2\cdot y(k-2) - \ldots - a_N\cdot y(k-N) $$

Based on that I have attempted to define a C++ template class which can be used for implementation of any IIR filter based on passing values of the coefficients of the difference equation comming from design in Matlab or Scilab software package. The class is intended to be used in digital signal processor (with floating point ALU) for filtering in real time.

#ifndef IIRFILTER_H
#define IIRFILTER_H

#include <cstdint>
#include <iostream>

/**
 * @brief Versatile digital filter with infinite impulse response
 * i.e. difference equation in the form:
 * 
 * $\f 
 * y(k) =  b_0\cdot x(k) + b_1\cdot x(k-1) + \ldots + b_M\cdot x(k-M) 
 *        -a_1\cdot y(k-1) - a_2\cdot y(k-2) - \ldots - a_N\cdot y(k-N)
 * $\f
 * 
 * where $\f x(k) \ldots x(k-M)$\f are the input samples,
 * $\f y(k) \ldots y(k-N)$\f are the output samples,
 * $\f b_0, b_1, \ldots , b_M $\f are the input coefficients and
 * $\f a_1, a_2, \ldots , a_N $\f are the output coefficients
 */
template <uint32_t NO_INPUT_COEFFICIENTS, uint32_t NO_OUTPUT_COEFFICIENTS>
class IirFilter
{
public:
  /**
   * @brief Constructor accepting coefficients of the difference equation.
   * The coefficients are expected in following order:
   * 
   * $\f b_0, b_1, \ldots , b_M, -a_1, -a_2, \ldots, a_N $\f 
   * 
   * i.e. at the beginning $\f M+1 $\f \emph{b} coefficients followed by $\f N $\f
   * \emph{a} coefficients with negative sign.
   */
  template <typename... Args>
  constexpr IirFilter(const Args &... args) :
    input_buffer{},
    output_buffer{},
    input_index{0},
    output_index{0},
    coefficients{args...}
  {
  }

  /**
   * @brief Method passes the filtered value into the filter
   * @param input filtered input
   */
  void setInput(float input)
  {
    input_buffer[input_index++] = input;
    if (input_index == NO_INPUT_COEFFICIENTS) {
      input_index = 0;
    }
  }

  /**
   * @brief Method calculates the filter.
   */
  void calculate()
  {
    // filter implemented in the direct form
    convolveInputs();
    convolveOutputs();
    calculateOutput();
    storeOutput();
  }

  /**
   * @brief Method returns output of the filter.
   * @return filter output
   */
  float getOutput() const
  {
    return output;
  }

private:
  static const uint32_t kFirstInputCoefficientIndex = 0;
  static const uint32_t kFirstOutputCoefficientIndex =
      kFirstInputCoefficientIndex + NO_INPUT_COEFFICIENTS;
  /**< Circular buffer for the input samples x(k) */
  float input_buffer[NO_INPUT_COEFFICIENTS];
  /**< Circular buffer for the output samples y(k) */
  float output_buffer[NO_OUTPUT_COEFFICIENTS];
  /**< Coefficients of the difference equation */
  const float coefficients[NO_INPUT_COEFFICIENTS + NO_OUTPUT_COEFFICIENTS];
  /**< Position of the current oldest sample of the input sequence in the
   * input circular buffer */
  uint32_t input_index;
  /**< Position of the current oldest sample of the output sequence in the
   * output circular buffer */
  uint32_t output_index;
  /**< Convolution of the last $\fM+1$\f input samples */
  float input_convolution;
  /**< Convolution of the last $\fN$\f input samples */
  float output_convolution;
  /**< Current output $\fy(k)$\f */
  float output;

  /**
   * @brief Method calculates the convolution of the last $\fM+1$\f
   * input samples i.e.
   * 
   * $\f
   * b_0\cdot x(k) + b_1\cdot x(k-1) + \ldots + b_M\cdot x(k-M)
   * $\f
   * 
   */
  void convolveInputs()
  {
    input_convolution =
        convolve(input_buffer, NO_INPUT_COEFFICIENTS, input_index, coefficients,
                 kFirstInputCoefficientIndex, NO_INPUT_COEFFICIENTS);
  }

  /**
   * @brief Method calculates the convolution of the last $\fN$\f
   * input samples i.e.
   * 
   * $\f
   * -a_1\cdot y(k-1) - a_2\cdot y(k-2) - \ldots - a_N\cdot y(k-N)
   * $\f
   * 
   */
  void convolveOutputs()
  {
    output_convolution = convolve(
        output_buffer, NO_OUTPUT_COEFFICIENTS, output_index, coefficients,
        kFirstOutputCoefficientIndex, NO_OUTPUT_COEFFICIENTS);
  }

  /**
   * @brief Method calculates current sample of the output sequence
   * $\f y(k) $\f
   */
  void calculateOutput()
  {
    output = input_convolution + output_convolution;
  }

  /**
   * @brief Method inserts current sample of the output sequence into
   * the output buffer 
   */
  void storeOutput()
  {
    output_buffer[output_index++] = output;
    if (output_index == NO_OUTPUT_COEFFICIENTS) {
      output_index = 0;
    }
  }

  /**
   * @brief Method calculates convolution of the sequence stored
   * in the circular buffer with given impulse response.
   * @param circular_buffer circular buffer where the input sequence is stored
   * @param circular_buffer_length length of the circular buffer
   * @param circular_buffer_oldest_sample_index index of the current oldest 
   * sample of the input sequence in the circular buffer
   * @param impulse_response array where the impulse response samples are stored
   * @param impulse_response_first_sample_index index of the first sample of the
   * impulse response
   * @param impulse_response_length number of samples per impulse response
   * @return convolution of the sequence stored in the circular buffer with given 
   * impulse response
   */
  float convolve(const float circular_buffer[],
                 const uint32_t circular_buffer_length,
                 const uint32_t circular_buffer_oldest_sample_index,
                 const float impulse_response[],
                 const uint32_t impulse_response_first_sample_index,
                 const uint32_t impulse_response_length)
  {
    float convolution = 0;
    // position of the last inserted i.e. newest sample of the input sequence
    int32_t j = circular_buffer_oldest_sample_index - 1;
    if (j < 0) {
      // last inserted sample is at the end of the buffer
      j = circular_buffer_length - 1;
    }
    // iterate over the last "impulse_response_length" samples of the input
    // sequence in direction from the "newest" to the "oldest" sample
    for (uint32_t i = impulse_response_first_sample_index;
         i < impulse_response_first_sample_index + impulse_response_length;
         i++) {
      convolution += impulse_response[i] * circular_buffer[j];
      if (--j < 0) {
        j = circular_buffer_length - 1;
      }
    }
    return convolution;
  }
};

#endif /* IIRFILTER_H */

Below is a code example documenting usage of the IirFilter class (namely calculation of a step response of a filter with following transfer function $$H(z) = \frac{0.0008663387 + 0.001732678\cdot z^{-1} + 0.0008663387\cdot z^{-2}}{1 - 1.919129\cdot z^{-1} + 0.9225943\cdot z^{-2}}$$)

int main(int argc, char** argv) {

    const uint32_t kInputLength = 64;        
    float x[kInputLength] = {1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f, 
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f, 
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f, 
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f, 
                             1.0f, 1.0f, 1.0f, 1.0f,
                             1.0f, 1.0f, 1.0f, 1.0f};

    IirFilter<3, 2> iir(0.0008663387f, 0.001732678f, 0.0008663387f, 
                        1.919129f, -0.9225943f);

    for (uint32_t i = 0; i < kInputLength; i++) {
        iir.setInput(x[i]);
        iir.calculate();
        std::cout << iir.getOutput() << std::endl;
    }
    
    return 0;
}

Answer 1

Doxygen errors

It's great that you are using Doxygen to document the code, and also that you include the formulas that your filter implements. However, you did not use the right opening and closing tags for the formulas. You have to use \f$...\f$ instead of $\f...$\f. Furthermore, you cannot use \$\LaTeX\$ commands in the normal text in Doxygen, use Doxygen's own commands instead. For example, you cannot write \emph{a} coefficients, you have to write \em a coefficients if you really just wanted emphasis, but I think it's better to write \f$a\f$-coefficients.

You are also missing @params for the constructor of IirFilter.

Make sure you validate that your Doxygen documentation is correct by creating a Doxygen config file if you haven't already, and set all WARN_* parameters, including WARN_AS_ERROR, to YES.

Let functions take parameters and return values

There is no reason to make calculate() a function that takes no parameters and returns no value, and then have separate functions setInput() and getOutput() to set the input to the calculation and get the result of the calculation. Instead, make it look like this:

float calculate(float input) {
    setInput(input);
    ...
    return getOutput();
}

Now also apply this principle to the private member functions. And once you've done this with everything, you'll note while before you needed to store all intermediate results in member variables to pass them from function to function, now you can avoid that in several places. Consider rewriting calculate() like so:

float calculate(float input) {
    setInput(input);
    float output = convolveInputs() + convolveOutputs();
    storeOutput(output);
    return output;
}

Make more member functions const

You made getOutput() const, which is great, but there are more functions that don't modify any of the member variables, and thus should be marked as being const, like convolve(). And if you modified the functions to take and return values as mentioned above, then you should be able to make convoleInputs() and convolveOutputs() const as well.

Consider implementing a circular buffer class

You don't have actual circular buffers in your code, you only have the array coefficients[]. It is convolve() that not only does the convolution, but also has to implement the circular buffer semantics. Ideally, you would have a circular buffer class that stores the input and output values, and that allows you to push new values into them, and that acts like a regular STL container, including providing iterators. Then convolve() would be much simpler. In fact, you wouldn't even need to implement convolve() yourself anymore, you could then just use std::inner_product(). Consider being able to write:

#include <numeric>

template <...>
class IirFilter {
public:
    ...
    float calculate(float input) {
         inputs.push(input);

         float input_convolution = std::inner_product(
             inputs.begin(), inputs.end(),
             coeffients.begin(),
             0.f);

         float input_convolution = std::inner_product(
             outputs.begin(), outputs.end(),
             coeffients.begin() + NO_INPUT_COEFFICIENTS,
             0.f);

         float output = input_convolution + output_convolution;
         outputs.push(output);
         return output;
    }
    ...
private:
    circular_buffer<float, NO_INPUT_COEFFICIENTS> inputs;
    circular_buffer<float, NO_OUTPUT_COEFFICIENTS> outputs;
    ...
};