本文想借此表达一下个人的一个观点，做开发如果遇到无法解决的难题，可以试着从数学的角度出发，看能否找到答案。

注：文中配图只为阅读轻松一点，本人数学也是半吊子，有错误帮忙指正。

4 j# n& G0 \0 y# Z/ \- M

是个啥坑？

一个项目中用到一个传感器测量一物理量，这里假定测量温度吧。需要判断其变化趋势，利用这个变化趋势去做一些应用。

那么要怎么判断一个物理量的变化趋势呢？我们能自然能想到去求取该随机序列的变化率。这里涉及到一些数序定义。随机序列有很多可能的来源，最为常见是我之前在<<模数转换知多少>>中介绍的模数采样。

这样将S(t)信号转换为离散信号序列S(n)，那么对于当前时刻其斜率怎么求取呢？(这里忽略中间的过度态，仅将其看为线段相连，当然现实应用中如果有更高要求，可以做曲线拟合)

但是如果只判断，斜率极容易误判，比如下面这样的情况：

其斜率一会儿正，一会儿负，但是其总体趋势又是在增加的，所以只考察斜率显然不可取，获取需要在代码在加各种复杂的条件或者限值去判断。即使加这么多条件系统仍然可能表现的非常不健壮。

对于模拟信号2而言，趋势又在不断变化。那么怎么做才能稳定呢？先卖个关子？

函数的凹凸性凹函数

9 X  Z& y  `3 G# B/ g9 A

凹函数是一个定义在某个向量空间的凸集C（区间）上的实值函数f。设f为定义在区间I上的函数，若对I上的任意两点x1<x2和任意的实数t属于(0,1),总有,

3 ~. z- U* }7 _4 M, ~) I

则称函数f为l上凹函数，有的书上也称为下凸函数。

如果把上述条件中的“≥”改成“>”，则叫做严格上凹函数，或叫做严格下凸函数。

上面是一维函数情况，这里来个2维函数的图，刚方便理解

凸函数
) c, m3 ?9 c5 H, ]2 h9 R% A' B9 C

设f为定义在区间I上的函数，若对I上的任意两点x1<x2和任意的实数t属于(0,1)，上面不等式变成大于等于，则在该区间为凸函数。

可见，凹凸是相对的，如f(x)在某区间为凹，则-f(x)则在该区间为凸。

性质

• 若一个函数在某区间二阶可导且大于0，则函数在该区间为凹函数
• 若一个函数在某区间二阶可导且小于0，则函数在该区间为凸函数
  

证明，这里就不推导了，可以利用拉格朗日中值定理可以推导出上面这个性质。

来看一下会动的图,加深一下理解:

函数从到切线为蓝色，曲线向上凹,绿色表示曲线是向下凹的，红色表示曲线的拐点。

sin(2x)的一阶导数为：

sin(2x)的二阶导数为：

装逼结束，也可能没装对～～～

回到坑里

通过上面装逼，是否可以利用离散序列的求导数来判断传感器的变化趋势。啥？导数？又要开始表演了？

前面说了一阶导数是这样的：

  T" [% T5 f3 R# X

那么二阶导数是哪样捏？

化简一下：

其中S[n]表示当前测量点，S[n-1]表示前一个测量点，S[n-2]表示前第2个测量点。应为+S[n-2]

: G# C0 ]" t  T" }) K* L* n) k

上代码
#include <stdio.h>

3 a4 M# l# B, j& N- ~; l2 m2 }& g#include <math.h>
#include <string.h>
6 @5 C6 i6 P3 I5 D: q. \typedef struct _T_2ND_DRV
% y* w5 O2 E% x! Y( ^{
    float xn1;
    float xn2;
}t_2ND_DRV;
/ {/ V  q- A6 |& h) Jtypedef struct _T_1ST_DRV
+ W  H% D( _( J$ X5 G! K" a. I{
    float xn1;
' |" k' R5 H1 W& i' N- e4 D0 X}t_1ST_DRV;
& u! o5 O4 l1 W6 y9 g
& x% _8 G( H; s2 C* Vvoid init_second_derivative(t_2ND_DRV *pSndDrv)
{
    pSndDrv->xn1 = 0;
    pSndDrv->xn2 = 0;
}
+ X2 O" Q1 q" a
$ ?9 D  T% S! O. Kfloat second_derivative(t_2ND_DRV *pSndDrv, float xn,float T)
0 V% n3 @0 s" M4 z1 M. t' h+ l{
# i4 Y) u- x& x$ Z     float result=0.0f;
     if(T<=0)
/ V( S- t7 y) [% x         return 0x7FBFFFFF; /*非法数据*/
     result = (xn-2*pSndDrv->xn1+pSndDrv->xn2)/T/T;
$ O  o. v+ C0 [     pSndDrv->xn2 = pSndDrv->xn1;
/ |8 [( t4 C# Y/ ~: B     pSndDrv->xn1 = xn;
  X! d( c9 }# O3 i   
+ e4 n+ q. Y. c# @3 Y4 a     return result;
. B  ~' M- c2 B& z( o# Y+ `  O}
! h- y( |. Z& Y& l  ~' Q' @; W
void init_fisrt_derivative(t_1ST_DRV *p1stDrv)
{
    p1stDrv->xn1 = 0;
$ K- @" J( H# h  B}

float fisrt_derivative(t_1ST_DRV *p1stDrv, float xn,float T)
/ I' F8 M1 m. R+ }5 e$ }3 B* g{
& _8 J# b! `) t8 j1 _8 F: }: W     float result=0.0f;
1 v" A; k8 `3 D% v     if(T<=0)
         return 0x7FBFFFFF; /*非法数据*/
% y' L4 R( d, }! }; T; f8 n" f  B     result = (xn-p1stDrv->xn1)/T;
     p1stDrv->xn1 = xn;
   
     return result;
) a$ U( C* B% s* q' l- z  `}
" D$ c: o" X; R1 `; J- ?#define PI 3.1415f
#define SAMPLE_RATE 500.0f
" v8 q# w2 k0 d( ~! A% u#define SAMPLE_T (1/SAMPLE_RATE)
; Q4 {9 y1 e! \$ Z: n  J  [#define SAMPLE_SIZE (100)
int main()
{
    float sim1[SAMPLE_SIZE];
    float sim2[SAMPLE_SIZE];
$ P- x- W5 w" r$ i6 i    float out1[SAMPLE_SIZE];
0 R. ~! z8 d0 H) S$ Y! B2 b0 X    float out2[SAMPLE_SIZE];
" _: D, E/ H; {" Q9 X    t_2ND_DRV sndDrv;
    t_1ST_DRV frtDrv;
    init_fisrt_derivative(&frtDrv);
    init_second_derivative(&sndDrv);
6 N- s1 D4 a4 M8 E% D% c   
: {: h. P  n( d( b! I    FILE *pFile=fopen("./simulationSin.csv","wt+");
    if(pFile==NULL)
$ S& q7 }* A$ E& T9 T    {
        printf("simulationSin.csv opened failed");
3 r/ r: z) o8 y/ c        return -1;
" d( J  `4 W8 @# X4 O  g    }
   
* j, \; \7 z$ [1 c- q$ [( |    for(int i=0;i<SAMPLE_SIZE;i++)
    {
2 A; c  H) q( q. H        sim1=10*sin(2*PI*10*i/500);
    }
    for(int i=0;i<SAMPLE_SIZE;i++)
. P: w4 S, D8 c* r6 L* I    {
        out1=fisrt_derivative(&frtDrv,sim1,SAMPLE_T);
3 [: p- Y, I5 ~3 s4 n6 D" R        out2=second_derivative(&sndDrv,sim1,SAMPLE_T);
        fprintf(pFile,"%f,%f,%f\n",sim1,out1,out2);
    }
# ?, N, B0 E2 l* `) ^- A
; s: d7 g0 s6 {9 m  K  V2 ^    fclose(pFile);
! V. m( A" s/ ]/ [' o
    return 0;
' t" i! V' Z; K0 T}
* y( z- N% H/ H( k

忽略前两个点，利用excel生成曲线：

从图中可看出：

• 一阶导数为正时，函数递增趋势；
• 一阶导数为负时，函数递减趋势；
• 二阶导数为0时，出现拐点，趋势改变；此时如果左右两侧的一阶导符号相反，则出现极值。
• 二阶导数为负时，其一阶导数也即原函数斜率规律单调减，二阶导数为正时，其一阶导数也即原函数斜率规律单调增。
  

再进一步：

一阶导数与二阶导数结合起来看，就可以看出测量值变化趋势的趋势，比如在前1/4周期，此区间变换趋势为增，也即一阶导数为正，而其二阶导数为负，也可以看出递增的趋势是逐渐减小到0的。

代码优化

如果只是做定性判断，上述函数，完全没必要与采样周期做除法，只需要考察其增量即可，代码可优化如下：

typedef struct _T_2ND_DRV
{
3 I! V; y- E. a    float xn1;
# J9 ^& G  p, Q) @# S2 d    float xn2;
( H' H5 N% @" Y. N* R}t_2ND_DRV;
typedef struct _T_1ST_DRV
. y* @+ X7 m5 c7 C5 l0 {# e{
" J) v$ u& ?. z2 k: B    float xn1;
& r# |& f8 v. g/ a6 O}t_1ST_DRV;
% P/ u% L0 O, w
void init_second_derivative(t_2ND_DRV *pSndDrv)
{
2 \2 }/ e( ^, V: E% a1 \. F    pSndDrv->xn1 = 0;
2 w( n7 i9 g7 X2 `; Q4 \    pSndDrv->xn2 = 0;
& X+ ~+ I. O. ]5 M1 q}

float second_derivative(t_2ND_DRV *pSndDrv, float xn)
# ^" [. x- A& @9 l{
     float result=0.0f;
     result = xn-2*pSndDrv->xn1+pSndDrv->xn2;
     pSndDrv->xn2 = pSndDrv->xn1;
     pSndDrv->xn1 = xn;
) A, P& I) U% ^  \" i9 K/ o( {   
     return result;
2 g0 f2 M1 F; X. |( \0 [# D; w2 C7 n}
# D% }9 n) r* P7 K4 l
void init_fisrt_derivative(t_1ST_DRV *p1stDrv)
{
% Y0 e/ T9 j1 |1 o    p1stDrv->xn1 = 0;
}

float fisrt_derivative(t_1ST_DRV *p1stDrv, float xn)
6 c- H  j" ]* F) v! a{
& V1 I9 p6 T2 ^  h+ S5 m) B8 \     float result=0.0f;
     result = xn-p1stDrv->xn1;
     p1stDrv->xn1 = xn;
3 S/ {: ?4 H( P   
& n* `# y* j7 f; @% S1 q     return result;
; z: w- a6 x; B* |# z( d( p}
意外收获

这里意外引入一个可能很多人没注意的知识点NaN，在计算中，NaN代表非数字，是数字数据类型的成员，可以将其解释为不确定的或无法表示的值，尤其是在浮点运算中。1985年，IEEE 754浮点标准引入了NaN的系统使用，并表示了其他无限量（如无穷大）。

前述函数返回0x7FBFFFFF,也就是表示无穷大。

不同的操作系统和编程语言可能具有NaN的不同字符串表示形式：

nan
NaN
# k) U( b  l$ e$ A" w! ?1 y* n* D NaN%
NAN
, ~$ }. @4 h% W1 ]% |, T" L NaNQ
NaNS
qNaN
sNaN
1.#SNAN
/ f0 ~" }9 a- e4 w; H( L- E6 Y 1.#QNAN
-1.#IND

实际上，由于编码的NaN具有符号，因此通常也可以在NaN的字符串表示中找到它们，例如：

-NaN
  NaN12345
-sNaN12300
-NaN(s1234)
工程应用

这里给出我的建议方案：

将传感器信号经由电路处理，模数采样，在进入前级数字滤波器，滤除不必要的噪声，在进行一阶/二阶求导。对于一阶和二阶求导再做一级移动平均滤波，最后在按照上面描述进行判别变化趋势，则个人认为基本就比较健壮了。实际移动均值滤波长度不宜选择过长，否则响应就比较滞后了。不能对传感器的变化趋势做出实时的判别。加了后级均值滤波器，则会消除由于波形忽上忽下的随机噪声干扰影响，使得系统判别更为健壮，实际滤波器长度需根据不同的场合进行调试优化。或者也可以选择别的IIR/FIR滤波器形式实现。

具体实现可参考(点击可阅读)：
6 H  ~) q) d) f% |! B

[color=var(--weui-LINK)]手把手教系列之移动平均滤波器C实现
5 e$ t# M- R) x

[color=var(--weui-LINK)]手把手教系列之IIR数字滤波器设计实现
: U- a3 q2 Z1 F! X9 w0 p

手把手教你系列之FIR滤波器设计实现
; N) H  u2 ^  t8 p+ ~) ]3 C* s+ {

总结一下

做为嵌入式er编程，有时候有必要去看看数学书，了解一下数学原理的背后故事，可能会给你带来意想不到的作用哦。

9 |' x9 f$ _1 G

谢谢分享