1.通过差分估算
已知同维度的x和y序列，则可使用diff(y)./diff(x)来估算。设x为n维向量，Dx=diff(x) 计算向量x的向前差分，DX(i)=X(i+1)-X(i)，0<i<n。
例一
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...0.68 0.30];
x=0:0.04:1.48;
plot(x,y)
dy=diff(y)./diff(x);
dx=0.04:0.04:1.48;
figure
plot(dx,dy)
2.通过梯度估算
已知同维度的x和y序列，则可使用gradient(y)./gradient(x)来估算 。梯度用的是中心点差分，diff()用的前后两点差分；所以从区间上看梯度用的范围比导数大一倍！所以梯度方式精度会更高一些！但是梯度法的边界可能会出现误差问题。
2.1梯度的含义与使用
[Fx,Fy]=gradient(F)，其中Fx为其水平方向上的梯度，Fy为其垂直方向上的梯度，Fx的第一列元素为原矩阵第二列与第一列元素之差，Fx的第二列元素为原矩阵第三列与第一列元素之差除以2，以此类推：Fx(i,j)=(F(i,j+1)-F(i,j-1))/2。最后一列则为最后两列之差。同理，可以得到Fy。
例二
x=[6,9,3,4,0;5,4,1,2,5;6,7,7,8,0;7,8,9,10,0];
[Fx,Fy]=gradient(x)运行结果：
Fx =3.0000 -1.5000 -2.5000 -1.5000 -4.0000-1.0000 -2.0000 -1.0000 2.0000 3.00001.0000 0.5000 0.5000 -3.5000 -8.00001.0000 1.0000 1.0000 -4.5000 -10.0000Fy =-1.0000 -5.0000 -2.0000 -2.0000 5.00000 -1.0000 2.0000 2.0000 01.0000 2.0000 4.0000 4.0000 -2.50001.0000 1.0000 2.0000 2.0000 02.2通过梯度估算
例三：差分与梯度对比
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...0.68 0.30];
x=0:0.04:1.48;
dy=gradient(y)./gradient(x);
dx=0:0.04:1.48;
dx1=0.04:0.04:1.48;
figure
plot(dx,dy,'r')
dy1=diff(y)./diff(x);
hold on
plot(dx1,dy1,'b')
legend('梯度','差分')运行结果
3. 数据拟合后求导
例四
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...0.68 0.30];
x=0:0.04:1.48;
p=polyfit(x,y,10); %使用10次多项式拟合
y3=polyval(p,x); %求出预测值
plot(x,y,'b',x,y3,'r')
legend('原始数据','拟合函数')
dy=diff(y)./diff(x);
dx=0.04:0.04:1.48;
% dy3=diff(y3)./diff(x);
p1=polyder(p);
dy3=polyval(p1,dx);
figure
plot(dx,dy,'b',dx,dy3,'r')
legend('差分求导','拟合多项式求导')运行结果
4.3点、4点、5点求导法
例五：使用5点求导法[5]
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...0.68 0.30]';
x=(0:0.04:1.48)';
y1 = fivePoint1Order(y,0.04);
figure
plot(x(2:end,1),y1,'r')
function [value] = fivePoint1Order(Sig,h)
% 使用五点法求一阶导数
% Sig被求导的列向量
% h步长，单位为秒
lengSig = size(Sig,1);
value = zeros(lengSig-1,1);
% 边缘点使用差分，参考matlab梯度函数的做法
value(1,1) = (Sig(2,1) - Sig(1,1))/h;for i = 1:lengSig-4
fivePoints = Sig(i:i+4,1);
f_2 = fivePoints(1,1);
f_1 = fivePoints(2,1);
f1 = fivePoints(4,1);
f2 = fivePoints(5,1);
value(i+1,1) = (f_2 - 8*f_1 + 8*f1 - f2)/(12*h);
end
value(i+2,1) = (Sig(i+3,1) - Sig(i+2,1))/h;
value(i+3,1) = (Sig(i+4,1) - Sig(i+3,1))/h;% 边缘点
end运行结果
三点、四点求导法请参考https://blog.csdn.net/weixin_30565327/article/details/95006050
参考文献
[1]https://zhidao.baidu.com/question/406266767.html
[2]http://blog.sina.com.cn/s/blog_53be544e0101cd8a.html
[3]https://www.jianshu.com/p/eb8f64a4bec4
[4]https://jingyan.baidu.com/article/59a015e3586c7ff79488650a.html
[5]RAFATI FARD M, SHARIAT MOHAYMANY A, SHAHRI M. A new methodology for vehicle trajectory reconstruction based on wavelet analysis [J]. Transportation Research Part C: Emerging Technologies, 2017, 74(150-67).
