# Simple low pass digital filter implementation for non-uniform sampling?

I have a system sampling an analog input and I want to filter that value using a simple, first order filter. EG something of the form:

Yn = A*Xn + B*Yn-1


However I have no control over the underlying process and the sampling period is non-uniform, and varies (randomlyish) somewhere between 25ms and 29ms. In addition I want to use a 30 minute time constant, and can't keep a running buffer of 72,000 sample values.

Given that a typical filter implementation requires a uniform sampling period, is there any simple algorithm that takes into account the time difference between Yn-1 and Xn? (which I guess really means a time varying version of my "constants" A and B)

Note that I'm inclined to simply take the average sample time over an extended period and using that to calculate my constants, however there is also the possibility of the final system not running at the same average rate as my test system, so I want to come up with a solution that is somewhat adaptable.

1. The filtered signal is being used as a part of a vibration monitoring system in some industrial equipment. It will be used for more of as a qualitative than quantitative purpose.
2. The customer requirements are 30 minute rolling average for the data.
3. I have zero control over the sample rate.
4. And the sample rate varies every second
5. And the sample rate defines the speed at which my code runs
6. The final system should run at about the same speed as my test system, however it could be somewhat slower, maybe down to a 35ms rate
• 1) What will the filtered signal be used for ? Non uniform sampling of the filtered output may cause problems down the line. 2) "In addition I want to use a 30 minute time constant," Can you elaborate more on this ? If your system time constant is so slow, do you really need to process the samples with an accuracy of 4 ms ? 3) Will the average sample time of the final system be fixed / known during operation of the final system ?
– AJN
Commented Apr 17, 2021 at 16:24
• Design a first order continuous time filter. Discretise it with sampling period Ts, keeping Ts as a variable and not a number. A and B will then be a function of Ts. Recalculate A and B for each sample if you know the Ts for that sample.
– AJN
Commented Apr 17, 2021 at 16:28
• S/N ratio improves by $\sqrt{n}$ or 10 log n in dB. Your desired to use a 1st order filter in 1/20 Hz Nyquist BW then then average it over 1800 samples then make it adaptable to the next average is contradictatory for S/N ratio . Pls specify your Noise spectrum and wandering average of final system in measurable terms 1st. Signal spectrum, then noise spectrum. It is better to oversample then use a 2nd order filter than a 6th order Nyquist filter and a 1st order digital filter is NG Commented Apr 17, 2021 at 16:28
• Why don't you resample to be uniform and at the same average rate? This can easily be done in realtime and for a few thousand ALU operations per second given your low sampling rate. Commented Apr 17, 2021 at 16:55
• @user1850479 I think I was thinking of oversampling as a physical rather than a computational operation. It looks like Andy's answer has set me straight on that Commented Apr 17, 2021 at 17:22

You say you have no control over the input sampling rate. That is understood but you are not considering the fact that you can resample once inside your digital system AND, I'm presuming you do have control over that else why would you be asking this question. So, resample your digital data like this: -

Now "real sample" 1 took 25 ms and "real sample" 2 took 29 ms but who cares if you internally resample that data? Yes, you have to process at a much higher rate but it's not like you are having to do major recalculations each time because $$\A*X_n\$$ is static for 25 to 29 samples.

• except that my code runs at the same speed as the sample rate. Or another way the sample rate is same speed as my code execution rate. I can't process data faster than the sample rate. Commented Apr 17, 2021 at 17:13
• Hmm .. however I could pretend that I run the code at a faster rate. EG design for a 1 ms rate and run X iterations at each code execution point, where X is the time between samples Commented Apr 17, 2021 at 17:15
• But you can execute hundreds (if not thousands) of lines of code at a much higher rate than the incoming flaky sample rate. Commented Apr 17, 2021 at 17:15
• I can only run code at the sample rate. They are in lock step Commented Apr 17, 2021 at 17:17
• Now I know. It's been a really long time since I've had to work in this area so I am not familiar with the terminology. But it looks like that is the solution to my problem Commented Apr 17, 2021 at 17:24

Consider the first order continuous time filter $$\\frac{dy(t)}{dt} = -k y(t) + k x(t)\$$

One of the many choices to discretise the above filter is

$$\frac{y_{(t)}-y_{(t-T_s)}}{Ts} = -k y_{(t-T_s)} + k x_{(t)}$$ $$y_{(t)} = k\cdot T_s x_{(t)} + (1 - k\cdot T_s)y_{(t- T_s)}$$

A and B to be used for each sample can be found from above expression.

The above is only an ad-hoc method. There may be proper algorithms for variable sampling rate filtering available in lierature with well defined mathematical properties.

A suggestion from From DSP.SE

Remember that filters are dynamical systems and solve the corresponding ordinary differential equation (ODE) In the second method, you basically use the retrieved time step as the time step in the ODE solver. Assuming you have a second order system:

Here comes the math. Let me assume, your low pass filter has a time constant $$\\tau\$$. It takes amplitude inputs values $$\a_n\$$ at points in time $$\t_n\$$ and it shall produce ouptut values $$\b_n\$$ for the same points in time. The starting condition is $$\a_0 = 0\$$ at $$\t_0 = 0\$$. So I initialize

$$\b_0 = 0\$$

All further output values can be computed as a weighted average of the previous output and the current input.

$$\b_n = w_n a_n + (1-w_n) b_{n-1}\$$

The weights $$\w_n\$$ can be computed from the well known low pass repsonse in time to the unit step function.

$$\w_n = 1 - exp(\frac{-\Delta t_n}{\tau} )\$$

where

$$\\Delta t_n = t_n - t_{n-1}\$$

So that is basically it. Here comes an example with $$\\tau = 32\$$. I have used a far higher non-uniformity for the sampling period as in the given question to make the effects more visible.

When working on a restricted processor such as a tiny microcontroller, calculating the weights nonlinearily as $$\w_n = 1 - exp(-\Delta t_n / \tau)\$$ may be costly. In those cases a linear approximation may serve well.

$$\w_n \approx min(1, \frac{\Delta t_n}{\tau})\$$

I tried it with the depicted example and the results cannot be distinguished visually.

There are ways to do what you ask for exactly... but you don't really need to do that. A 30 minute running average will erase all of your hard work.

Start by converting to a regular 1-second sampling rate just by taking the average of all the sample values within each second.

Then apply regular-sampling methods on the result.

• What you are suggesting is basically what AJN suggests. However in your case you are suggesting converting the data into 1 second blocks vs oversampling in order to treat the data as 1ms samples. Programmatically both approaches will end up having a similar amount of instructions. Computationally your approach could be more efficient, but that isn't a concern in my system. Numerically your approach may require variables of a smaller bit width, but again I have sufficient overhead that I can absorb AJNs approach. And finally I've already implemented it and I'm not redoing it! LOL Commented Apr 20, 2021 at 14:20
• lol at "could" :). I'm sure your time is more valuable than the CPU you would save, but oversampling by 30 times in order to decimate by 60000 will get you disapproving glances at all the DSP parties. Commented Apr 20, 2021 at 15:16
• The task was implemented in a manner that works. However my time has a dollar cost associated with it, and with all the other tasks performed on that project the purchase order was totally drained. So arbitrarily changing one working methodology to another one would mean working without a PO - effectively working for free. That is something I only do when there is a really compelling reason, and this doesn't even come close. And I also don't go to those sorts of parties :P Commented Apr 20, 2021 at 15:47