2
$\begingroup$

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.

slope noise

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?

$\endgroup$
  • $\begingroup$ Have you tried point 1. with different steps, i.e. $(x(i)-x(i-T))/(T+1)$ and then average over all $T$? $\endgroup$ – firion Jul 12 '17 at 14:13
  • 1
    $\begingroup$ Welcome to SE.DSP! Interesting question. How high a computational cost is too high? $\endgroup$ – Peter K. Jul 12 '17 at 14:29
3
$\begingroup$

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.

enter image description here

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.

enter image description here


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')
$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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