I've implemented something similar to calculate the rate/concentration of particles in an air stream. They come at random events (probably poisson distributed) and I want to know the average rate in particles per time. My approach is as follows (taken from the docstring of my code):
Event timestamps are placed in a buffer of a fixed maximum size.
Whenever a concentration estimate is requested, all timestamps up to
a set maximum age are filtered out. If the remaining number of
recent timestamps is below some threshold, the reported
concentration is zero. Otherwise, the concentration is calculated
as the number of remaining events divided by the time difference
between the oldest timestamp in the buffer and now.
I'll attach the Python implementation below for reference. It is a part of a much larger project, so I had to slightly modify it to get rid of references to external code:
particle_concentration_estimator_params.py
#!/usr/bin/env python3
"""This module implements a definition of parameters for the particle
concentration estimator.
"""
__author__ = "bup"
__email__ = "bup@swisens.ch"
__copyright__ = "Copyright 2021, Swisens AG"
__license__ = "GNU General Public License Version 3 (GPLv3)"
__all__ = ['ParticleConcentrationEstimatorParams']
import dataclasses
@dataclasses.dataclass
class ParticleConcentrationEstimatorParams:
"""This provides storage for the parameters used for the particle
concentration estimator.
"""
timestamp_buffer_max_size: int
"""Maximum size of the buffer that is used to keep track of event
timestamps. The size of this buffer mainly affects the filtering
of the reported data.
Unit: - (count)
"""
timestamp_buffer_max_age: float
"""Maximum age of events in the timestamp buffer which are
considered for the concentration calculation. This value is a
tradeoff between a smooth filtered value and the dynamic response
to a changed concentration.
Unit: s
"""
min_number_of_timestamps: int
"""Minimum number of timestamps to use for the concentration
estimation. If less timestamps are available, the concentration is
reported as zero.
Unit: - (count)
"""
particle_concentration_estimator.py
#!/usr/bin/env python3
"""This module implements the particle concentration estimation.
"""
__author__ = "bup"
__email__ = "bup@swisens.ch"
__copyright__ = "Copyright 2021, Swisens AG"
__license__ = "GNU General Public License Version 3 (GPLv3)"
__all__ = ['ParticleConcentrationEstimator']
import logging
import time
from typing import Optional
import numpy as np
from .particle_concentration_estimator_params import ParticleConcentrationEstimatorParams
logger = logging.getLogger(__name__)
class ParticleConcentrationEstimator:
"""An object of this class implements the Poleno particle
concentration estimator. Particle concentration is basically just
a number of particles per time unit. But since the particle events
arrive irregularly, there are various ways to filter the result, to
avoid too much noise especially when the concentration is low. This
class implements the following approach:
Event timestamps are placed in a buffer of a fixed maximum size.
Whenever a concentration estimate is requested, all timestamps up to
a set maximum age are filtered out. If the remaining number of
recent timestamps is below some threshold, the reported
concentration is zero. Otherwise, the concentration is calculated
as the number of remaining events divided by the time difference
between the oldest timestamp in the buffer and now.
"""
def __init__(self, params: ParticleConcentrationEstimatorParams):
"""Initializes the object with no events.
Args:
est_params: Initialized PolenoParams object which includes
information describing the estimator's behaviour.
"""
self.params = params
"""Stored params for the object."""
n_rows = self.params.timestamp_buffer_max_size
self._rb = np.full((n_rows, 2), -1e12)
self._rb_wp = 0
self._num_timestamps = 0
self._concentration_value = 0.0
self._concentration_value_no_mult = 0.0
def tick(self, now: float) -> float:
"""Recalculates the current concentration value.
Args:
now: Current timestamp to use to filter out old entries
in the buffer.
Returns:
The updated concentration value, which is also returned
using the concentration attribute.
"""
min_ts = now - self.params.timestamp_buffer_max_age
min_num = self.params.min_number_of_timestamps
used_rows = self._rb[:, 0] >= min_ts
filt_ts = self._rb[used_rows]
num_ts = round(np.sum(filt_ts[:, 1]))
self._num_timestamps = num_ts
num_ts_no_mult = round(np.sum(used_rows))
if num_ts < min_num:
self._concentration_value = 0.0
self._concentration_value_no_mult = 0.0
else:
t_diff = now - np.min(filt_ts[:, 0])
if t_diff >= 1e-3:
# Do not change the reported value if all events in the
# buffer have the same timestamp.
self._concentration_value = num_ts / t_diff
self._concentration_value_no_mult = num_ts_no_mult / t_diff
return self._concentration_value
def got_timestamp(self,
ts: Optional[float] = None,
multiplier: float = 1.0) -> None:
"""Passes in the most recent event timestamp. Timestamps need
not be ordered.
Calling this method does not immediately update the
concentration value, this is deferred to the tick() method.
Args:
ts: Event timestamp to use. If None, the current time is
used.
multiplier: Optional multiplier, which makes ts to be
counted that many times.
"""
if ts is None:
ts = time.time()
self._rb[self._rb_wp] = (ts, float(multiplier))
self._rb_wp = (self._rb_wp + 1) % self._rb.shape[0]
self._num_timestamps += round(multiplier)
@property
def concentration(self) -> float:
"""The calculated concentration value.
Unit: 1/s
"""
return self._concentration_value
@property
def concentration_no_mult(self) -> float:
"""The calculated concentration value without taking into
account the timestamp multipliers, i.e. as if all timestamps
were given with the default multiplier of 1.
Unit: 1/s
"""
return self._concentration_value_no_mult
@property
def num_timestamps(self) -> int:
"""Returns the number of timestamps which currently are in
the internal buffer.
"""
return self._num_timestamps
