Find smoothed first derivative from signal with noisy slope

Question

How can I get filtered first derivative from a noisy signal that has slowly changing slope in form of y=kx+b? k can slowly change in time, and I want to estimate its value.

I have tried 3 different approaches:

1. Take derivative as dx(i) = (x(i)-x(i-99))/100
2. Smooth with sliding mean window = 100, then take derivative as dx(i) = (x(i)-x(i-99))/100
3. Simple IIF filters (e.g. y(i) = 0.99*y(i-1) + 0.01*x(i), then take derivative as dx(i) = y(i)-y(i-1) and again filter with similar IIR, e.g. dy(i) = 0.95*dx(i-1) + 0.05*dx(i)

Problems:

1. Least-squares, regression and FIR filters (except rectangular window) have high computational cost, since I have to translate it to micro-controller with no DSP. That is why I can use only rectangular windows and IIR filters (they have low order).
2. If I find first derivative first, then smooth, it will be very noisy. So, I should smooth the original signal first, then find derivative from smoothed signal (and perhaps smooth the derivative again!).
3. I should play with filter parameters manually and it is hard to understand the frequency response of the whole system.

Question:

Maybe there is a single special (optimal?) IIR filter for this specific problem - finding smoothed first derivative from signal with noisy slope?

Answer 1

I think least squares is going to be the best approach, and that's not going to be that computationally expensive (I think! Please correct me if I'm wrong).

The gradient can be estimated from a sliding window of your data using: $$ \hat{k} = \frac{\sum (x_n - \bar{x})(y_n - \bar{y})}{\sum (x_n - \bar{x})^2 } $$ where the sum over $n$ is taken over the window.

If I implement an example in R (code below) then it seems to do the right thing. The picture shows the noiseless data (red) and noisy data (black) with a change in $k$ at half-way. The window length I've chosen is 20 points, and the output shifts at about the right point.

This estimate may still be too noisy for you; in that case, increase the window length. This picture shows what happens when window_length is increased to 100.

---

R Code Only Below

#42364

T <- 1000
k1 <- 0.006998
k2 <- 0.0188728493
k_true <- c(k1*rep(1,T/2), k2*rep(1,T/2))
c_true <- c(2.90238432*rep(1,T/2), 2.90238432*rep(1,T/2) + (k1-k2)*(T/2+1))
t <- seq(1,T)

y_true <- k_true*t + c_true

y_noisy <- y_true + rnorm(T)*0.1

k_hat <- rep(0,T)


window_length <- 20
for (t_idx in t)
{
  oldest_idx <- max(t_idx-window_length, 1)
  window_idx <- seq(oldest_idx, t_idx)
  x_bar <- sum(window_idx)/length(window_idx)
  y_bar <- sum(y_noisy[window_idx])/length(window_idx)

  xy <- sum((window_idx - x_bar)*(y_noisy[window_idx] - y_bar))
  xx <- sum((window_idx - x_bar)*(window_idx - x_bar))

  k_hat[t_idx] <- xy /xx
}

par(mfrow=c(2,1))
plot(t, y_noisy, type='l')
lines(t,y_true, col='red') #, lwd=3

plot(t,k_hat)
lines(t,c(k1*rep(1,T/2), k2*rep(1,T/2)), col='red')

Answer 2

You might try an $\alpha$ $\beta$ filter.

You have two "states" position which is the filtered x value and velocity which is a filtered derivative. Very simple 3 recursions , and a difference.

https://en.m.wikipedia.org/wiki/Alpha_beta_filter

The $\alpha$ and $\beta$ parameters require some tuning, and the article talks a bit about that.

The filter can also be upgraded to include another state or an adaptive gain rule if needed.

Answer 3

If anyone here is curious in a computationally efficient solution (say you are working on a small microcontroller like Arduino), you can run two alpha filters (Exponential Smoothing) with slightly different 'average ages' (similar to time constant).

'Average Age' Definition: https://en.wikipedia.org/w/index.php?title=Exponential_smoothing&oldid=1162913370#Average_Age

Here is an implementation of mine:

    % Required Parameters
    % dt = time since last reading
    % tau1 = desired average age of filter 1
    % tau2 = desired average age of filter 2

    % calculate alphas from desired average age: 
    alpha1 = 1/(tau1/dt + 1); % alpha for filter 1
    alpha2 = 1/(tau2/dt + 1); % alpha for filter 2

    % calculate filtered values
    value_filter1 = value_filter1*(1-alpha1) + new_value*alpha1;
    value_filter2 = value_filter2*(1-alpha2) + new_value*alpha2;

    % calculate derivative
    derivative = (error_filter1 - error_filter2) / (tau1 - tau2);

Answer 4

Since the time of question, I found this regression filter to be useful: http://www.claysturner.com/dsp/FIR_Regression.pdf

Filter coefficients for window size N:

$$\beta_i =\frac{12 \cdot i - 6(N-1)}{N(N^2-1)}$$

Code in Julia:

function slope_filter_coef(n::Int)
    s1 = n * (n+1) ÷ 2 # sum(1:N)
    s2 = n * (n+1) * (2n+1) ÷ 6 # sum((1:N).^2)
    b = map(n:-1:1) do i
        n*i - s1
    end
    gain = (n * s2 - s1^2)
    delay = (n-1) / 2
    return b, 1, gain, delay
end

using DSP, Plots
x = randn(500)
fs = 1000; fcut = 50; nr = 3
b, a, gain, delay = slope_filter_coef(fs, fcut, nr)
b = b ./ gain
y = DSP.filt(b, a, x)
y[1:end-delay] = y[1+delay:end]
y[end-delay+1:end] .= NaN
plot([x, y])

Also because coefs are linear, it can be optimized with recurrent form: $$ dy = y_{t} - y_{t-1} = \sum_{i=0}^{N-1}\beta_i \cdot x_{t-i} - \sum_{i=0}^{N-1}\beta_i \cdot x_{t-1-i} = \beta_0 (x_t + x_{t-N}) + d\beta \sum_{i=1}^{N-1} x_{t-i} $$

So we can update this filter O(1), independent from window length (storing buffer for N last points, accumulator for xs and previous y):

$$ y_t = y_{t-1} + \beta_0 (x_t + x_{t-N}) + d\beta \sum_{i=1}^{N-1} x_{t-i} $$ where $$\beta_i - \beta_{i-1} = d\beta = \frac{-12}{N(N^2-1)} = const $$ $$\beta_0 = - \beta_{N-1} = \frac{6}{N(N+1)} $$