來源 | 嵌入式客棧

今天遇到一個(gè)網(wǎng)友問一個(gè)問題，他有一個(gè)傳感器測(cè)量一個(gè)物理量，需要判斷其變化趨勢(shì)，我給了一些建議，這里將這個(gè)建議展開做些深入分析，并分享給大家。

本文想借此表達(dá)一下個(gè)人的一個(gè)觀點(diǎn)，做開發(fā)如果遇到無法解決的難題，可以試著從數(shù)學(xué)的角度出發(fā)，看能否找到答案。

是個(gè)啥坑？

一個(gè)項(xiàng)目中用到一個(gè)傳感器測(cè)量一物理量，這里假定測(cè)量溫度吧。需要判斷其變化趨勢(shì)，利用這個(gè)變化趨勢(shì)去做一些應(yīng)用。

那么要怎么判斷一個(gè)物理量的變化趨勢(shì)呢？我們能自然能想到去求取該隨機(jī)序列的變化率。這里涉及到一些數(shù)序定義。隨機(jī)序列有很多可能的來源，最為常見是模數(shù)采樣。

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

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

其斜率一會(huì)兒正，一會(huì)兒負(fù)，但是其總體趨勢(shì)又是在增加的，所以只考察斜率顯然不可取，獲取需要在代碼在加各種復(fù)雜的條件或者限值去判斷。即使加這么多條件系統(tǒng)仍然可能表現(xiàn)的非常不健壯。

對(duì)于模擬信號(hào)2而言，趨勢(shì)又在不斷變化。那么怎么做才能穩(wěn)定呢？先賣個(gè)關(guān)子？

函數(shù)的凹凸

• 凹函數(shù)
  

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

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

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

上面是一維函數(shù)情況，這里來個(gè)2維函數(shù)的圖，剛方便理解

• 凸函數(shù)

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

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

• 性質(zhì)

若一個(gè)函數(shù)在某區(qū)間二階可導(dǎo)且大于0，則函數(shù)在該區(qū)間為凹函數(shù)

若一個(gè)函數(shù)在某區(qū)間二階可導(dǎo)且小于0，則函數(shù)在該區(qū)間為凸函數(shù)

證明，這里就不推導(dǎo)了，可以利用拉格朗日中值定理可以推導(dǎo)出上面這個(gè)性質(zhì)。

來看一下會(huì)動(dòng)的圖,加深一下理解:

函數(shù)從到切線為藍(lán)色，曲線向上凹,綠色表示曲線是向下凹的，紅色表示曲線的拐點(diǎn)。

sin(2x)的一階導(dǎo)數(shù)為：

sin(2x)的二階導(dǎo)數(shù)為：

裝逼結(jié)束，也可能沒裝對(duì)～～～

• 回到坑里

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

前面說了一階導(dǎo)數(shù)是這樣的：

那么二階導(dǎo)數(shù)是哪樣捏？

化簡一下：

其中S[n]表示當(dāng)前測(cè)量點(diǎn)，S[n-1]表示前一個(gè)測(cè)量點(diǎn)，S[n-2]表示前第2個(gè)測(cè)量點(diǎn)。應(yīng)為+S[n-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;?/*非法數(shù)據(jù)*/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;?/*非法數(shù)據(jù)*/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; }

忽略前兩個(gè)點(diǎn)，利用excel生成曲線：

從圖中可看出：

• 一階導(dǎo)數(shù)為正時(shí)，函數(shù)遞增趨勢(shì)；

• 一階導(dǎo)數(shù)為負(fù)時(shí)，函數(shù)遞減趨勢(shì)；

• 二階導(dǎo)數(shù)為0時(shí)，出現(xiàn)拐點(diǎn)，趨勢(shì)改變；此時(shí)如果左右兩側(cè)的一階導(dǎo)符號(hào)相反，則出現(xiàn)極值。

• 二階導(dǎo)數(shù)為負(fù)時(shí)，其一階導(dǎo)數(shù)也即原函數(shù)斜率規(guī)律單調(diào)減，二階導(dǎo)數(shù)為正時(shí)，其一階導(dǎo)數(shù)也即原函數(shù)斜率規(guī)律單調(diào)增。

再進(jìn)一步：

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

• 代碼優(yōu)化

如果只是做定性判斷，上述函數(shù)，完全沒必要與采樣周期做除法，只需要考察其增量即可，代碼可優(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; }

• 意外收獲

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

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

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

nanNaNNaN%NANNaNQNaNSqNaNsNaN1.#SNAN1.#QNAN-1.#IND

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

?-NaNNaN12345-sNaN12300-NaN(s1234)

• 工程應(yīng)用

這里給出我的建議方案：

將傳感器信號(hào)經(jīng)由電路處理，模數(shù)采樣，在進(jìn)入前級(jí)數(shù)字濾波器，濾除不必要的噪聲，在進(jìn)行一階/二階求導(dǎo)。對(duì)于一階和二階求導(dǎo)再做一級(jí)移動(dòng)平均濾波，最后在按照上面描述進(jìn)行判別變化趨勢(shì)，則個(gè)人認(rèn)為基本就比較健壯了。實(shí)際移動(dòng)均值濾波長度不宜選擇過長，否則響應(yīng)就比較滯后了。不能對(duì)傳感器的變化趨勢(shì)做出實(shí)時(shí)的判別。加了后級(jí)均值濾波器，則會(huì)消除由于波形忽上忽下的隨機(jī)噪聲干擾影響，使得系統(tǒng)判別更為健壯，實(shí)際濾波器長度需根據(jù)不同的場(chǎng)合進(jìn)行調(diào)試優(yōu)化。或者也可以選擇別的IIR/FIR濾波器形式實(shí)現(xiàn)。

總結(jié)

做為嵌入式er編程，有時(shí)候有必要去看看數(shù)學(xué)書，了解一下數(shù)學(xué)原理的背后故事，可能會(huì)給你帶來意想不到的作用哦。

推薦閱讀

• 一行代碼引來的安全漏洞，就讓我們丟失了整個(gè)服務(wù)器的控制權(quán)

• 怒肝 8 個(gè)月源碼，我成為了 Spring 開源貢獻(xiàn)者

• 讀取模式錯(cuò)誤，計(jì)算引擎操作復(fù)雜……面對(duì)Hadoop這些問題該如何應(yīng)對(duì)？
  

• 從地?cái)偪丛朴?jì)算：規(guī)模產(chǎn)業(yè)歷程大揭秘

• 這位超級(jí)電腦之父，年近九旬仍不愿退休
  

• 阿里云科學(xué)家入選計(jì)算機(jī)頂會(huì)HPCA名人堂，他是什么來頭？
  

• 中本聰真的轉(zhuǎn)移了2009年挖出的比特幣嗎？或許不一定！
  

真香，朕在看了！

總結(jié)

以上是生活随笔為你收集整理的数学之美：嵌入式编程凹凸性之妙用（附C代码）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。