基于DSP的三角函数快速计算

算 该方法简单 精度高 速度快且不需要额外的存储空间 任一角度 θ 的正弦和余弦组成复数 cosθ + i sin θ 设等 差角度数列的公差为 ∆ ∆ 的正弦和余弦组成复数 cos∆ + i sin ∆ 二者相乘的结果为 cos (θ + ∆ ) + i sin (θ + ∆ ) 因而可通过复数相乘得到角度 θ + ∆ 的正余弦 迭代多次就 可以得到整个序列的正余弦 该算法的缺点是由于舍入误差的存在 迭代次数越多 运算结果的累积误差越大 具体的误差大小可以由公式 (1)进 行有效的估计
x + x1 x + x1 sin x 2 ( 2 − x1 ) + sin x1 x 2 − 2 2 2
x 2 − x1
x
x + x1 = sin 2 ⋅ 1 − cos 2
2
− x1 2
从上面的结果可以看出
整个周期的最大误差出现在
由此定义造表间距 D = xi+1 − xi
有时我们会采用插值的方法 ( 常采用线性插值 ) 来得到介 于关键点之间点的函数值 4.2 算法分析 查表法主要有表的大小 精度分析和快速索引 3 个问题 需要考虑 由于存储空间的限制 查找表不可能无限大 因此查找 法适合计算定义域有限的函数或者周期性函数 可只在一个 周期内构造查找表 表的大小是造表区间与造表间距的比 值 即 L = R / D 从下面的分析可看出表的大小对精度影响 带插值的查表法与不带插值的查表法的精度分析结果 不同 对于不带插值的查表法 由于介于关键点之间点的函 数值必须用两个关键点中的某个点的函数值来近似 因此在 函数斜率最大的点误差最大 对于 sin x 来说 由于它的斜 率最大为 cos (0 ) = 1 因此误差上限是 ∆ x
将公式 (5)绝对值号中的式子求导并令之为 0 sin x 2 − sin x1 cos x − =0 x 2 − x1 由上式 可得 x + x2 x = 1 2 所以每个插值段中最大误差出现在中点处
E ( x ) = sin x −
可得 (6)
其中最大误差为
4 查表法
4.1 算法介绍 查表法适用于各种计算起来很耗时的函数 其想法是将 某些点的函数值构造成表以便需要时查找 一般只在函数定义域的部分区间内造表 称这个选取的 部分区间为造表区间 ( 记为 R ), 用关键点 x 0 , x1 , L , x n −1 , x n 将 造表区间划分成若干份 x 0 , x 1 , L , x n −1 , x n 的 函 数 值 分 查找表中记录的只是点 常采取的划分方法是等距划
将对 2π 的求模转换为对 2 的求模 而对 2 的求模可以转换 为移位运算 从而减少了运算量 这种方法对于定点数特别 13
有效 (3)第 3 种方法是浮点运算中常用的方法 该方法的思想 是将浮点的整数部分与小数部分分离 将整数部分定点化后 移位 然后将结果浮点化后与小数部分相加 4.3 DSP 实现与结果分析 若存储器的访问速度与处理器内核速度差异很大 查表 法相比其它方法没有优势 但是对于具有片内存储器的 DSP 由于访问片内存储器的速度仅需要几个周期 查表法的优越 性很明显 如果对于正余弦在 (0 , 2π ) 间构造表 并且样点间距 I 为
(中科院计算技术研究所系统结构室 北京 100080) 摘 要 分析了常用三角函数(主要是正余弦)的各种近似计算方法 包括迭代法 级数法 查表法以及 CORDIC 算法 给出了常用算法的 综合各种算法优点的实现比一般的库函数快 3~5 倍而且相对精度很高 正余弦 查表法 CORDIC 算法 DSP
误差特性 误差范围以及时空效率 依据现代 DSP 的流水 并行( SIMD)和片内存储器等特点 对各种算法进行了优化调整 提高了它们 的并行性 关键词 三角函数
则
y 的误差为
∑
n
i =1
δ y (x ) ∆xi 或 ∆y ≤ δxi
∑ max
i =1
n
δ y (x ) ∆xi δxi
(1)
Fra Baidu bibliotek
时间效率可以用周期数或程序执行时间来衡量 而空间 效率主要是指算法所消耗的存储空间以及寄存器个数
2 迭代法
迭代法适合于构成等差数列的一系列角度的正余弦计 12
作者简介 马士超(1978 ) 男 博士生 主研方向 安全芯片设计 SOC 以及嵌入式系统 王贞松 教授 博导 收稿日期 2004-09-30 E-mail scma@ict.ac.cn
2
x + x1 最大的点 也就是曲率最大 斜率最小 的点 sin 2 2 最大值为 ∆x E max = 1 − cos 2 图 2 分别是表的长度为 32 与 64 的带插值查表法计算正 弦函数的误差曲线 (y 轴为计算误差 x 轴为正弦函数的因变 量 ) 从图中可看出表的长度对于精度的影响
3.2 DSP 实现 另外由于该算法的实现过程中不需要跳转指令 因此可 以很好地利用 DSP 的 SIMD 特性同时得到多个数据的运算结 果 算法中的多项式可以采用秦九韶算法来减少计算的次数 和提高计算的并行度 举例来说 对于多项式 (2) y = ax 4 + bx 3 + cx 2 + dx + e 以秦九韶的方式可以表示为
(5)
(3)
计算式 (2)一般需要 20 条指令 若将 x 的幂次方作为中间结 果保留并调整运算顺序最少也需要 11 条指令 而式 (3)仅需 要 8 条指令即可完成计算 为了避免 DSP 中产生指令相关 可以将 (3)式进一步调整为 (4) y = ( ax + b ) x 3 + ( cx + d ) x + e 来增加算法的并行度 该调整方式对于高阶多项式尤为有效
cos x 具有类
(a) (b) 图 2 表长为 32 与 64 的带插值查表法正弦误差曲线
似的结果 图 1 分别是用表长为 32 与 64 的不带插值查表法 计算正弦函数的误差曲线 (y 轴为计算误差 x 轴为正弦函数 的因变量 ) 从图中可以看出表的长度与精度的关系
0.1 0.05 0 -.005 -0.1 0 0.1 0.05 0 -.005 -0.1 0
y = ((( ax + b ) x + c ) x + d ) x + e
sin x 在 x 点的近似值为 sin x2 ( x − x1 ) + sin x1 (x2 − x ) sin x = x2 − x1
误差可以表示为
E ( x ) = sin x − sin x 2 ( x − x 1 ) + sin x 1 ( x 2 − x ) x 2 − x1
第 31 卷 Vol.31
第 22 期 22
文章编号
计 算 机 工 程 Computer Engineering
1000 3428(2005)22 0012 03 文献标识码 A
2005 年 11 月 November 2005
中图分类号 TP301.6
博士论文
基于 ＤＳＰ 的三角函数快速计算
马士超 王贞松
3 级数法
3.1 算法介绍与分析 在不具备片内存储器的处理器上 一种通用的并且精度 相对可以很高的方法是多项式 (级数 ) 逼近 对于正弦函数 文献 [2] 中有
sin( x) = ∑ (− 1)
n =1 ∞ n
1 性能衡量标准
一般可以从数值精度 时间效率和空间效率等角度来衡 量各种算法的优劣 误差是数值计算精度的范畴 主要有绝 对误差 (Absolute error) 相对误差 (Relative error) 截断误差 舍入误差等 其中绝对误差定义为
Rapid Computation of Trigonometric Function on DSP
MA Shichao, WANG Zhensong
(Architecture Lab, Institute of Computing Technology, CAS, Beijing 100080) Abstract This paper analyses several kinds of computation method of Trigonometric function, including iterative algorithm, power serials algorithm, lookup-table algorithm and CORDIC algorithm, and gives out their error character, error volume, time efficiency and memory efficiency. Based on the pipelining, parallel and on-chip memory of DSP, the paper optimizes each algorithm to improve their parallelism. Combined all the virtue of every algorithm, the implementation runs faster than most of the commercial software library function by 3~5 times with adequate precision. Key words Trigonometric function; Cosine and sine; Lookup-table algorithm; CORDIC algorithm; DSP
error abs = f actrual − f approx
x 2 n +1 (2n + 1)!
由于计算时间以及寄存器数目的限制 只能展开级数的 有限项 若将级数展到第 N 项 那么它的理论误差就是
对于舍入误差目前一般只能就四则运算和函数运算给出 其 误 差 界 的 估 计 [2] 若 y = y( x1 , x2 ,L, xn ) xi 的误差为
(2 n + 1)!
从误差公式可以看出当 x 很小时 只需展开较少的几项 就可以得到很高的精度 而当 x 很大时 需要 N 很大才能 将 x 化归 保证级数收敛 因此一般利用正弦函数的周期性 到 ( − π , π ) 的范围后再进行计算
x
2 n +1
∆x i (i = 1,2, L , n )
∆y ≈
比较图 1 与图 2 可以发现 不带插值的查表法与带线性 插值的查表法的误差有两点不同
(1) 误差特性不同 不带插值的查表法最大误差出现在斜率最 大的点 而带插值的查表法最大误差在出现曲率最大的点 (2) 误差范围不同 在表大小相同的情况下 带插值的查表法相 对于不带插值的查表法具有更好的误差结果
如何快速地索引是查表法性能的关键 正余弦函数周期 为 2π 我们需要用求模的方法将整个定义域上的点映射到 (0 , 2π ) 之间 下面给出实现中采用的几种方法 (1) 若 模 为 A 那 么 任 意 数 x 对 A 求 模 可 以 通 过 式 y = x − floor ( x / A) * A 完 成 可 以 事 先 计 算 出 1/ A 记 为
常用三角函数 ( 尤其是正余弦 ) 在各种信号处理系统中有 着广泛的应用 且一般有实时性要求 因此有必要考虑这些 函数的快速计算 常用三角函数的近似计算方法主要有迭代 法 级数法 查表法以及 CORDIC 法 这些算法的提出和应 用有着悠久的历史 但是近年来集成电路与计算机体系结构 的飞速发展 使得各种算法具有了与以往不同的特性与结论 需要依据体系结构做适当的优化调整 本文的目的是探讨在目前的半导体工艺技术和体系结构 条件下 如何计算三角函数 ( 以正余弦为例 ) 使得它们满足特 定的精度 时间与存储要求 虽然本文的实现是针对可编程 的通用 DSP 但是对算法精度 时间效率和空间利用率的分 析与改进具有通用性
A
2
4
6
2
4
6
(a) (b) 图 1 表长为 32 与 64 的不带插值查表法正弦函数误差曲线
从而将除法运算转化为乘法运算 (2) 利用公式 x = y * π ⇒ x mod (2 * π ) = ( y mod 2 ) * π
y = x − floor( x * A) * A
对于带插值的情况 仅考虑线性插值 插值公式是 f (x 2 ) − f (x1 ) f (x ) = f (x1 ) + (x − x1 ) x 2 − x1