在嵌入式设备中用多项式快速计算三角函数和方根

惯性传感器的倾角计算要用到三角函数.

在 MCS-51, Cortex M0, M3 之类的芯片上编程时, 能使用的资源是非常有限, 通常只有两位数KB的Flash, 个位数KB的RAM. 如果要使用三角函数和开方就要引入 math.h, 会消耗掉10KB以上的Flash空间. 在很多情况下受硬件资源限制无法使用 math.h, 这时候使用简化的方法进行三角函数和开方运算就非常有意义, OlliW’s Bastelseiten在2014年的一篇文章里, 提供了几个实用的计算方法. 下面介绍其计算方法和代码实现.

快速正弦余弦(Sin, Cos)计算

将角度 x ∈ [ 0 , π 2 ] x \in [0, \frac{\pi}{2}] x[0,2π]通过下面的式子转换到 $ \alpha \in [-\frac{1}{2}, \frac{1}{2}]$ 区间

α = 2 π x − 1 2 \alpha = \frac{2}{\pi} x - \frac{1}{2} α=π2x21

于是, 对应 α \alpha α 的多项式近似计算为
sin ⁡ α = a 0 − b 1 α + a 2 α 2 − b 3 α 3 + a 4 α 4 − b 5 α 5 + a 6 α 6 cos ⁡ α = a 0 + b 1 α + a 2 α 2 + b 3 α 3 + a 4 α 4 + b 5 α 5 + a 6 α 6 \sin\alpha = a_0 - b_1\alpha + a_2\alpha^2 - b_3\alpha^3 + a_4\alpha^4 - b_5\alpha^5 + a_6\alpha^6 \\ \cos\alpha = a_0 + b_1\alpha + a_2\alpha^2 + b_3\alpha^3 + a_4\alpha^4 + b_5\alpha^5 + a_6\alpha^6 sinα=a0b1α+a2α2b3α3+a4α4b5α5+a6α6cosα=a0+b1α+a2α2+b3α3+a4α4+b5α5+a6α6
如果将上面的符号固定项和变化项分成 A A A B B B两部分
A = a 0 + a 2 α 2 + a 4 α 4 + a 6 α 6 B = b 1 α + b 3 α 3 + b 5 α 5 A = a_0 + a_2\alpha^2 + a_4\alpha^4 + a_6\alpha^6 \\ B = b_1\alpha + b_3\alpha^3 + b_5\alpha^5 A=a0+a2α2+a4α4+a6α6B=b1α+b3α3+b5α5
sin ⁡ α \sin\alpha sinα cos ⁡ α \cos\alpha cosα 可以通过 A 和 B 的值表达
sin ⁡ α = A − B cos ⁡ α = A + B \sin\alpha = A - B \\ \cos\alpha = A + B sinα=ABcosα=A+B

对应的各项系数值

a 0 = 0.707106781187 a 2 = − 0.872348075361 a 4 = 0.179251759526 a 6 = − 0.0142718282624 b 1 = − 1.110670322264 b 3 = 0.4561589075945 b 5 = − 0.0539104694791 a_0 = 0.707106781187 \\ a_2 = -0.872348075361 \\ a_4 = 0.179251759526 \\ a_6 = -0.0142718282624 \\ \\ b_1 = -1.110670322264 \\ b_3 = 0.4561589075945 \\ b_5 = -0.0539104694791 a0=0.707106781187a2=0.872348075361a4=0.179251759526a6=0.0142718282624b1=1.110670322264b3=0.4561589075945b5=0.0539104694791

使用上面的计算方式, 结果绝对误差小于 6.5 × 1 0 − 6 6.5 \times 10^{-6} 6.5×106, 并且 cos ⁡ 2 x + sin ⁡ 2 x \cos^2 x + \sin^2 x cos2x+sin2x 不会超过 1. 计算过程只需要7次乘法和7次加法.

C语言实现

const   float coeff[7] = {// a0 ~ a6           b1 ~ b50.707106781187,  -1.110670322264,-0.872348075361,   0.4561589075945,0.179251759526,  -0.0539104694791,-0.0142718282624
};/*** @param alpha: value between 0 and 0.5
*/
void sincos_normalized(float alpha, float *sin, float *cos)
{int i;float alpha_exp = 1.0, part_a = 0, part_b = 0;for (i = 0; i < 7; i++){if (i % 2 == 0){part_a = part_a + (coeff[i] * alpha_exp);}else{part_b = part_b + (coeff[i] * alpha_exp);}alpha_exp = alpha_exp * alpha;}*sin = part_a - part_b;*cos = part_a + part_b;
}float calculate(float degree_in)
{int quadrant, multi;float degree = degree_in, alpha, cos, sin, c, s;multi = (int)(degree / 90.0);degree = degree - (multi * 90.0);alpha = (degree / 90) - 0.5;sincos_normalized(alpha, &s, &c);multi = multi % 4;if (multi == 0){sin = s;cos = c;}else if (multi == 1){sin = c;cos = -s;}else if (multi == 2){sin = -s;cos = -c;}else if (multi == 3){sin = -c;cos = s;}printf("d_in:%5.0f d:%5.0f a:%10.5f  sin:%10.5f  cos:%10.5f\r\n", degree_in, degree, alpha, sin, cos);
}

计算的结果和 math.h 的 sin cos 函数对比, 数值几乎一样, 仅在个别数值的小数点后第五位会有 ± 1 \pm1 ±1的差异.

平方根倒数计算

对于1附近的数值, 平方根倒数可以使用牛顿迭代法计算, 实际上非常简单,因为它只涉及加法和乘法,而不涉及除法, 对于 x ∈ [ 0.6 , 1.4 ] x \in [0.6, 1.4] x[0.6,1.4], 计算式为
y 0 = 1 y n + 1 = y n ( 1.5 − 0.5 x y n 2 ) y_0 = 1 \\ y_{n+1} = y_n (1.5 - 0.5 x {y_n}^2) \\ y0=1yn+1=yn(1.50.5xyn2)

计算两次牛顿迭代需要3次乘法, 而二阶泰勒级数只需要2次, 但是牛顿迭代法精度更高, 甚至比三阶泰勒级数的精度更高. 如果执行三次牛顿迭代则需要6次乘法, 在 0.6 < x < 1.4 0.6 < x < 1.4 0.6<x<1.4的范围内结果精度优于 1 × 1 0 − 4 1 \times 10^{-4} 1×104, 注意 x x x的取值范围, 因为近似是以1为中心展开的, 所以离1越远差异越大, 在这个范围之外例如 x = 0.5 x = 0.5 x=0.5的时候, 三次迭代就达不到这个精度. 在实际应用中, 可以将要计算的数值提一个系数转换到 x ∈ [ 0.6 , 1.4 ] x \in [0.6, 1.4] x[0.6,1.4]区间进行计算.

C语言实现

float inverse_sqrt(int interates, float x)
{float y;y = 1.5 - (x / 2);while (interates--){y = y * (1.5 - 0.5 * x * y * y);}return y;
}// 使用 0.5 ~ 2.1 之间的数字测试, 分别迭代5次
int main(int argc, char *const argv[])
{int i, j;for (i = 0; i < 17; i++){printf("%4.1f ", i*0.1 + 0.5);for (j = 0; j < 5; j++){printf("%11.9f ", inverse_sqrt(j, i*0.1 + 0.5));}printf("\r\n");}return 0;
}

快速反正弦(Arcsin)计算

原文最终选择的是多项式近似, 避免了取绝对值等复杂处理, 只是在 x = ± 1 x = \pm 1 x=±1 附近的绝对精度较差, 输出值规范化为 π \pi π,即定义 arcsin ⁡ ( x ) = π × a s i n ( x ) \arcsin(x) = \pi \times asin(x) arcsin(x)=π×asin(x). 计算式为
a s i n ( x ) = x 2 × a 0 + a 2 x 2 + a 4 x 4 + a 6 x 6 b 0 + b 2 x 2 + b 4 x 4 + b 6 x 6 asin(x) = \frac{x}{2} \times \frac{a_0 + a_2x^2 + a_4x^4 + a_6x^6}{b_0 + b_2x^2 + b_4x^4 + b_6x^6} asin(x)=2x×b0+b2x2+b4x4+b6x6a0+a2x2+a4x4+a6x6
对应的系数数值为
a 0 = 0.318309886 a 2 = − 0.5182875 a 4 = 0.222375 a 6 = − 0.016850156 b 0 = 0.5 b 2 = − 0.89745875 b 4 = 0.46138125 b 6 = − 0.058377188 a_0 = 0.318309886 \\ a_2 = -0.5182875 \\ a_4 = 0.222375 \\ a_6 = -0.016850156 \\ \\ b_0 = 0.5 \\ b_2 = -0.89745875 \\ b_4 = 0.46138125 \\ b_6 = -0.058377188 a0=0.318309886a2=0.5182875a4=0.222375a6=0.016850156b0=0.5b2=0.89745875b4=0.46138125b6=0.058377188

∣ x ∣ < 0.75 |x|<0.75 x<0.75时, 绝对误差小于 1 × 1 0 − 5 1 \times 10^{-5} 1×105, 当 ∣ x ∣ < 0.91 |x|<0.91 x<0.91时, 绝对误差小于 4.2 × 1 0 − 5 4.2 \times 10^{-5} 4.2×105, 在 x ≈ 0.997 x \approx 0.997 x0.997时, 最大误差为 0.011 0.011 0.011.

C语言实现

const float coeffa[4] = {// a0 ~ a60.318309886,-0.5182875,0.222375,-0.016850156
};const float coeffb[4] = {0.5,-0.89745875,0.46138125,-0.058377188
};const float pi = 3.14159265358979;float arcsin(float x)
{int i;float x2 = 1, a = 0, b = 0;for (i = 0; i < 4; i ++){a = a + coeffa[i] * x2;b = b + coeffb[i] * x2;x2 = x2 * x * x;}return (x * pi / 2) * (a / b);
}int main(int argc, char *const argv[])
{int i;float x, alpha, expect;for (i = 0; i < 20; i++){x = 0.01 + (i * 0.05);alpha = arcsin(x);expect= asin(x);printf("x:%4.2f  a:%9.6f e:%9.6f\r\n", x, alpha, expect);}
}

计算结果, 最右侧一列为 math.h 的 asin() 函数, 用于对比

x:0.01  a: 0.010000 e: 0.010000
x:0.06  a: 0.060036 e: 0.060036
x:0.11  a: 0.110223 e: 0.110223
x:0.16  a: 0.160691 e: 0.160691
x:0.21  a: 0.211575 e: 0.211575
x:0.26  a: 0.263022 e: 0.263022
x:0.31  a: 0.315193 e: 0.315193
x:0.36  a: 0.368268 e: 0.368268
x:0.41  a: 0.422454 e: 0.422454
x:0.46  a: 0.477996 e: 0.477995
x:0.51  a: 0.535185 e: 0.535185
x:0.56  a: 0.594386 e: 0.594386
x:0.61  a: 0.656060 e: 0.656061
x:0.66  a: 0.720815 e: 0.720819
x:0.71  a: 0.789485 e: 0.789498
x:0.76  a: 0.863278 e: 0.863313
x:0.81  a: 0.944073 e: 0.944152
x:0.86  a: 1.035139 e: 1.035270
x:0.91  a: 1.143404 e: 1.143284
x:0.96  a: 1.291451 e: 1.287002

将0.9附近细分一下

x:0.90  a: 1.119760 e: 1.119769
x:0.91  a: 1.143404 e: 1.143284
x:0.92  a: 1.168431 e: 1.168081
x:0.93  a: 1.195150 e: 1.194413
x:0.94  a: 1.224027 e: 1.222630
x:0.95  a: 1.255752 e: 1.253236
x:0.96  a: 1.291451 e: 1.287002
x:0.97  a: 1.333107 e: 1.325231
x:0.98  a: 1.384628 e: 1.370462
x:0.99  a: 1.455034 e: 1.429257

0 < x < 0.6 0 < x < 0.6 0<x<0.6区间的精度最高, 在 0.6 < x < 0.9 0.6 < x < 0.9 0.6<x<0.9区间结果数值偏小, 在 0.9 < x < 1 0.9 < x < 1 0.9<x<1区间结果数值偏大. 越接近1精度越差, 实际使用时在大于 0.97 0.97 0.97时需要做一些补偿.

参考

  • 用多项式快速计算三角函数等 https://www.olliw.eu/2014/fast-functions/

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/716331.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【 10X summary report】怎么看?详细解读笔记

报告内容 在开始正式的分析之前&#xff0c;需要查看在对齐和计数过程中生成的任何总结统计信息。下图是由Cell Ranger工具创建的10X总结报告&#xff0c;在从10X scRNA-seq实验生成计数矩阵时会生成。 The left half of the report describes sequencing and mapping statist…

卖wordpress网站模板的网站

WP模板牛 http://www.wpniu.com 上面有很多免费wordpress模板资源的网站&#xff0c;除了免费模板&#xff0c;还有付费模板。 My模板(我的模板) http://www.mymoban.com 老牌网站模板资源站&#xff0c;上面有wordpress模板、帝国CMS模板、WooCommerce模板可以直接免费下载…

Linux whois命令教程:查询域名所有者信息(附案例详解和注意事项)

Linux whois命令介绍 whois命令是一个用于查询域名所有者信息的工具。它可以直接从命令行进行查询&#xff0c;这对于没有图形用户界面的系统或者需要在shell脚本中进行查询的情况非常有用。 Linux whois命令适用的Linux版本 whois命令在大多数Linux发行版中都可以使用&…

C++之stack

1、stack简介 stack是实现的一个先进后出&#xff0c;后进先出的容器。它只有一个出口&#xff0c;只能操作最顶端元素。 2、stack库函数 &#xff08;1&#xff09;push() //向栈压入一个元素 &#xff08;2&#xff09;pop() //移除栈顶元素 &#xff08;3…

基于springboot+vue的中国陕西民俗网

博主主页&#xff1a;猫头鹰源码 博主简介&#xff1a;Java领域优质创作者、CSDN博客专家、阿里云专家博主、公司架构师、全网粉丝5万、专注Java技术领域和毕业设计项目实战&#xff0c;欢迎高校老师\讲师\同行交流合作 ​主要内容&#xff1a;毕业设计(Javaweb项目|小程序|Pyt…

在 Angular 中使用 Renderer2

Renderer2 类 Renderer2 类是 Angular 提供的一个抽象服务&#xff0c;允许在不直接操作 DOM 的情况下操纵应用程序的元素。这是推荐的方法&#xff0c;因为它使得更容易开发可以在没有 DOM 访问权限的环境中渲染的应用程序&#xff0c;比如在服务器上、在 Web Worker 中或在原…

Java如何剪切视频

背景&#xff1a;如何使用Java批量切割视频 FFmpeg 是一个强大的开源多媒体处理工具&#xff0c;被广泛应用于音视频的录制、转码、编辑等方面。它支持几乎所有主流的音视频格式&#xff0c;能够在各种操作系统平台上运行&#xff0c;包括 Windows、macOS 和 Linux。FFmpeg 提…

nginx,php-fpm

一&#xff0c;Nginx是异步非阻塞多进程&#xff0c;io多路复用 1、master进程&#xff1a;管理进程 master进程主要用来管理worker进程&#xff0c;具体包括如下4个主要功能&#xff1a; &#xff08;1&#xff09;接收来自外界的信号。 &#xff08;2&#xff09;向各worker进…

SAP PP学习笔记04 - BOM2 -通过Serial来做简单的BOM变式配置,副明细,BOM状态,BOM明细状态,项目种类,递归BOM

本章继续讲BOM。 本章讲通过Serial来做简单的BOM变式配置。还讲了BOM的相关概念&#xff1a;副明细&#xff0c;BOM状态&#xff0c;BOM明细状态&#xff0c;项目种类&#xff0c;递归BOM 等。 1&#xff0c;通过Serial&#xff08;序列号&#xff09;来做简单的 VC&#xff0…

spring自定义注解之-ElementType.METHOD方法级注解声明

自定义注解类型和常用场景 可以参考之前的文章 &#xff1a; ElementType.FIELD字段级注解声明 如果在项目中&#xff0c;多处地方都需调用到同一个方法进行逻辑处理&#xff0c;且与方法的业务逻辑无关&#xff0c;比如监控&#xff0c;日志等&#xff0c;则可用自定义的方法…

【JavaSE】面向对象——继承性

继承性 继承性的概念 所谓继承&#xff0c;就是程序猿在保持原有类特性的基础上进行扩展&#xff0c;增加新功能&#xff0c;这样的类被称为派生类或者子类&#xff0c;原有类被称为超类或者基类。 在对于继承性概念进行书写前&#xff0c;我曾查阅许多资料来保证对其表达的…

Some collections -- 2024.3

一、TensorFlow Android (dataset: Mnist) We used TensorFlow to define and train our machine learning model, which can recognize handwritten numbers, called a number classifier model in machine learning terminology. We transform the trained TensorFlow mod…

C++学习第五天(内存管理)

1、内存分布 int globalVar 1; static int staticGlobalVar 1; void Test() {static int staticVar 1;int localVar 1;int num1[10] { 1, 2, 3, 4 };char char2[] "abcd";const char* pChar3 "abcd";int* ptr1 (int*)malloc(sizeof(int) * 4);int…

2024.03.01作业

1. 基于UDP的TFTP文件传输 #include "test.h"#define SER_IP "192.168.1.104" #define SER_PORT 69 #define IP "192.168.191.128" #define PORT 9999enum mode {TFTP_READ 1,TFTP_WRITE 2,TFTP_DATA 3,TFTP_ACK 4,TFTP_ERR 5 };void get_…

高维中介数据:基于交替方向乘子法(ADMM)的高维度单模态中介模型的参数估计(入门+实操)

全文摘要 用于高维度单模态中介模型的参数估计&#xff0c;采用交替方向乘子法&#xff08;ADMM&#xff09;进行计算。该包提供了确切独立筛选&#xff08;SIS&#xff09;功能来提高中介效应的敏感性和特异性&#xff0c;并支持Lasso、弹性网络、路径Lasso和网络约束惩罚等不…

npm 镜像源切换与设置

项目背景 依赖安装中断或响应特别慢。 可以看到当前所用的镜像是 https://registry.npmjs.org 。 切换淘宝镜像之后总算能够安装下来 命令行模式 查看当前镜像源 # 查看当前镜像源 npm config get registry 可以看到默认情况下是官方默认全局镜像 https://registry.npmjs.o…

竞争加剧下,登顶后的瑞幸该做什么?

瑞幸咖啡仅用短短18个月时间从品牌创立到纳斯达克上市&#xff0c;刷新全球最快上市记录。2020年因交易造假事件被勒令退市股价暴跌80%&#xff0c;有人说这个创造了赴美IPO奇迹的“巨婴”将是下一个倒下的ofo。2023年瑞幸咖啡以逆势超速增长领跑咖啡赛道有力回应了市场的质疑&…

Vector中的begin和end函数是左闭右开的区间

vector::end() 函数的语法 vector::end(); 参数&#xff1a; none——它什么都不接受。 返回值&#xff1a; iterator– 它返回一个指向向量的 past-the-end 元素的迭代器。 实际上Vector中的begin和end函数是左闭右开的区间。 例&#xff1a; Input: vector<int>…

Java多线程实现发布和订阅

目录 简介 步骤 1: 定义消息类 步骤 2: 创建发布者 步骤 3: 创建订阅者 步骤 4: 实现发布-订阅模型 前言-与正文无关 生活远不止眼前的苦劳与奔波&#xff0c;它还充满了无数值得我们去体验和珍惜的美好事物。在这个快节奏的世界中&#xff0c;我们往往容易陷入工作的漩涡…

棋牌室计时计费管理系统的灯控器连接教程

棋牌室计时计费管理系统的灯控器连接教程 一、前言 以下教程以 佳易王棋牌室计时计费管理系统软件V18.0为例说明 软件文件下载可以点击最下方官网卡片——软件下载——试用版软件下载 如上图&#xff0c;计时计费软件在开始计时的时候&#xff0c;点击 开始计时 如果连接了…