查看: 739|回复: 0
收起左侧

单片机如何能运行如飞?一种高效实现数学函数的方式!

[复制链接]

  离线 

  • TA的每日心情
    奋斗
    2022-6-21 08:23
  • 签到天数: 2 天

    [LV.1]

    发表于 2022-2-13 09:13:42 | 显示全部楼层 |阅读模式

    有人预言,RISC-V或将是继Intel和Arm之后的第三大主流处理器体系。欢迎访问全球首家只专注于RISC-V单片机行业应用的中文网站

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x
    本帖最后由 塞巴斯蒂安 于 2022-2-13 09:13 编辑

    大家好,我是小麦,今天给大家分享一下如何在资源紧张,算力较低的单片机上实现三角函数的算法。
    之前发过一篇关于IQMath的文章,这个是ti公司平台上的一个数学运算库,里面封装了很多高效的数学运算方法。

    例如在不具备浮点运算器的定点处理器使用定点运算,以前写过一篇Q格式的文章,有简单介绍过这些知识。

    那么问题来了,有一个读者朋友的硬件平台无法使用IQMath,但是他要进行一些三角函数的运算,那么该如何自己动手实现呢?
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(1)

    下面我们来简单介绍一下整体的思路吧,因为硬件平台的资源比较紧张;
    • RAM比较少;
    • ROM比较少;
    • CPU处理速度比较慢;

    所以这里比较常用的方法就是通过空间换时间,预先将sincos的值存储到数组中,需要用的时候,访问数组就可以得到具体的数据。这也就是我们经常会提到的查表法


    下面我们来详细介绍一下。

    一、正弦表

    这个正弦函数表达式是这样的,
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(2)
    具体如下图所示:
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(3)
    正弦波
    首先我们来简单分析一下这个波形:
    • 在蓝色框内是一整个周期的波形;
    • 在红色框内是四分之一个周期的波形;

    其实不难发现,我们只要表示出这四分之一个波形的数据,其余剩下的波形都可以通过换算表示出来。


    这样做就大大节省了查表法所需要的空间。

    下面我们来介绍一下具体如何实现;

    首先我们得搞清楚一个点,就是量纲,统一用归一化的形式来做。
    • y的范围是 [-1, 1]
    • x的范围是[0, 2π],当然,x的范围[-π,π]也是没问题的,下面会继续介绍;

    而在实际的程序中,我们是无法这样去做的,这些数值我们期望通过整形类型去访问,所以我们要做到几点:

    • 尽量避免使用浮点运算;
    • 尽量避免除法;
    • 尽量避免乘法;

    所以这里有必要先了解一下Q格式,用左移和右移去代替乘法和除法,提高运算效率;


    对于X轴的数据,于是可以将[0, 2π]细分成 128 ,256,512或者 1024 等等;

    这里我们先细分成1024等份,正如前面提到的,只需要选择前四分之一周期的内容即可;
    1.     #define POINT_NUM  256
    2. #define PI          3.141592f
    3. for (int i = 0; i < POINT_NUM; i++) {

    4.         printf( "[%03d:%1.4f]\t", i , (sin(i*PI/2 / POINT_NUM)));
    5.         if((i+1) % 8 == 0){
    6.             printf("\r\n");
    7.         }
    8.     }
    复制代码
    打印的输出结果如下:
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(4)
    浮点类型的正弦表

    这里我们可以简单取几个特殊点验证一下,发现整体还是可以接受的;
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(5)
    matlab输出的波形
    下一步就是将浮点数据y转化为Q1.15格式哈,
    1. #define POINT_NUM  256
    2. #define PI          3.141592f
    3. printf("sin=============================================\r\n");
    4.     for (int i = 0; i < POINT_NUM; i++) {

    5.         printf("[ %d:\t0x%04X ]", i, (int16_t)(sin(i*PI/2 / POINT_NUM) * 32768));
    6.         if((i+1) % 8 == 0){
    7.             printf("\r\n");
    8.         }
    9.     }
    复制代码
    最终输出结果如下所示;
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(6)
    Q格式正弦表
    二、源码部分

    下面这部分代码是参考STmcsdk中的一个实例,下面我们会依次分析每个部分的作用,整体的代码具体如下所示;
    1. #define SIN_COS_TABLE {\
    2.     0x0000,0x00C9,0x0192,0x025B,0x0324,0x03ED,0x04B6,0x057F,\
    3.     0x0648,0x0711,0x07D9,0x08A2,0x096A,0x0A33,0x0AFB,0x0BC4,\
    4.     0x0C8C,0x0D54,0x0E1C,0x0EE3,0x0FAB,0x1072,0x113A,0x1201,\
    5.     0x12C8,0x138F,0x1455,0x151C,0x15E2,0x16A8,0x176E,0x1833,\
    6.     0x18F9,0x19BE,0x1A82,0x1B47,0x1C0B,0x1CCF,0x1D93,0x1E57,\
    7.     0x1F1A,0x1FDD,0x209F,0x2161,0x2223,0x22E5,0x23A6,0x2467,\
    8.     0x2528,0x25E8,0x26A8,0x2767,0x2826,0x28E5,0x29A3,0x2A61,\
    9.     0x2B1F,0x2BDC,0x2C99,0x2D55,0x2E11,0x2ECC,0x2F87,0x3041,\
    10.     0x30FB,0x31B5,0x326E,0x3326,0x33DF,0x3496,0x354D,0x3604,\
    11.     0x36BA,0x376F,0x3824,0x38D9,0x398C,0x3A40,0x3AF2,0x3BA5,\
    12.     0x3C56,0x3D07,0x3DB8,0x3E68,0x3F17,0x3FC5,0x4073,0x4121,\
    13.     0x41CE,0x427A,0x4325,0x43D0,0x447A,0x4524,0x45CD,0x4675,\
    14.     0x471C,0x47C3,0x4869,0x490F,0x49B4,0x4A58,0x4AFB,0x4B9D,\
    15.     0x4C3F,0x4CE0,0x4D81,0x4E20,0x4EBF,0x4F5D,0x4FFB,0x5097,\
    16.     0x5133,0x51CE,0x5268,0x5302,0x539B,0x5432,0x54C9,0x5560,\
    17.     0x55F5,0x568A,0x571D,0x57B0,0x5842,0x58D3,0x5964,0x59F3,\
    18.     0x5A82,0x5B0F,0x5B9C,0x5C28,0x5CB3,0x5D3E,0x5DC7,0x5E4F,\
    19.     0x5ED7,0x5F5D,0x5FE3,0x6068,0x60EB,0x616E,0x61F0,0x6271,\
    20.     0x62F1,0x6370,0x63EE,0x646C,0x64E8,0x6563,0x65DD,0x6656,\
    21.     0x66CF,0x6746,0x67BC,0x6832,0x68A6,0x6919,0x698B,0x69FD,\
    22.     0x6A6D,0x6ADC,0x6B4A,0x6BB7,0x6C23,0x6C8E,0x6CF8,0x6D61,\
    23.     0x6DC9,0x6E30,0x6E96,0x6EFB,0x6F5E,0x6FC1,0x7022,0x7083,\
    24.     0x70E2,0x7140,0x719D,0x71F9,0x7254,0x72AE,0x7307,0x735E,\
    25.     0x73B5,0x740A,0x745F,0x74B2,0x7504,0x7555,0x75A5,0x75F3,\
    26.     0x7641,0x768D,0x76D8,0x7722,0x776B,0x77B3,0x77FA,0x783F,\
    27.     0x7884,0x78C7,0x7909,0x794A,0x7989,0x79C8,0x7A05,0x7A41,\
    28.     0x7A7C,0x7AB6,0x7AEE,0x7B26,0x7B5C,0x7B91,0x7BC5,0x7BF8,\
    29.     0x7C29,0x7C59,0x7C88,0x7CB6,0x7CE3,0x7D0E,0x7D39,0x7D62,\
    30.     0x7D89,0x7DB0,0x7DD5,0x7DFA,0x7E1D,0x7E3E,0x7E5F,0x7E7E,\
    31.     0x7E9C,0x7EB9,0x7ED5,0x7EEF,0x7F09,0x7F21,0x7F37,0x7F4D,\
    32.     0x7F61,0x7F74,0x7F86,0x7F97,0x7FA6,0x7FB4,0x7FC1,0x7FCD,\
    33.     0x7FD8,0x7FE1,0x7FE9,0x7FF0,0x7FF5,0x7FF9,0x7FFD,0x7FFE}

    34. const int16_t hSin_Cos_Table[256] = SIN_COS_TABLE;

    35. typedef struct
    36. {
    37.   int16_t hCos;
    38.   int16_t hSin;
    39. } Trig_Components;

    40. /**
    41.   * @brief  This function returns cosine and sine functions of the angle fed in
    42.   *         input
    43.   * @param  hAngle: angle in q1.15 format (-1~0.9999)
    44.   * @retval Trig_Components Cos(angle) and Sin(angle) in Trig_Components format
    45.   */
    46. Trig_Components trig_functions( int16_t hAngle )
    47. {
    48. int32_t shindex;
    49. uint16_t uhindex;

    50. Trig_Components Local_Components;

    51. /* 10 bit index computation  */
    52. shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
    53. uhindex = ( uint16_t )shindex;
    54. //uhindex /= ( uint16_t )64;
    55.     uhindex = uhindex >> 6;
    56. /**
    57.   | hAngle    | angle  | std   |
    58.   | (0,16384]   | U0_90  | (0,0.5] |
    59.   | (16384,32767]  | U90_180  | (0.5,0.99]|
    60.   | (-16384,-1]   | U270_360  | (0,-0.5]  |
    61.   | (-16384,-32768] | U180_270  | (-0.5,-1) |
    62. */
    63. //SIN_MASK        0x0300u
    64. switch ( ( uint16_t )( uhindex ) & SIN_MASK )
    65. {         
    66. //0x0200u
    67.    case U0_90:
    68.   Local_Components.hSin =
    69.             hSin_Cos_Table[( uint8_t )( uhindex )];
    70.   Local_Components.hCos =
    71.             hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    72.   break;
    73. //0x0300u
    74.    case U90_180:
    75.   Local_Components.hSin =
    76.             hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    77.   Local_Components.hCos =
    78.             -hSin_Cos_Table[( uint8_t )( uhindex )];
    79.   break;
    80. //0x0000u
    81.    case U180_270:
    82.   Local_Components.hSin =
    83.             -hSin_Cos_Table[( uint8_t )( uhindex )];
    84.   Local_Components.hCos =
    85.             -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    86.   break;
    87. //0x0100u
    88.    case U270_360:
    89.   Local_Components.hSin =  
    90.             -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    91.   Local_Components.hCos =  
    92.             hSin_Cos_Table[( uint8_t )( uhindex )];
    93.   break;
    94.    default:
    95.   break;
    96. }
    97. return ( Local_Components );
    98. }
    复制代码
    由于输入的hAngleQ1.15格式,所以这里可以简单画个图;下面是角度hAngle0x0000~0xFFFF的示意图,如下所示;
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(7)
    角度值
    这里注意,负数是以补码形式进行保存的,正数的补码等于他本身;
    负数的补码是除了符号位外,其他位取反,然后加上1;
    所以可以算一下 0xFFFF表示-1
    0x8000表示 -32768

    因为Q格式中有无符号的范围和带符号的范围,所以这里的hAngle充分利用这个16 bit的数据,并且兼容了传入参数可以是有符号int16或者是无符号uint16,这里比较绕,先看下面这张图片;
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(8)
    有符号和无符号 对比
    上图中;
    • 左边是有符号int16,右边是无符号数uint16
    • 两个圆形分别表示int16uint16的数值范围;
    • 左边绿色框内的波形相对应,橙色框内的波形相对应;

    这里有几点我们要注意一下,无论是有符号和无符号,他们的周期都是相同的;

    • 有符号整数 int16 :-32768 ~ 32765 ,
    • 无符号整数 uint16 :0 ~ 65535,

    所以这两者都使用 65536个数来表示正弦的一个周期,也就是


    这里是比较关键的地方,因此对于 0x8000 这个关键点,有符号和无符号所表示的数值是不同的;
    • 有符号整数 int16 :0x8000 表示为 -32768;
    • 无符号整数 uint16 :0x8000 表示为 32768;

    因此这他们刚好相差了一个周期 65536,所以表示的正弦数值y是相同的,正如上图中蓝色箭头①和②所示。


    三、内部实现

    由于有符号整数 int16 的最高位是符号位,所以这里我们先把它转化成无符号整形;

    前面用 int32类型是为了防止数据溢出,这里加上32768,相当于对正弦波平移了半个周期,所以在下面y和x的映射关系需要根据实际情况来修改;
    1. /* 10 bit index computation  */
    2. shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
    3. uhindex = ( uint16_t )shindex;
    4. //uhindex /= ( uint16_t )64;
    5. uhindex = uhindex >> 6;
    复制代码
    因为前面提高过正弦表的四分之一是256个数据,所以整个正弦周期应该是 1024 个细分数据,那也就是2的10次,就需要 10 bit
    • 10 bit的数据范围是 0~1023;
    • 16 bit的数据范围是 0~65535;

    为了获取有效的高10 bit数据,对数据右移 6 bit,具体如下所示;

    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(9)
    所以,我们又可以得到以下这个数据的范围 0 ~ 10230 ~ 0x400
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(10)
    因此我们在程序中引入四个掩码,作为正弦波形落在哪个象限的标识位,这样也避免了使用除法运算,提高了效率,具体如下所示;
    1. #define SIN_MASK        0x0300u
    2. #define U0_90           0x0200u
    3. #define U90_180         0x0300u                  
    4. #define U180_270        0x0000u
    5. #define U270_360        0x0100u
    复制代码
    其中,U0_90表示 0° ~ 90°,以此类推;

    那为什么是这个映射关系呢?

    0~90°不应该是从 0x000u~0x100u吗?这里我们再简单解释一下;

    前面有一个这样的操作,具体如下;
    1. shindex = ( ( int32_t )32768 + ( int32_t )hAngle );
    2. uhindex = ( uint16_t )shindex;
    复制代码
    这里的hAngle加上32768,相当于加了一个π正弦波形向左移动了半个周期;因此整体的映射关系要和原始的数据对应起来,具体如下所示;
    国外芯片技术交流-单片机如何能运行如飞?一种高效实现数学函数的方式!risc-v单片机中文社区(11)
    最后,既然我们已经知道波形在哪个象限了,就可以根据当前象限和我们正弦表的关系来得到新的波形,这里有中心对称,关于y轴对称,简单做一下变换就可以得到正弦值和余弦值;
    1. //SIN_MASK        0x0300u
    2. switch ( ( uint16_t )( uhindex ) & SIN_MASK )
    3. {         
    4. //0x0200u
    5.    case U0_90:
    6.   Local_Components.hSin =
    7.             hSin_Cos_Table[( uint8_t )( uhindex )];
    8.   Local_Components.hCos =
    9.             hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    10.   break;
    11. //0x0300u
    12.    case U90_180:
    13.   Local_Components.hSin =
    14.             hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    15.   Local_Components.hCos =
    16.             -hSin_Cos_Table[( uint8_t )( uhindex )];
    17.   break;
    18. //0x0000u
    19.    case U180_270:
    20.   Local_Components.hSin =
    21.             -hSin_Cos_Table[( uint8_t )( uhindex )];
    22.   Local_Components.hCos =
    23.             -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    24.   break;
    25. //0x0100u
    26.    case U270_360:
    27.   Local_Components.hSin =  
    28.             -hSin_Cos_Table[( uint8_t )( 0xFFu - ( uint8_t )( uhindex ) )];
    29.   Local_Components.hCos =  
    30.             hSin_Cos_Table[( uint8_t )( uhindex )];
    31.   break;
    32.    default:
    33.   break;
    34. }
    复制代码





    上一篇:国产WiFi 6芯片走向何方?
    下一篇:树莓派要洗牌MCU市场 - 0.7美元的价格、批量直销、2千万颗的锁定产能
    RISCV作者优文
    全球首家只专注于RISC-V单片机行业应用的中文网站
    回复

    使用道具 举报

    高级模式
    B Color Image Link Quote Code Smilies

    本版积分规则

    关闭

    RISC-V单片机中文网上一条 /2 下一条



    版权及免责声明|RISC-V单片机中文网 |网站地图

    GMT+8, 2024-3-29 02:23 , Processed in 0.506648 second(s), 48 queries .

    快速回复 返回顶部 返回列表