Basic digital RC approximation filter in python (Micropython)

Question

I have written a simple RC Filter approximation in python that I have a feeling could be more concise. Are there any obvious improvements?

class EWMA:
    def __init__(self, coeff: list, initialValue: float):
        ##
        # ewma3 states for coefficient optimization
        ##
        self.last_ewma3 = [0.0, 0.0, 0.0]

        ##
        # ewma6 states for coefficient optimization
        ##
        self.last_ewma6 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

        ##
        # default coefficients to 1.0 so the order can be from 0 - 6
        # since cascade elements will pass input signal to output with a=1
        ##
        self.coeff = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

        for c in range(0, len(coeff)):
            if(c >= len(self.coeff)):
                print(f'EWMA Coefficients Length Mismatch! len(coeff) = {len(coeff)}, max is 6')
                break
            self.coeff[c] = coeff[c]
        ##
        # realtime filter states
        ##
        self.states = [0, 0, 0, 0, 0, 0, 0]
        self.states[0] = initialValue
        self.states[1] = initialValue
        self.states[2] = initialValue
        self.states[3] = initialValue
        self.states[4] = initialValue
        self.states[5] = initialValue
        self.states[6] = initialValue

    def preload(self, value: float):
        self.states[0] = value
        self.states[1] = value
        self.states[2] = value
        self.states[3] = value
        self.states[4] = value
        self.states[5] = value
        self.states[6] = value

    ##
    # @brief      calculate single EWMA element
    ##
    # @param      self   The object
    # @param      alpha  filter coefficient
    # @param      this   current input sample
    # @param      last   last output sample from this stage (feedback)
    ##
    # @return     EWMA result
    ##
    def ewma(self, alpha: float, this: float, last: float) -> float:
        return (float(alpha)*float(this)) + ((1.0-float(alpha))*float(last))

    ##
    # @brief      calculate 6th order cascade ewma
    ##
    # @param      self        The object
    # @param      inputValue  Raw input sample
    ##
    # @return     output of 6th cascade element
    ##
    def calculate(self, inputValue: float) -> float:
        result = 0.0
        self.states[0] = float(inputValue)
        self.states[1] = self.ewma(float(self.coeff[0]), self.states[0], self.states[1])
        self.states[2] = self.ewma(float(self.coeff[1]), self.states[1], self.states[2])
        self.states[3] = self.ewma(float(self.coeff[2]), self.states[2], self.states[3])
        self.states[4] = self.ewma(float(self.coeff[3]), self.states[3], self.states[4])
        self.states[5] = self.ewma(float(self.coeff[4]), self.states[4], self.states[5])
        self.states[6] = self.ewma(float(self.coeff[5]), self.states[5], self.states[6])
        return self.states[6]

    def get_last_output(self) -> float:
        return self.states[6]

    def model_ewma3_preload(self, v: float):
        self.last_ewma3[0] = v
        self.last_ewma3[1] = v
        self.last_ewma3[2] = v

    ##
    # @brief      ewma 3rd order for IIR Model Fitting via SciPy Optimize
    ##
    # @param      self  The object
    # @param      y0    The input value
    # @param      a     coeff a
    # @param      b     coeff b
    # @param      c     coeff c
    ##
    # @return     IIR output
    ##
    def model_ewma3(self, y0, a, b, c):
        y1 = self.ewma(a, y0, self.last_ewma3[0])
        y2 = self.ewma(b, y1, self.last_ewma3[1])
        y3 = self.ewma(c, y2, self.last_ewma3[2])
        self.last_ewma3[0] = y1
        self.last_ewma3[1] = y2
        self.last_ewma3[2] = y3
        return y3

    def model_ewma6_preload(self, v):
        self.last_ewma6[0] = v
        self.last_ewma6[1] = v
        self.last_ewma6[2] = v
        self.last_ewma6[3] = v
        self.last_ewma6[4] = v
        self.last_ewma6[5] = v

    ##
    # @brief      ewma 6th order for IIR Model Fitting via SciPy Optimize
    ##
    # @param      self  The object
    # @param      y0    The Input Value
    # @param      a     coeff a
    # @param      b     coeff b
    # @param      c     coeff c
    # @param      d     coeff d
    # @param      e     coeff e
    # @param      f     coeff f
    ##
    # @return     { description_of_the_return_value }
    ##
    def model_ewma6(self, y0, a, b, c, d, e, f):
        y1 = self.ewma(a, y0, self.last_ewma3[0])
        y2 = self.ewma(b, y1, self.last_ewma3[1])
        y3 = self.ewma(c, y2, self.last_ewma3[2])
        y4 = self.ewma(d, y3, self.last_ewma3[3])
        y5 = self.ewma(e, y4, self.last_ewma3[4])
        y6 = self.ewma(f, y5, self.last_ewma3[5])
        self.last_ewma6[0] = y1
        self.last_ewma6[1] = y2
        self.last_ewma6[2] = y3
        self.last_ewma6[3] = y4
        self.last_ewma6[4] = y5
        self.last_ewma6[5] = y6
        return y6

    def get_cutoff(self, Fs: float=1.0) -> float:
        x = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
        try:
            x[0] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[0], 2)/(2.0*(1.0 - self.coeff[0]))))
            print(f"Tap 1 {x[0]}")
        except:
            print("filter tap 1 not initialized")
        try:
            x[1] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[1], 2)/(2.0*(1.0 - self.coeff[1]))))
            print(f"Tap 2 {x[1]}")
        except:
            print("filter tap 2 not initialized")
        try:
            x[2] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[2], 2)/(2.0*(1.0 - self.coeff[2]))))
            print(f"Tap 3 {x[2]}")
        except:
            print("filter tap 3 not initialized")
        try:
            x[3] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[3], 2)/(2.0*(1.0 - self.coeff[3]))))
            print(f"Tap 4 {x[3]}")
        except:
            print("filter tap 4 not initialized")
        try:
            x[4] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[4], 2)/(2.0*(1.0 - self.coeff[4]))))
            print(f"Tap 5 {x[4]}")
        except:
            print("filter tap 5 not initialized")
        try:
            x[5] = (Fs/2*math.pi)*math.acos(1.0 - (math.pow(self.coeff[5], 2)/(2.0*(1.0 - self.coeff[5]))))
            print(f"Tap 6 {x[5]}")
        except:
            print("filter tap 6 not initialized")

        return x

    def apply_to_data(self, data: list) -> list:
        output = []
        for d in data:
            output.append(self.calculate(d))
        return output

it is worth it to note that the following class functions are only used when optimizing the filter coefficients for system identification purposes, so its not the end of the world if they are sloppy to me. The remaining functions will be used in real-time:

def model_ewma3_preload(self, v: float):
def model_ewma3(self, y0, a, b, c):
def model_ewma6_preload(self, v):
def model_ewma6(self, y0, a, b, c, d, e, f):
def get_cutoff(self, Fs: float = 1.0) -> float:

usage is akin to this:

# called on a timer at ~100Hz
sensor_reading = sensor.read()
system_approximation = filter.calculate(sensor_reading)

Welcome to Code Review! Your indentation was broken, probably when pasting the code here. The easiest way to avoid that is to paste the code directly from your editor, select all of it and then either click the {} button at the top or press Ctrl+k. — AlexV, Commented Sep 20, 2019 at 16:13
Great, Thank you for the corrections! I will remember next time. — Luke Gary, Commented Sep 20, 2019 at 16:19
For this and future questions, it'll be important to mention that this is not standard Python, but in fact MicroPython. — Reinderien, Commented Sep 20, 2019 at 16:36
Can you share some code that exercises this class in an application-typical way? — Reinderien, Commented Sep 20, 2019 at 17:19

AlexV · Accepted Answer · 2019-09-20 16:50:10Z

Style

Python comes with an "official" Style Guide for Python Code (often just called PEP8). Among others, it lists conventions regarding function documentation. In Python they're usually called docstrings (further detailed in PEP257) and written in """triple quotes""" after the function definition with def. For example:

def ewma(self, alpha: float, this: float, last: float) -> float:
    """Calculate single EWMA element

    Calculate EWMA with alpha being the filter coefficient, this the current
    input sample, and last the previous output value
    """

Of course this is only, if you don't have other conventions to follow. Seems like your code has some kind of doxygen-like syntax, but IIRC doxygen support for Python is not terribly well. If you're looking for a more structured approach that is better supported, numpydoc in conjunction with Sphinx might be an option to consider. The same example using numpydoc:

def ewma(self, alpha: float, this: float, last: float) -> float:
    """Calculate single EWMA element

    Parameters
    ----------
    alpha : float
        filter coefficient
    this : float
        current sample
    last : float
        the previous sample

    Returns
    -------
    ewma_result : float
        the result of the EMWA computation
    """

This combination is especially used in the "scientific Python stack" (numpy, scipy, ...).

The code itself

Use ~~numpy~~ (seems like you don't want to AND @Reinderien has beaten me ;-))! It would make your code a lot easier to read, and is also likely faster. But let's focus on what you have already written:

Also you're doing a lot of work repeatedly. Python is a little bit "dumb", i.e. it will happily compute whatever you write (very likely) without realizing that the same computation happened just a few lines ago. An example of this would be (Fs/2*math.pi) inget_cutoff`. Compute it once, put it into a variable and reuse it.

Python also has some tricks up its sleeves that make working with lists a little bit easier, e.g. where you have:

def apply_to_data(self, data: list) -> list:
    output = []
    for d in data:
        output.append(self.calculate(d))
    return output

You could instead write

def apply_to_data(self, data: list) -> list:
    return [self.calculate(d) for d in data]

Depending on the actual implementation, list comprehensions might also be faster than manually appending.

Lists can also be built by "multiplication" which allows you to write something like self.states = [value] * 7.

Using try: ... catch: ... without specifying an exception will likely give you some headache, because it will simply catch all exceptions, including that triggered by pressing Ctrl+C. So you will never know for sure whether the error was something that you expected or not.

In terms of performance, [func(d) for d in data] is less efficient than [*map(func, data)], if func is a defined function. However, the former may be a bit more readable for some people. Also note that if func(d) is a simple expression like d*d, using an explicit for-loop is normally faster than the map counterpart where the expression needs to be wrapped into a lambda expression. — GZ0, Commented Sep 20, 2019 at 18:20
@GZ0: In the end it boils down to "it depends" and the code should be carefully profiled to find out if that part is the actual bottleneck. — AlexV, Commented Sep 20, 2019 at 19:18

Reinderien · Accepted Answer · 2019-09-20 18:11:31Z

Copy a list into a slice

This code:

for c in range(0, len(coeff)):
    if(c >= len(self.coeff)):
        print(f'EWMA Coefficients Length Mismatch! len(coeff) = {len(coeff)}, max is 6')
        break
    self.coeff[c] = coeff[c]

has a couple of problems. It's fine to validate the length of coeff, but it shouldn't be done like this. Also, don't do an element-by-element copy. Instead:

N = 6
self.coeff = [1]*N

if len(coeff) > N:
   raise ValueError(f'EWMA Coefficients Length Mismatch! {len(coeff)} > {N}')
self.coeff[:len(coeff)] = coeff

This:


    self.states = [0, 0, 0, 0, 0, 0, 0]
    self.states[0] = initialValue
    self.states[1] = initialValue
    self.states[2] = initialValue
    self.states[3] = initialValue
    self.states[4] = initialValue
    self.states[5] = initialValue
    self.states[6] = initialValue

should be

self.states = [initialValue] * (N+1)

And so on.

Docstrings

Move your function documentation into """triple quotes""" at the first line inside of your function.

Never `try/except`

This breaks Ctrl+C quit, and is too broad to be useful. Narrow your caught exception type.

Don't repeat yourself

get_cutoff needs to be rewritten as a loop over N values.

Don't cast unnecessarily

This:

def ewma(self, alpha: float, this: float, last: float) -> float:
    return (float(alpha)*float(this)) + ((1.0-float(alpha))*float(last))

already assumes that the inputs are floats - so don't call float again. Drop all of your casts.

Probable bug

def model_ewma6(self, y0, a, b, c, d, e, f):
    y1 = self.ewma(a, y0, self.last_ewma3[0])

Seems that you're using the wrong array here. Also - why are you hard-coding for 3rd- or 6th-order filters? Can you not just accept N as a parameter?

General

Once you've cleaned up your usage of lists, you should really consider moving to ~~numpy~~ the array module. It'll execute more quickly.

Suggested

This comes with a lot of caveats. Since you don't have test usage, I haven't been able to test it, so I don't know whether it's valid. You're going to want to develop a numerical test suite to ensure that it's calculating the right thing. I also assumed that there's no need to hard-code for 3rd- or 6th-order filters, so just added an n. Finally: I don't have micropython, so this is written naively, assuming that standard Python usage is valid.

from array import array
from math import pi, acos
from typing import Sequence, Iterable


class EWMA:
    def __init__(self, coeff: Sequence[float], initial_value: float, n: int = None):
        nc = len(coeff)
        if n:
            self.n = n
            if nc > n:
                raise ValueError(f'len(coeff)={nc} > n={n}')
        else:
            self.n = nc  # default to the length of the coefficients

        # ewma states for coefficient optimization
        self.last_ewma = array('f', (0 for _ in range(self.n)))

        # default coefficients to 1.0 so the order can be from 0 - n
        # since cascade elements will pass input signal to output with a=1
        self.coeff = array('f', (1 for _ in range(self.n)))
        self.coeff[:nc] = coeff

        self.states = array('f', (0 for _ in range(1 + self.n)))

    def preload(self, value: float):
        self.states[:] = value

    @staticmethod
    def ewma(alpha: float, this: float, last: float) -> float:
        """
        calculate single EWMA element
        :param alpha: filter coefficient
        :param this: current input sample
        :param last: last output sample from this stage (feedback)
        :return: EWMA result
        """
        return alpha*this + (1 - alpha)*last

    def calculate(self, input_value: float) -> float:
        """
        calculate nth order cascade ewma
        :param input_value: Raw input sample
        :return: output of nth cascade element
        """
        self.states[0] = input_value
        for i, (c, s) in enumerate(zip(self.coeff, self.states[:-1])):
            self.states[i + 1] = self.ewma(c, s, self.states[i + 1])

        return self.get_last_output()

    def get_last_output(self) -> float:
        return self.states[-1]

    def model_ewma_preload(self, v: float):
        self.last_ewma[:] = v

    def model_ewma(self, y0: float, coeffs: Sequence[float]) -> float:
        """
        ewma nth order for IIR Model Fitting via SciPy Optimize
        :param y0: The input value
        :param coeffs: Sequence of coefficients
        :return: IIR output
        """
        prev = y0
        for i, (c, e) in enumerate(zip(coeffs, self.last_ewma)):
            new = self.ewma(c, prev, e)
            self.last_ewma[i] = new
            prev = new
        return prev

    def get_cutoff(self, fs: float = 1) -> array:
        return array(
            'f',
            (
                fs*pi/2 * acos(1 - c**2/2/(1 - c))
                for c in self.coeff
            )
        )

    def apply_to_data(self, data: Iterable[float]) -> array:
        return array('f', (self.calculate(d) for d in data))

Thank you! I will make these edits and give it a shot. I am intentionally trying to stay away from numpy as I would like to run this on micropython (which does not support numpy). — Luke Gary, Commented Sep 20, 2019 at 16:21
Fair enough. Does micropython support docs.python.org/3/library/array.html ? — Reinderien, Commented Sep 20, 2019 at 16:22
Looks like this is the case! docs.micropython.org/en/latest/library/array.html — Luke Gary, Commented Sep 20, 2019 at 16:23
OK, great. I think particularly in an embedded environment it'll be important to efficiently represent your data, so array will probably improve your memory usage, and may (?) improve your runtime performance. I encourage you to attempt an implementation using it, and post a new question with your updated code. — Reinderien, Commented Sep 20, 2019 at 16:33

Stack Exchange Network

Basic digital RC approximation filter in python (Micropython)

2 Answers 2

Style

The code itself

Copy a list into a slice

Docstrings

Never `try/except`

Don't repeat yourself

Don't cast unnecessarily

Probable bug

General

Suggested

Your Answer

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

Hot Network Questions

Basic digital RC approximation filter in python (Micropython)

2 Answers 2

Style

The code itself

Copy a list into a slice

Docstrings

Never try/except

Don't repeat yourself

Don't cast unnecessarily

Probable bug

General

Suggested

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

Never `try/except`

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