Find smoothed first derivative from signal with noisy slope 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:


Take derivative as dx(i) = (x(i)-x(i-99))/100
Smooth with sliding mean window = 100, then take derivative as dx(i) = (x(i)-x(i-99))/100
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: 


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).
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!).
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?
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 &lt;- 1000
k1 &lt;- 0.006998
k2 &lt;- 0.0188728493
k_true &lt;- c(k1*rep(1,T/2), k2*rep(1,T/2))
c_true &lt;- c(2.90238432*rep(1,T/2), 2.90238432*rep(1,T/2) + (k1-k2)*(T/2+1))
t &lt;- seq(1,T)

y_true &lt;- k_true*t + c_true

y_noisy &lt;- y_true + rnorm(T)*0.1

k_hat &lt;- rep(0,T)


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

  xy &lt;- sum((window_idx - x_bar)*(y_noisy[window_idx] - y_bar))
  xx &lt;- sum((window_idx - x_bar)*(window_idx - x_bar))

  k_hat[t_idx] &lt;- 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')