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.0000
    1.0000    0.5000    0.5000   -3.5000   -8.0000
    1.0000    1.0000    1.0000   -4.5000  -10.0000


Fy =

   -1.0000   -5.0000   -2.0000   -2.0000    5.0000
         0   -1.0000    2.0000    2.0000         0
    1.0000    2.0000    4.0000    4.0000   -2.5000
    1.0000    1.0000    2.0000    2.0000         0

 2.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).