zoukankan      html  css  js  c++  java
  • 数学之美:嵌入式编程凹凸性之妙用(附C代码)

    注、星标嵌入式客栈,干货及时送达

    [导读] 咦,你已被成功吸引进来了,不是你想的那样哈~~~

    皮一下哈,言归正传,今天遇到一个网友问一个问题,他有一个传感器测量一个物理量,需要判断其变化趋势,我给了一些建议,这里将这个建议展开做些深入分析,并分享给大家。

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

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

    是个啥坑?

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

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

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

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

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

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

    函数的凹凸性

    凹函数

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

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

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

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

    凸函数

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

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

    性质

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

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

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

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

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

    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
    ",sim1[i],out1[i],out2[i]);
        }
    
        fclose(pFile);
    
        return 0;
    }
    

    利用excel生成曲线:

    从图中可看出:
    • 一阶导数为正时,函数递增趋势;

    • 一阶导数为负时,函数递减趋势;

    • 二阶导数为0时,出现拐点,趋势改变;此时如果左右两侧的一阶导符号相反,则出现极值。

    • 二阶导数为负时,其一阶导数也即原函数斜率规律单调减,二阶导数为正时,其一阶导数也即原函数斜率规律单调增。

    再进一步:

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

    代码优化

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

    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的系统使用,并表示了其他无限量(如无穷大)。

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

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

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

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

     -NaN
      NaN12345
     -sNaN12300
     -NaN(s1234)
    

    工程应用

    这里给出我的建议方案:

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

    具体实现可参考(点击可阅读):

    手把手教系列之移动平均滤波器C实现

    手把手教系列之IIR数字滤波器设计实现

    手把手教你系列之FIR滤波器设计实现

    总结一下

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

    本文辛苦原创,也是为了帮助网友,如有错误之处,也请联系我指出。希望大家不吝点个赞,随心打个赏也是可以的,哈哈。也请帮忙转发分享支持一下~~~

    END

    往期精彩推荐,点击即可阅读

    ▲Linux内核中I2C总线及设备长啥样? [墙裂推荐]

    学习AI之机器学习概念篇

    手把手教系列之IIR数字滤波器设计实现

    为节约时间,文章同步自公众号(首发) 如需学习资料,请用微信关注文中公众号二维码,后台发送"领取",可免费获取海量学习资料,涵盖单片机技术,信号处理,人工智能,嵌入式linux,C/C++编程,数据结构与算法
  • 相关阅读:
    统一Windows Azure和一般web应用之间的文件操作代码(转+译)
    Windows Azure真实案例Lokad 公司通过软件+服务高效提供先进的预测服务
    Windows Azure Marketplace入门教学利用TabLeau Public构建可视化DataMarket应用
    SQL Azure 一款强大的管理工具 Houston CTP 1(转+译)
    Windows Azure真实案例:Invensys Operations Management 公司使用Windows Azure AppFabric 实现动态节能的智能电网
    开始Azure之旅,参加深度培训 (转)
    Windows Azure Marketplace入门教学通过代码操作DataMarket数据源
    Windows Azure真实案例:Infosys Technologies 使用SQL Data Services(现为SQL Azure)为汽车经销商创建了基于云的方案
    Windows Azure Marketplace入门教学 DataMarket for Excel插件
    Windows Azure AppFabric Caching入门简介
  • 原文地址:https://www.cnblogs.com/embInn/p/14038134.html
Copyright © 2011-2022 走看看