How I can produce the double differentiation of digital signal data without knowing the function?

Hello all,

I have 5120 displacement data points. I want to double differentiate the data points to get the acceleration.

Example of my data points:x = [0 0.002 0.0039 0.0059 0.0078 0.0098 0.0117]; time(s) y = [-0.5570 0.1310 0.4010 0.0722 -0.5210 -0.8090 -0.7230]; displacement (m/s)

How I can produce the double differentiation of digital signal data without knowing the function?

Please help me!

Please advice.

Thank you for your time.

Regards, SMY

Answer 1

[Star Strider]

The del2 function can do this, but you have to be careful. Differentiation significantly amplifies noise.

[SMY]

Thank you Start Strider! :)

Yup. The differentiation will results a lot of noise.

Remember my last problem regarding the integration? after I integrate my data from acceleration to displacement using trapezium methods, now I want to check either I can get the same results when I do backward (displacement-acceleration).

I try d2ydx2 = diff(y(:),2)./diff(x(:),2) but seems not working. I dont why.

Just now, I try your suggestion using del2, but the answer totally different from original (acceleration)

[Star Strider]

My pleasure.

I’m not certain what we did before. It may have involved filtering, so you need to compare it with the filtered signal, not the original.

Using the diff function will result in the output being (in this instance) 2 elements shorter than the original, and only takes the differences of the elements. The del2 and gradient functions are much more mathematically sophisticated, and will produce output that is the same size as the input.

I would use the gradient function twice for the reconstruction:

h = mean(diff(x));                          % Step Size
acc = gradient(gradient(pos, h), h);

giving acceleration for position ‘pos’ and step size ‘h’. That should be reasonably accurate.

[SMY]

Thank you Star Strider,

Yes, You are right. Used diff will make the output shorter.

I already tried your suggestion. It works!!

You really help me.

Thank you for your kindness and time.

Regards, SMY

[Star Strider]

As always, my pleasure!

[SMY]

The differentiation works, but seems now I need to filter the signal to remove the noise. Because I need to achieved closer RMS between original signal (acceleration) and the signal that I got from double differentiation. Do you have any ideas what are the best differentiation filter for gradient?

[Star Strider]

There are several options. I suspect the noise in the gradient is broadband noise, so if you have the Signal Processing Toolbox, my first approach would be to use the Savitzky-Golay filter function, sgolayfilt. You will have to experiment with it, but that is so for everything in signal processing.

Otherwise, first do a Fourier transform fft on the differentiated signal to see if the noise is band-limited. If so, you can use a frequency selective filter. There are several ways to design filters in MATLAB, such as designfilt and dfilt. My IIR filter design procedure is here: How to design a lowpass filter for ocean wave data in Matlab?

If you have low-frequency baseline variation that you want to eliminate, as well as high-frequency noise, use a bandpass filter.

[SMY]

Thank you Star.

Can I have your email for further discussion. I really needs help. I really3 needs help.

[Star Strider]

My pleasure.

It’s easier to work here, especially since you can upload your data (as a .mat file) easily and I can work with it. It’s also easy for me to post code here.

I only need a representative sample of your data, not all of it, if uploading all of it would be a problem. Other upload sites aren’t always available to me.

I will provide what help I can.

[SMY]

Thank you Star,

&#1I upload my coding (matlabwork.m) for your references, comments and help. I don't know how to improve it anymore

&#1It work well to get my velocity but not for acceleration. The target error RMS for the differentiation signal need to be <= 10% of original signal

&#1I really appreciated your time

&#1Regards

&#1SMY

[Star Strider]

I do not entirely understand what you are doing. I used your code for guidance on what the .mat files contained. I did not run it because it is difficult for me to follow it.

One problem that I see immediately is that your interpretation of your frequencies may not be correct. When I plotted the PSD, I found the predominant frequency at 32 Hz with two smaller peaks at 40 and 48 Hz and a second harmonic at 65 Hz. (The FFT I calculated and plotted shows a similar spectrum.) If you are centring the passband of your filter at 60 Hz, you are eliminating almost all of your signal energy. (See figure(3).)

My code:

dv = load('SMY y1x1.mat');                                  % Velocity
v = dv.y1x1;
da = load('SMY y2x2.mat');                                  % Acceleration
a = da.y2x2;
dx = load('SMY yx.mat');                                    % Position
x = dx.yx;
dt = load('SMY t.mat');                                     % Time
t = dt.t;

figure(1)
subplot(3,1,1)
plot(t,x)
grid
title('Position')
axis([0  1    ylim])

subplot(3,1,2)
plot(t,v)
grid
title('Velocity')
axis([0  1    ylim])

subplot(3,1,3)
plot(t,a)
grid
title('Acceleration')
xlabel('Time (sec)')
axis([0  1    ylim])

Ts = mean(diff(t));                                         % Sampling Time Interval
Fs = 1/Ts;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
L = length(t);

ft_x = fft(x)/L;                                            % FFT Of ‘x’
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:length (Fv);                                         % Index Vector

figure(2)
plot(Fv, abs(ft_x(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('Fourier Transform (Position)')

[pxx,fp] = pwelch(x, [], [], [], Fs);
[pks,frqs] = findpeaks(pxx,fp, 'MinPeakDist',5, 'MinPeakHeight',2.0E-5);

txtcell = regexp(sprintf('%.0f Hz\n', frqs), '\n', 'split');

figure(3)
plot(fp, pxx)
hold on
plot(frqs, pks, '^r')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Power')
title('PSD (Position)')
text(frqs, pks+5.0E-6, txtcell(1:end-1))

[SMY]

&#1Thank you Star&#12

&#1Maybe I do not understand about the concept of the passband filter. After I plotting your coding, did you meant I need to passband the filter at 65Hz or 32Hz? or instead of using bandpass, I should used highpass filter to keep the significant energy?

&#1

[Star Strider]

My pleasure.

I would use a bandpass filter. You need to adjust the frequencies to include those you want to keep. I would definitely include the energy at 32 Hz, so I would begin with a passband [30 60]/Fn.

Use the freqz function with your filter design to be certain it has the passband characteristics you want it to.

Filter design usually requires some experimentation, so you will probably have to adjust the filter until it produces the output you want.

[SMY]

&#12Thank you Star for the suggestion&#12

&#12But still work to get velocity but not for acceleration. I feeling bad seriously.

&#12If you dont mind, can you show me how to get the Sine Wave, Square Wave and triangular wave based on my data?

[Star Strider]

‘If you dont mind, can you show me how to get the Sine Wave, Square Wave and triangular wave based on my data?’

I don’t mind, but I don’t see how you can get that from your data. When I plotted it (in figure(1) in my latest code), I just saw what is essentially an amplitude-modulated sine wave. I did not see square or triangle waves at all.

[SMY]

&#12What in my mind now, I would like to try differentiated the signal manual. So I would like to find the sine wave equation and from there I can do the differentiation manual. &#1

&#1I need to do all possible way to solve this problem. Maybe from the manual differentiation, I can see what are the problem with my coding.&#1

[Star Strider]

The Fourier transform gives you the frequencies, amplitudes, and phases of the sine wave components of your signal. If you have the Symbolic Math Toolbox, the differentiation will be easier. I calculated it in my code for the displacement, but only plotted the magnitudes. You can get the phase angles with the angle function.

The Wikipedia article on Product-to-sum and_sum-to-product identities and related pages will be helpful.

I’m still not certain what you want, since you have the position, velocity, and acceleration.

[SMY]

&#1 The velocity and acceleration data that I have now are my target signal. I need to differentiated the displacement to get the velocity and acceleration (measured signal). The RMS for my both measured signal need to achieved 10% or less from RMS my target signal&#1

[Star Strider]

When I add this code and modify the subplot calls slightly:

Dx = gradient(x, Ts);
D2x = gradient(Dx, Ts);

RMS_x = sqrt(mean(x.^2));
% RMS_Dx = sqrt(mean((Dx' - v).^2));
% RMS_D2x = sqrt(mean((D2x' - a).^2));
RMS_Dx = sqrt(mean((Dx').^2));
RMS_D2x = sqrt(mean((D2x').^2));

figure(1)
subplot(3,1,1)
plot(t,x)
grid
title('Position')
axis([0  1    ylim])

subplot(3,1,2)
plot(t,v,    t,Dx/RMS_Dx,'-')
grid
title('Velocity')
axis([0  1    ylim])

subplot(3,1,3)
plot(t,a,    t,D2x/RMS_D2x,'-')
grid
title('Acceleration')
xlabel('Time (sec)')
axis([0  1    ylim])

the acceleration signal is a decent fit but the velocity is not. Since acceleration is calculated from velocity (and both from position), I find that mystifying. My only conclusion is that the ‘target’ velocity is not correct, for some reason.

[SMY]

&#1Thank you Star. I should check all one by one. I will update to you again to see how this gonna be happen. I really appreciated your time and kindness to help me

[Star Strider]

My pleasure.

You need to check the velocity data. Something is wrong.

[SMY]

&#1I have a good news to share with you which I already solved the problem. I used a wrong variable previously that’s why I cannot achieved the target that I want to.From here I learned something. Every coding need to check line by line carefully and never underestimate anything. Because coding is so sensitive things. :D &#1

&#1 Thank a lot Star for your time, your kindness. I really appreciated. Talking to people can give a lot benefits.

[Star Strider]

As always, my pleasure!

I’m glad you got it sorted!

‘Every coding need to check line by line carefully and never underestimate anything.’

I could not possibly have said it better! It brings to mind a poem posted on the Stanford University Computer Bulletin Board in the late 1970s:

I really hate this damned machine, 
I wish that they would sell it!  
It won’t do what I want it to, 
But only what I tell it!