The del2 function can do this, but you have to be careful. Differentiation significantly amplifies noise.
How I can produce the double differentiation of digital signal data without knowing the function?
1 次查看（过去 30 天）
显示 更早的评论
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
采纳的回答
Star Strider
2016-7-8
22 个评论
SMY
2016-7-8
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
2016-7-8
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
2016-7-8
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
SMY
2016-7-18
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
2016-7-18
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
2016-7-18
Thank you Star.
Can I have your email for further discussion. I really needs help. I really3 needs help.
Star Strider
2016-7-18
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
2016-7-18
编辑：SMY
2016-7-18
Thank you Star,
I upload my coding (matlabwork.m) for your references, comments and help. I don't know how to improve it anymore
It 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
I really appreciated your time
Regards
SMY
Star Strider
2016-7-19
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
2016-7-19
编辑：SMY
2016-7-19
Thank you Star
Maybe 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?

Star Strider
2016-7-19
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
2016-7-19
Thank you Star for the suggestion
But still work to get velocity but not for acceleration. I feeling bad seriously.
If you dont mind, can you show me how to get the Sine Wave, Square Wave and triangular wave based on my data?
Star Strider
2016-7-19
‘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
2016-7-19
What 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. 
I 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.
Star Strider
2016-7-19
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
2016-7-19
 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
Star Strider
2016-7-20
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
2016-7-20
Thank 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
SMY
2016-7-20
I 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 
 Thank a lot Star for your time, your kindness. I really appreciated. Talking to people can give a lot benefits.
Star Strider
2016-7-20
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!
更多回答（0 个）
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)