2
$\begingroup$

I generate random spike data that represent the ouput from a rotary encoder. Here the output of my algorithm: enter image description here The first plot is position, then velocity and then acceleration. Even though I downsample my position vector and pass it to a average filter of 20 points, the derivative are very sensitive as you can see. I check the following question: have position, want to calculate velocity and acceleration. Few option seem to be available , but I don't get the answer , the Savitzky-Golay Filters is just a smoothing function, I don't get how he get the velocity, also what are the orther alternatives.

$\endgroup$
6
  • 4
    $\begingroup$ there appears to be no difference between velocity and acceleration. are you sure the third plot is correct? that said, differentiating noisy data get even noisier, usually because noise on data is more high-frequency in content (low-frequency noise is often considered a "trend", not noise). so the noisy thing is normal. but the accleration plot looks just like the velocity plot and that can't be right. $\endgroup$ Oct 5, 2015 at 22:10
  • $\begingroup$ The Savitzky-Golay filters can give you a smoothed derivative as well, which is what's used in that answer. Look at the parameters for the sg_filter() method in the code provided (i.e. "k = Which derivative"). $\endgroup$
    – datageist
    Oct 5, 2015 at 23:24
  • $\begingroup$ What is "my algorithm" ? i.e. how are you converting from rotary encoder pulses to position, velocity, and acceleration estimates. If you use a Kalman filter to do this conversion, you may get better results --- or be better able to tweak the smoothing parameters. $\endgroup$
    – Peter K.
    Oct 5, 2015 at 23:50
  • $\begingroup$ @robertbristow-johnson well my velocity is garbage so doing a diff function from matlab on that will most likely end up a garbage result $\endgroup$
    – MathieuL
    Oct 6, 2015 at 0:22
  • $\begingroup$ @PeterK. I got 2 channel on my encoder so knowing which pulse happen first on a channel I know the rotation of the shaft. Therefore I add the distance that has been travel by the shaft to the previous position. $\endgroup$
    – MathieuL
    Oct 6, 2015 at 0:25

1 Answer 1

3
$\begingroup$

There are several ways to calculate the derivative of a stochastic signal. I will show you two methods.

First of all, adjust the output of a second order transfer function to your data: $${\omega_n^2}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ Once you are satisfied with the fit, put an additional $s$ in the numerator and you will get the derivative of your signal: $${\omega_n^2\times s}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ You can use a bilinear transformation and convert this solution into a discrete one, to implement in a software easily, if you like.

Another solution is to use the spectrum and calculate the fft of your signal. Then filter it vanishing the small values in the higher frequencies. Now, take the result and multiply it by $j\omega$. This is taken from the Fourier Transform Properties Table to calculate the derivative of a signal. Finally, redo the signal with ifft.

I added here a code to exemplify these aforementioned algorithms. The output picture is added as well at the end of the post. I hope it helps.

% \brief: this code exemplifies two techniques suitable for
%         calculating the first-order derivative of an
%         stochastic signal.
%
% \algorithm 1: LTI transfer function integrated with Tustin.
% \algorithm 2: FFT derivative property.
%
% \author:     Luciano Augusto Kruk
% \web:        www.kruk.eng.br

% some constants:
fs = 1000; % [Hz] sample rate
T  = 1/fs;
t  = 0:T:0.5;
nT = length(t);
% signal + colored noise:
[B,A] = butter(10,250*T*2);
s     = sin(2*pi*20*t)+filter(B,A,0.1*randn(1,nT));

figure(1); clf;
subplot(3,1,1);
plot(t,s);
title('signal s(t)');
grid on;

% first technique: second order transfer function
qsi  = 0.7;
wn   = 2*pi*60;
wnT2 = (wn*T)^2;
b    = 2*T*wn*wn;
a1   = 4 + (4*qsi*wn*T) + wnT2;
a2   = (2*wnT2)-8;
a3   = wnT2 - (4*qsi*wn*T) + 4;
dsdt = zeros(size(s));
for i = 3:nT
    dsdt(i) = (1/a1) * ...
        ((-a2*dsdt(i-1)) - (a3*dsdt(i-2)) + (b * (s(i) - s(i-2))));
end
subplot(3,1,2);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 1');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')

% second technique: fft/ifft properties
s_fft = fft(s);
idx   = 1:(nT/2);
f_fft = (idx-1) * (fs/2) * (1/length(idx));
idx   = round((nT/2) + ((-nT/2.5):(nT/2.5)));
s_fft(idx) = 0;
w     = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
dsdt  = ifft(s_fft .* sqrt(-1) .* ifftshift(w));

subplot(3,1,3);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 2');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')

derivative estimate of stochastic signal

$\endgroup$
3
  • 1
    $\begingroup$ The fft approch seem to be very nice. I will look into that $\endgroup$
    – MathieuL
    Oct 6, 2015 at 0:23
  • 1
    $\begingroup$ Can you please supply your position samples? I would like to test an idea. $\endgroup$ Oct 6, 2015 at 16:49
  • $\begingroup$ I misread the question. I understood to derive position from acceleration... $\endgroup$
    – Brethlosze
    Dec 30, 2016 at 23:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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