本文想借此表達一下個人的一個觀點，做開發(fā)如果遇到無法解決的難題，可以試著從數序的角度出發(fā)，看能否找到答案。

注：文中配圖只為閱讀輕松一點，本人數學也是半吊子，有錯誤幫忙指正。

是個啥坑？

一個項目中用到一個傳感器測量一物理量，這里假定測量溫度吧。需要判斷其變化趨勢，利用這個變化趨勢去做一些應用。

那么要怎么判斷一個物理量的變化趨勢呢？我們能自然能想到去求取該隨機序列的變化率。這里涉及到一些數序定義。隨機序列有很多可能的來源，最為常見是我之前在<<模數轉換知多少>>中介紹的模數采樣。

這樣將S(t)信號轉換為離散信號序列S(n)，那么對于當前時刻其斜率怎么求取呢？(這里忽略中間的過度態(tài)，僅將其看為線段相連，當然現實應用中如果有更高要求，可以做曲線擬合)

但是如果只判斷，斜率極容易誤判，比如下面這樣的情況：

函數的凹凸性

凹函數

凹函數是一個定義在某個向量空間的凸集C（區(qū)間）上的實值函數f。設f為定義在區(qū)間I上的函數，若對I上的任意兩點x1<x2和任意的實數t屬于(0,1),總有,

則稱函數f為l上凹函數，有的書上也稱為下凸函數。

如果把上述條件中的“≥”改成“>”，則叫做嚴格上凹函數，或叫做嚴格下凸函數。

上面是一維函數情況，這里來個2維函數的圖，剛方便理解

凸函數

設f為定義在區(qū)間I上的函數，若對I上的任意兩點x1<x2和任意的實數t屬于(0,1)，上面不等式變成大于等于，則在該區(qū)間為凸函數。

可見，凹凸是相對的，如f(x)在某區(qū)間為凹，則-f(x)則在該區(qū)間為凸。

性質

• 若一個函數在某區(qū)間二階可導且大于0，則函數在該區(qū)間為凹函數
• 若一個函數在某區(qū)間二階可導且小于0，則函數在該區(qū)間為凸函數

證明，這里就不推導了，可以利用拉格朗日中值定理可以推導出上面這個性質。

來看一下會動的圖,加深一下理解:

函數 從 到 切線為藍色，曲線向上凹,綠色表示曲線是向下凹的，紅色表示曲線的拐點。

sin(2x)的一階導數為：

sin(2x)的二階導數為：

裝逼結束，也可能沒裝對～～～

回到坑里

通過上面裝逼，是否可以利用離散序列的求導數來判斷傳感器的變化趨勢。啥？導數？又要開始表演了？

前面說了一階導數是這樣的：

那么二階導數是哪樣捏？

化簡一下：

其中S[n]表示當前測量點，S[n-1]表示前一個測量點，S[n-2]表示前第2個測量點。

上代碼

#include <stdio.h>
#include <math.h>
#include <string.h>
typedef struct _T_2ND_DRV
{
    float xn1;
    float xn2;
}t_2ND_DRV;
typedef struct _T_1ST_DRV
{
    float xn1;
}t_1ST_DRV;

void init_second_derivative(t_2ND_DRV *pSndDrv)
{
    pSndDrv->xn1 = 0;
    pSndDrv->xn2 = 0;
}

float second_derivative(t_2ND_DRV *pSndDrv, float xn,float T)
{
     float result=0.0f;
     if(T<=0)
         return 0x7FBFFFFF; /*非法數據*/
     result = (xn-2*pSndDrv->xn1-pSndDrv->xn2)/T/T;
     pSndDrv->xn2 = pSndDrv->xn1;
     pSndDrv->xn1 = xn;
    
     return result;
}

void init_fisrt_derivative(t_1ST_DRV *p1stDrv)
{
    p1stDrv->xn1 = 0;
}

float fisrt_derivative(t_1ST_DRV *p1stDrv, float xn,float T)
{
     float result=0.0f;
     if(T<=0)
         return 0x7FBFFFFF; /*非法數據*/
     result = (xn-p1stDrv->xn1)/T; 
     p1stDrv->xn1 = xn;
    
     return result;
}
#define PI 3.1415f
#define SAMPLE_RATE 500.0f
#define SAMPLE_T (1/SAMPLE_RATE)
#define SAMPLE_SIZE (100)
int main()
{
    float sim1[SAMPLE_SIZE];
    float sim2[SAMPLE_SIZE];
    float out1[SAMPLE_SIZE];
    float out2[SAMPLE_SIZE];
    t_2ND_DRV sndDrv;
    t_1ST_DRV frtDrv;
    init_fisrt_derivative(&frtDrv);
    init_second_derivative(&sndDrv);
    
    FILE *pFile=fopen("./simulationSin.csv","wt+");
    if(pFile==NULL)
    {
        printf("simulationSin.csv opened failed");
        return -1;
    }
    
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        sim1[i]=10*sin(2*PI*10*i/500);
    } 
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        out1[i]=fisrt_derivative(&frtDrv,sim1[i],SAMPLE_T);
        out2[i]=second_derivative(&sndDrv,sim1[i],SAMPLE_T);
        fprintf(pFile,"%f,%f,%f\n",sim1[i],out1[i],out2[i]);
    }

    fclose(pFile);

    return 0;
}

利用excel生成曲線：

• 一階導數為正時，函數遞增趨勢；
• 一階導數為負時，函數遞減趨勢；
• 二階導數為0時，出現拐點，趨勢改變；此時如果左右兩側的一階導符號相反，則出現極值。
• 二階導數為負時，其一階導數也即原函數斜率規(guī)律單調減，二階導數為正時，其一階導數也即原函數斜率規(guī)律單調增。

再進一步：

一階導數與二階導數結合起來看，就可以看出測量值變化趨勢的趨勢，比如在前1/4周期，此區(qū)間變換趨勢為增，也即一階導數為正，而其二階導數為負，也可以看出遞增的趨勢是逐漸減小到0的。

代碼優(yōu)化

如果只是做定性判斷，上述函數，完全沒必要與采樣周期做除法，只需要考察其增量即可，代碼可優(yōu)化如下：

typedef struct _T_2ND_DRV
{
    float xn1;
    float xn2;
}t_2ND_DRV;
typedef struct _T_1ST_DRV
{
    float xn1;
}t_1ST_DRV;

void init_second_derivative(t_2ND_DRV *pSndDrv)
{
    pSndDrv->xn1 = 0;
    pSndDrv->xn2 = 0;
}

float second_derivative(t_2ND_DRV *pSndDrv, float xn)
{
     float result=0.0f;
     result = xn-2*pSndDrv->xn1-pSndDrv->xn2;
     pSndDrv->xn2 = pSndDrv->xn1;
     pSndDrv->xn1 = xn;
    
     return result;
}

void init_fisrt_derivative(t_1ST_DRV *p1stDrv)
{
    p1stDrv->xn1 = 0;
}

float fisrt_derivative(t_1ST_DRV *p1stDrv, float xn)
{
     float result=0.0f;
     result = xn-p1stDrv->xn1; 
     p1stDrv->xn1 = xn;
    
     return result;
}

意外收獲

這里意外引入一個可能很多人沒注意的知識點NaN，在計算中，NaN代表非數字，是數字數據類型的成員，可以將其解釋為不確定的或無法表示的值，尤其是在浮點運算中。1985年，IEEE 754浮點標準引入了NaN的系統(tǒng)使用，并表示了其他無限量（如無窮大）。

前述函數返回0x7FBFFFFF,也就是表示無窮大。

不同的操作系統(tǒng)和編程語言可能具有NaN的不同字符串表示形式：

nan
 NaN
 NaN%
 NAN
 NaNQ
 NaNS
 qNaN
 sNaN
 1.#SNAN
 1.#QNAN
 -1.#IND

實際上，由于編碼的NaN具有符號，因此通常也可以在NaN的字符串表示中找到它們，例如：

 -NaN
  NaN12345
 -sNaN12300
 -NaN(s1234)

工程應用

這里給出我的建議方案：

將傳感器信號經由電路處理，模數采樣，在進入前級數字濾波器，濾除不必要的噪聲，在進行一階/二階求導。對于一階和二階求導再做一級移動平均濾波，最后在按照上面描述進行判別變化趨勢，則個人認為基本就比較健壯了。實際移動均值濾波長度不宜選擇過長，否則響應就比較滯后了。不能對傳感器的變化趨勢做出實時的判別。加了后級均值濾波器，則會消除由于波形忽上忽下的隨機噪聲干擾影響，使得系統(tǒng)判別更為健壯，實際濾波器長度需根據不同的場合進行調試優(yōu)化?；蛘咭部梢赃x擇別的IIR/FIR濾波器形式實現。

-END-