做FPGA的人迟早会遇到CORDIC。不管是电机控制里实时算转子的反正切角度还是信号处理里生成正弦波又或者是图像旋转、数字下变频几乎都绕不开这个老算法。我第一次正经接触CORDIC是在一个永磁同步电机的FOC控制项目里要在一个小规模FPGA上同时处理电流环和角度估算DSP资源紧张得很查找表又觉得精度和面积的平衡不好把握后来才把CORDIC真正捡起来用。这个算法最吸引人的地方在于实现三角函数、反正切、极坐标转换这些看起来需要乘法器和ROM的东西本质上只需要加减法和移位这在FPGA上简直是天作之合。这篇博文想聊的就是基于FPGA的CORDIC算法实现用Verilog写从原理推导开始到定点数设计再到可综合的流水线代码最后附上仿真验证和实际工程里的坑。不管你是刚入门FPGA想找一个练手项目还是已经在项目里被各种角度计算折磨这篇文章都能给你一套直接能用的思路。1. CORDIC为什么在FPGA里这么吃香1.1 先搞清楚CORDIC到底在算什么事CORDIC全称是Coordinate Rotation Digital Computer坐标旋转数字计算方法是1959年由J.E. Volder提出的。它的核心思路非常暴力任何一个角度旋转都可以拆成一系列固定小角度的旋转组合。每次旋转只做一件事——把坐标系的x和y按照某个特定角度旋转一下而这个特定角度经过精心设计使得旋转操作只涉及移位和加法。举个例子如果要把向量(x, y)旋转角度θ传统做法是x x·cosθ - y·sinθ y y·cosθ x·sinθ这需要至少两个乘法器和一个sin/cos查找表在FPGA上要么费DSP资源要么费BRAM。CORDIC的做法是把θ拆成θ d0·α0 d1·α1 ... dn·αn其中d_i是±1表示第i次旋转的方向α_i是预先选好的角度序列。关键来了如果令α_i arctan(2^(-i))那么第i次旋转的公式就变成了x_{i1} x_i - d_i·y_i·2^(-i) y_{i1} y_i d_i·x_i·2^(-i) z_{i1} z_i - d_i·α_i乘上2^(-i)在二进制里就是右移i位硬件上直接拉线就行一个乘法器都不用。所谓“坐标旋转数字计算机”就是用数字逻辑里的移位和加减法模拟坐标旋转本质消耗的是逻辑单元LUT和寄存器而不是DSP和BRAM。1.2 为什么不用查找表或者直接调DSP我知道很多人第一反应是算个sin/cos而已我直接用ROM存一个查找表多简单或者干脆用乘法器搭FPGA里DSP资源也没那么稀缺。这话对了一半但分场景。查找表方案在小范围、低精度场景确实简单粗暴。比如只算0到360度、1度步进的正弦值360个点存进ROM查表输出逻辑简单到不行。但一旦遇到高精度或者大动态范围查找表的面积会爆炸。你想要0.001度的分辨率就得存36万个点BRAM吃得厉害。而且查找表只能解决单输入单输出的问题遇到需要实时计算atan(y/x)这种双输入函数查找表根本没法直接做。DSP方案也一样。用DSP48算乘法确实快但一个项目中同时要用FIR滤波器、PID控制器、FFT蝶形运算DSP资源往往是第一个被抢光的。CORDIC用逻辑资源换DSP资源在一些DSP吃紧的工程里是救命稻草。还有一个常被忽略的点CORDIC迭代结构天然适合流水线化。做一个N级流水线CORDIC每级只做一次移位和一次加减一级组合逻辑延迟很小时钟频率很容易拉高。吞吐率能做到每时钟周期输出一个结果延迟固定为N个周期这对信号处理链路里的实时性要求太友好了。查找表虽然可以做到单周期出结果但面积和精度之间的矛盾解决不了。所以CORDIC在FPGA领域几十年不过时不是因为它有情怀而是它在资源、速度、精度的三角权衡里找到了一个很好的平衡点。2. 动手前必须想清楚的几个设计决策2.1 旋转模式和向量模式两种基本姿势CORDIC有两种基本工作模式用途完全不一样很多人刚接触时容易混。旋转模式Rotation Mode是求三角函数用的。输入一个角度z初始坐标(x0, y0)通常设为(1/K, 0)每级迭代根据z的符号决定旋转方向不断逼近目标角度。迭代完成后x_n ≈ cos(z0)y_n ≈ sin(z0)。这里的K是CORDIC迭代带来的模长增益迭代无限次时收敛到1.64676所以初始值放1/K是为了把增益抵消掉。向量模式Vectoring Mode是求反正切和模长用的。输入一个向量(x0, y0)初始z0通常设为0每级迭代根据y的符号决定旋转方向让y不断逼近0。迭代完成后z_n ≈ arctan(y0/x0)x_n ≈ √(x0² y0²)。这两种模式的区别只在于判断旋转方向的依据旋转模式看z的符号向量模式看y的符号。代码实现上几乎一样很多模块只需要改一个判断条件就能复用。我用一个parameter来控制模式选择一套代码走天下。2.2 迭代次数、数据位宽和角度量化怎么确定这三个参数直接决定算法精度和资源消耗是CORDIC工程化最关键的一步。迭代次数N决定角度逼近的理论精度。CORDIC每次迭代的残差大约在2^(-N)量级想做到10^(-4)精度至少要14次迭代做到10^(-5)级别需要17次左右。但迭代不是越多越好超过一定次数后定点数的截断误差会盖过迭代误差再加级数也没意义。工程上我一般取16到20级配合适当的位宽扩展。数据位宽是另一个关键。角度量化和坐标数据必须独立设计。角度如果范围是±π用二进制补码表示需要2个整数位能表示到±2π剩下的都是小数位。坐标数据要考虑迭代过程中的模长增益最大绝对值可能接近1.647所以至少需要2个整数位再加一个小数位作冗余防止中间运算溢出。角度表arctan(2^(-i))的量化也要注意。每个角度值都要乘以2^FF是小数位宽取整存成定点数。这里建议用四舍五入而不是直接截断能减少大约半个LSB的误差。角度表可以用case语句硬编码也可以用逻辑自动生成我后面会给出Python生成脚本。关于位宽选择还有一个经验之谈CORDIC内部运算的数据位宽最好比输入输出位宽大3到4位用来吸收中间迭代的截断误差。比如输入输出是16位内部寄存器至少干到20位。这个余量看起来“浪费”实际对精度提升非常明显尤其是在迭代级数比较深的场景。2.3 增益补偿最容易翻车的一个环节CORDIC每次旋转的旋转矩阵并不是纯正交的它带的缩放因子是√(12^(-2i))。全部N次迭代完成后向量模长会乘上一个总增益K ∏√(12^(-2i))i从0到N-1当N足够大时K ≈ 1.64676。如果不做补偿旋转模式算出来的cos/sin会整体偏大1.64676倍向量模式算出来的模长也会偏大这个倍数。这是个非常典型的新手错误仿真波形看起来形状是对的幅值却不对查了半天发现忘了处理增益。补偿方法有两种一是迭代结束后把结果乘以1/K≈0.60725需要一次乘法二是更常用的小技巧直接把初始x0设为1/K的定点值让每次迭代都从被缩小的初始值出发最后得到的结果天然就是正确幅值。第二种方法在旋转模式下特别好用因为初始x0是固定常量不增加任何额外逻辑。在向量模式下初始x0、y0是外部输入的数据没法预置但模长的补偿可以在后续处理里乘回来或者干脆放在系统的下一个处理级里。有一点要注意如果用预置初始值的方法x0 1/K必须量化到位宽范围内。比如18位定点2整数15小数1/K ≈ 0.60725量化后是19900左右远小于32767不会溢出。但如果你把数据位宽压缩得很小比如8位定点0.60725量化后大约是78也能放下不过精度会损失不少所以位宽别太小。3. Verilog实现从接口到流水线逐级拆解3.1 模块接口定义与整体架构我先给一个典型的旋转模式CORDIC模块接口输入是相位角输出是cos和sin。这个模块可以直接当数控振荡器NCO用也可以作为底层IP集成到更大的系统里。module cordic_rot #( parameter PHASE_WIDTH 18, // 输入角度位宽1符号2整数15小数 parameter DATA_WIDTH 20, // 内部数据位宽 parameter ITER_NUM 16 // 迭代级数 )( input wire clk, input wire rst_n, input wire en, // 输入有效 input wire signed [PHASE_WIDTH-1:0] phase_in, // 弧度定点数 output reg signed [DATA_WIDTH-1:0] cos_out, output reg signed [DATA_WIDTH-1:0] sin_out, output reg valid // 输出有效 );内部采用16级流水线结构每级有两组寄存器一组是旋转后的x、y坐标一组是剩余角度z。第i级的旋转方向由z_i的符号决定。输入使能信号en在第一级打入然后跟随流水线逐级传递同时用移位寄存器产生valid信号这样就知道输出数据哪一拍是有效的。这里有个很重要的设计点valid信号的延迟必须和数据处理路径完全对齐。我见过不少工程因为valid延迟没对齐导致后续模块采到无效数据整个链路输出乱跳。最简单的方法就是写一个和流水线深度相同的移位寄存器en信号进valid出。3.2 角度常量表的生成与量化角度表是CORDIC的“灵魂”每个角度值对应一次迭代的旋转步长。16级迭代需要16个角度arctan(2^0)、arctan(2^-1)、arctan(2^-2)……一直到arctan(2^-15)。这些角度是固定常量可以在Verilog里用localparam硬编码也可以用case语句生成。我习惯用Python先把角度值算好再生成Verilog代码片段避免手算出错。量化位数和数据位宽挂钩比如PHASE_WIDTH18位去掉1位符号和2位整数小数位F15角度值就乘以2^15四舍五入取整。import math F 15 # 小数位宽 ITER 16 for i in range(ITER): angle math.atan(2 ** (-i)) q int(round(angle * (2 ** F))) print(f15d{q}, // i{i}, atan{angle:.6f})生成的定点角度表可以直接放到一个case语句里或者放进一个只读数组。注意这些值是正数不需要符号位。每次迭代根据z的符号决定加上还是减去对应角度。在Verilog里case语句的形式如下function [PHASE_WIDTH-2:0] get_atan_table; input [3:0] idx; begin case (idx) 4d0: get_atan_table 15d25736; // atan(1) 4d1: get_atan_table 15d15193; // atan(0.5) 4d2: get_atan_table 15d8026; // atan(0.25) // ... 其他迭代级 default: get_atan_table 15d0; endcase end endfunction具体量化值会因为F取15略有不同实际工程里用Python重新生成一份最稳妥。3.3 流水线迭代级代码详解下面是核心的单级流水线逻辑。以第i级为例输入是x_in、y_in、z_in输出是x_out、y_out、z_out。旋转方向由z_in的最高位决定z_in为正符号位为0逆时针旋转z_in为负符号位为1顺时针旋转。wire z_neg z_in[PHASE_WIDTH-1]; // 符号位1表示负数 wire [DATA_WIDTH-1:0] x_shift x_in i; // 算术右移保留符号 wire [DATA_WIDTH-1:0] y_shift y_in i; always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_out 0; y_out 0; z_out 0; end else begin if (z_neg) begin // 顺时针旋转x加y的移位项y减x的移位项z加角度 x_out x_in y_shift; y_out y_in - x_shift; z_out z_in angle_i; end else begin // 逆时针旋转x减y的移位项y加x的移位项z减角度 x_out x_in - y_shift; y_out y_in x_shift; z_out z_in - angle_i; end end end有两点必须强调。第一右移必须用算术右移因为x_in和y_in是有符号数如果用逻辑右移负数的高位会补0整个数据会变成巨大的正数结果彻底错误。第二每级移位量i是不同的第0级右移0位第1级右移1位第2级右移2位。流水线展开写时可以直接用generate-for循环参数化避免手工复制16遍代码。使用generate-for的完整流水线结构大致是这样genvar i; generate for (i 0; i ITER_NUM; i i 1) begin: cordic_stage cordic_stage #( .DATA_WIDTH(DATA_WIDTH), .PHASE_WIDTH(PHASE_WIDTH), .SHIFT(i), .ANGLE(get_atan_table(i)) ) u_stage ( .clk(clk), .rst_n(rst_n), .x_in(i 0 ? x0 : stage_x[i-1]), .y_in(i 0 ? y0 : stage_y[i-1]), .z_in(i 0 ? phase_abs : stage_z[i-1]), .x_out(stage_x[i]), .y_out(stage_y[i]), .z_out(stage_z[i]) ); end endgenerate实际上get_atan_table用function在generate里调用会有综合限制更稳妥的做法是把角度表定义成常量数组然后按循环变量i索引localparam [PHASE_WIDTH-2:0] atan_table [0:ITER_NUM-1] { 18d25736, 18d15193, // 必须和实际位宽对齐 // ... }; // generate块内 wire signed [PHASE_WIDTH-1:0] angle_i {{1{1b0}}, atan_table[i]};总之流水线展开的关键是每级用不同的移位量和不同的角度值其他逻辑完全一样。综合工具会复制16份相同的逻辑资源消耗大约是串行结构的16倍但换来的是每个时钟周期都能出结果。3.4 象限映射预处理与初始值注入之前说CORDIC迭代的收敛范围大约在±99.88度超过这个范围算法不收敛。但实际输入角度可能覆盖±π甚至更大所以必须在进入迭代之前做象限映射。常用的做法是把所有角度折叠到第一象限0到π/2用符号标志记录折叠方式迭代完成后根据标志恢复cos/sin的符号。象限映射逻辑不难但要仔细。我写过一版输入phase_in先取绝对值判断它落在哪个象限然后映射到第一象限的角度。映射规则是第二象限π-x映射成x第三象限x-π映射成x第四象限-x映射成x。cos和sin的最终符号根据原始象限决定。初始值注入方面旋转模式把x0预置为1/K的定点值y0为0。1/K的量化和数据位宽相关比如DATA_WIDTH20取2整数17小数1/K×2^177962320位有符号数最大是524287没问题。如果采用不预置而事后补偿的方案初始x0就直接赋1最后再加一个乘法器对结果乘0.60725。两种方案各有优劣预置法省逻辑但限制初始值事后补偿法灵活但要消耗乘法器。4. 仿真验证与实测误差到底能压到多低4.1 Testbench怎么写才能真正发现问题很多人写CORDIC的testbench就是给几个特定角度跑一遍看波形对不对。这样做不是不行但对精度评估来说信息量太少了。我建议至少做两类测试定点对比测试和全范围扫描测试。定点对比测试是拿CORDIC输出和Matlab/Python计算的浮点参考值对比。把输入角度和输出数据都导出成文本然后在Python里统一做误差统计。可以算最大绝对误差、均方根误差还可以看误差的分布是随机噪声还是系统性偏差。系统性偏差通常来自增益补偿不精确或者角度表量化误差随机噪声主要来自迭代截断。区分这两类误差对定位问题很有帮助。全范围扫描测试是把输入角度从-π到π按一定步长扫描把所有输出结果记录下来。重点检查两个东西一是角度在边界0、±π/2、±π附近时输出是否正确这些地方最容易出现符号翻转错误二是检查整个范围内误差是否均匀如果某个区间误差突然变大很可能象限映射逻辑在那个边界出问题了。Testbench里还可以加一个随机角度激励用$urandom生成随机相位持续跑几千个周期统计valid输出和采样时刻是否稳定。这种随机测试特别适合找时序对齐的问题。4.2 实测误差分析迭代次数和位宽对精度的影响我实际做过一组对比实验固定输入输出位宽18位分别用12级、16级、20级迭代扫描全范围角度测得最大绝对误差大致如下迭代级数最大绝对误差cos最大绝对误差sin备注12约3.2e-3约3.5e-3迭代误差主导16约2.1e-4约2.3e-4逼近理论精度20约1.9e-4约2.0e-4受限于数据位宽不再明显改善从数据能看到两个规律迭代从12级加到16级误差降低了一个数量级但从16级加到20级误差几乎没有变化。原因就是内部数据位宽20位截断误差已经盖过了迭代误差。想进一步提升精度光加迭代级数没用得同时加数据位宽。另一个观察是误差在角度接近±π/2时略微偏大。这是因为cos在0附近比较平坦对角度误差不敏感sin在π/2附近导数接近0同样不敏感但CORDIC在迭代过程中z的残差总会留一点表现在输出上就是靠近某些角度时会出现小的“台阶”。这个现象在要求高精度信号合成的场景需要注意如果完全不能接受可以后级加一个小的补偿修正或者提高定点数位宽。4.3 资源占用与性能实测在Cyclone V上一块中等规模的FPGA上16级流水线、数据位宽20位的CORDIC模块实测资源占用大约是不到700个Logic ElementLE基本不用DSP和BRAM。最高时钟频率可以跑到200MHz以上每周期输出一个结果。这个资源占用和性能指标对绝大多数中低端FPGA都是很轻松的负担。作为对比如果是串行迭代结构资源占用会降到200个LE以内但计算一次需要16个时钟周期吞吐率只有流水线的1/16。所以架构选择完全看项目需求如果正弦波发生器后面接的DAC采样率只有几十kHz串行结构绰绰有余如果要用在高速解调链路里数据率是几十兆甚至上百兆流水线结构才是正解。5. 常见问题速查与避坑经验5.1 精度上不去先查这几个地方遇到CORDIC精度不如预期我总结了一套排查顺序。第一查数据位宽。内部位宽如果和输出位宽一样没有余量那么每个中间级的截断误差一路累积到最后精度通常会比理论值差很多。把内部位宽扩到输出位宽加3到4位精度立刻改善。第二查角度表量化。角度表如果直接截断而不是四舍五入每个角度值会损失约0.5个LSB累积起来是系统性的负偏。把角度表全部换成四舍五入量化能明显减少系统性误差。第三查增益补偿。如果算出来的cos/sin波形形状对但幅值整体偏大约1.647倍那就是增益没补偿。用预置初始值1/K的方法最省事但注意1/K量化后要检查是否超出位宽范围。第四查象限映射。如果误差只在特定角度区间突然恶化大概率是象限映射的边界情况处理错了。比如角度正好等于π/2时符号位的判断和绝对值处理要特别小心。5.2 时序收敛和资源优化经验CORDIC流水线每一级的组合逻辑很浅就是一次移位和一次加减所以时序收敛一般都挺顺利。但如果你把16级全部写成一个大的组合逻辑块中间不加寄存器组合逻辑链会非常长时序必挂。记得每一级之间要有一级寄存器采用真正意义上的流水线结构。如果资源很紧张可以考虑半流水设计把16级迭代拆成两组8级的反馈结构每个时钟周期算半级4个时钟周期算完一级迭代总面积大约只有全流水的1/4吞吐率降低4倍。这个方案在很多中低速率应用里性价比很高。另一个优化技巧是角度表的实现。16个宽度18位的小常量用LUT就能放下不需要用BRAM。但如果你把它写成一个ROM有些综合工具会为了优化把它放进BRAM反而浪费了宝贵的存储资源。可以在综合属性里把角度表声明为luts来避免这种情况。5.3 流水线深度带来的对齐问题使用流水线CORDIC时数据从输入到输出需要固定N个时钟周期的延迟。如果你的系统里还有其他并行路径必须确保这些路径的延迟和CORDIC一致否则后续模块会采到时间上没对齐的数据。这个问题的隐蔽性很高因为单看CORDIC本身功能正确但整个链路一配合就出错。解决办法是给每个路径加上匹配延迟。用SRL16或者移位寄存器做延迟对齐是最常见的方式。工程上还要注意延迟匹配要按周期精确对齐哪怕差一个周期输出的星座图或者波形都会有明显的相位跳变。6. 项目实战两个可以直接套用的应用场景6.1 用CORDIC做高精度数控振荡器NCO数控振荡器NCO在通信系统里几乎是标配。传统NCO用相位累加器加查找表实现查找表的深度直接决定频率分辨率和杂散性能。用CORDIC替代查找表后一个很大的优势是可以把ROM省掉而且输出正弦/余弦是逐点实时计算出来的不存在查找表的相位截断杂散。我在一个数字解调项目里用过这个结构。相位累加器输出32位相位取高18位送给CORDICCORDIC输出cos和sin给后级混频器。实测无杂散动态范围比同参数的查找表方案高了不少关键是查表方案为了压低杂散需要存两万个点的ROM而CORDIC方案几乎不占BRAM。把省下来的BRAM留给后面的滤波器整体链路性能提升明显。6.2 用向量模式做电机控制里的角度与速度估算永磁同步电机的无感FOC控制里经常需要根据αβ轴的反电动势估算转子角度和速度。这个本质上就是算atan(eβ/eα)。用CORDIC向量模式把eα和eβ作为输入x0、y0迭代完成后z输出就是转子角度x输出近似幅值一举两得。实际工程里还有个细节CORDIC向量模式在x0为负时需要特殊处理否则求出来的角度范围会不对。我的做法是先判断x0的符号如果为负就把输入向量旋转180度得到一个正x0的等效向量迭代完成后再把角度加回π。这样整个角度范围可以达到-π到π满足FOC控制的全部需求。类似的处理在很多通信和伺服应用里都会遇到可以作为一个通用技巧记住。写到这里关于CORDIC算法的原理、定点化设计、Verilog实现和工程实践经验基本都讲到了。我个人实际做下来的体会是CORDIC这种“用移位和加减法做数学运算”的思路在现代FPGA资源日益丰富的今天依然有着不可替代的价值。它不追求算法本身的花哨胜在结构规整、时序友好、资源开销可控。如果你手头正好有角度计算、三角函数或者坐标变换的需求不妨按这篇文章的思路搭一个出来试试跑一遍仿真看看误差曲线再放到板子上测一测实际效果整个流程走通之后你会对FPGA里的数值运算有一种完全不同的掌控感。
FPGA上CORDIC算法实现：移位与加法搞定三角函数与反正切
做FPGA的人迟早会遇到CORDIC。不管是电机控制里实时算转子的反正切角度还是信号处理里生成正弦波又或者是图像旋转、数字下变频几乎都绕不开这个老算法。我第一次正经接触CORDIC是在一个永磁同步电机的FOC控制项目里要在一个小规模FPGA上同时处理电流环和角度估算DSP资源紧张得很查找表又觉得精度和面积的平衡不好把握后来才把CORDIC真正捡起来用。这个算法最吸引人的地方在于实现三角函数、反正切、极坐标转换这些看起来需要乘法器和ROM的东西本质上只需要加减法和移位这在FPGA上简直是天作之合。这篇博文想聊的就是基于FPGA的CORDIC算法实现用Verilog写从原理推导开始到定点数设计再到可综合的流水线代码最后附上仿真验证和实际工程里的坑。不管你是刚入门FPGA想找一个练手项目还是已经在项目里被各种角度计算折磨这篇文章都能给你一套直接能用的思路。1. CORDIC为什么在FPGA里这么吃香1.1 先搞清楚CORDIC到底在算什么事CORDIC全称是Coordinate Rotation Digital Computer坐标旋转数字计算方法是1959年由J.E. Volder提出的。它的核心思路非常暴力任何一个角度旋转都可以拆成一系列固定小角度的旋转组合。每次旋转只做一件事——把坐标系的x和y按照某个特定角度旋转一下而这个特定角度经过精心设计使得旋转操作只涉及移位和加法。举个例子如果要把向量(x, y)旋转角度θ传统做法是x x·cosθ - y·sinθ y y·cosθ x·sinθ这需要至少两个乘法器和一个sin/cos查找表在FPGA上要么费DSP资源要么费BRAM。CORDIC的做法是把θ拆成θ d0·α0 d1·α1 ... dn·αn其中d_i是±1表示第i次旋转的方向α_i是预先选好的角度序列。关键来了如果令α_i arctan(2^(-i))那么第i次旋转的公式就变成了x_{i1} x_i - d_i·y_i·2^(-i) y_{i1} y_i d_i·x_i·2^(-i) z_{i1} z_i - d_i·α_i乘上2^(-i)在二进制里就是右移i位硬件上直接拉线就行一个乘法器都不用。所谓“坐标旋转数字计算机”就是用数字逻辑里的移位和加减法模拟坐标旋转本质消耗的是逻辑单元LUT和寄存器而不是DSP和BRAM。1.2 为什么不用查找表或者直接调DSP我知道很多人第一反应是算个sin/cos而已我直接用ROM存一个查找表多简单或者干脆用乘法器搭FPGA里DSP资源也没那么稀缺。这话对了一半但分场景。查找表方案在小范围、低精度场景确实简单粗暴。比如只算0到360度、1度步进的正弦值360个点存进ROM查表输出逻辑简单到不行。但一旦遇到高精度或者大动态范围查找表的面积会爆炸。你想要0.001度的分辨率就得存36万个点BRAM吃得厉害。而且查找表只能解决单输入单输出的问题遇到需要实时计算atan(y/x)这种双输入函数查找表根本没法直接做。DSP方案也一样。用DSP48算乘法确实快但一个项目中同时要用FIR滤波器、PID控制器、FFT蝶形运算DSP资源往往是第一个被抢光的。CORDIC用逻辑资源换DSP资源在一些DSP吃紧的工程里是救命稻草。还有一个常被忽略的点CORDIC迭代结构天然适合流水线化。做一个N级流水线CORDIC每级只做一次移位和一次加减一级组合逻辑延迟很小时钟频率很容易拉高。吞吐率能做到每时钟周期输出一个结果延迟固定为N个周期这对信号处理链路里的实时性要求太友好了。查找表虽然可以做到单周期出结果但面积和精度之间的矛盾解决不了。所以CORDIC在FPGA领域几十年不过时不是因为它有情怀而是它在资源、速度、精度的三角权衡里找到了一个很好的平衡点。2. 动手前必须想清楚的几个设计决策2.1 旋转模式和向量模式两种基本姿势CORDIC有两种基本工作模式用途完全不一样很多人刚接触时容易混。旋转模式Rotation Mode是求三角函数用的。输入一个角度z初始坐标(x0, y0)通常设为(1/K, 0)每级迭代根据z的符号决定旋转方向不断逼近目标角度。迭代完成后x_n ≈ cos(z0)y_n ≈ sin(z0)。这里的K是CORDIC迭代带来的模长增益迭代无限次时收敛到1.64676所以初始值放1/K是为了把增益抵消掉。向量模式Vectoring Mode是求反正切和模长用的。输入一个向量(x0, y0)初始z0通常设为0每级迭代根据y的符号决定旋转方向让y不断逼近0。迭代完成后z_n ≈ arctan(y0/x0)x_n ≈ √(x0² y0²)。这两种模式的区别只在于判断旋转方向的依据旋转模式看z的符号向量模式看y的符号。代码实现上几乎一样很多模块只需要改一个判断条件就能复用。我用一个parameter来控制模式选择一套代码走天下。2.2 迭代次数、数据位宽和角度量化怎么确定这三个参数直接决定算法精度和资源消耗是CORDIC工程化最关键的一步。迭代次数N决定角度逼近的理论精度。CORDIC每次迭代的残差大约在2^(-N)量级想做到10^(-4)精度至少要14次迭代做到10^(-5)级别需要17次左右。但迭代不是越多越好超过一定次数后定点数的截断误差会盖过迭代误差再加级数也没意义。工程上我一般取16到20级配合适当的位宽扩展。数据位宽是另一个关键。角度量化和坐标数据必须独立设计。角度如果范围是±π用二进制补码表示需要2个整数位能表示到±2π剩下的都是小数位。坐标数据要考虑迭代过程中的模长增益最大绝对值可能接近1.647所以至少需要2个整数位再加一个小数位作冗余防止中间运算溢出。角度表arctan(2^(-i))的量化也要注意。每个角度值都要乘以2^FF是小数位宽取整存成定点数。这里建议用四舍五入而不是直接截断能减少大约半个LSB的误差。角度表可以用case语句硬编码也可以用逻辑自动生成我后面会给出Python生成脚本。关于位宽选择还有一个经验之谈CORDIC内部运算的数据位宽最好比输入输出位宽大3到4位用来吸收中间迭代的截断误差。比如输入输出是16位内部寄存器至少干到20位。这个余量看起来“浪费”实际对精度提升非常明显尤其是在迭代级数比较深的场景。2.3 增益补偿最容易翻车的一个环节CORDIC每次旋转的旋转矩阵并不是纯正交的它带的缩放因子是√(12^(-2i))。全部N次迭代完成后向量模长会乘上一个总增益K ∏√(12^(-2i))i从0到N-1当N足够大时K ≈ 1.64676。如果不做补偿旋转模式算出来的cos/sin会整体偏大1.64676倍向量模式算出来的模长也会偏大这个倍数。这是个非常典型的新手错误仿真波形看起来形状是对的幅值却不对查了半天发现忘了处理增益。补偿方法有两种一是迭代结束后把结果乘以1/K≈0.60725需要一次乘法二是更常用的小技巧直接把初始x0设为1/K的定点值让每次迭代都从被缩小的初始值出发最后得到的结果天然就是正确幅值。第二种方法在旋转模式下特别好用因为初始x0是固定常量不增加任何额外逻辑。在向量模式下初始x0、y0是外部输入的数据没法预置但模长的补偿可以在后续处理里乘回来或者干脆放在系统的下一个处理级里。有一点要注意如果用预置初始值的方法x0 1/K必须量化到位宽范围内。比如18位定点2整数15小数1/K ≈ 0.60725量化后是19900左右远小于32767不会溢出。但如果你把数据位宽压缩得很小比如8位定点0.60725量化后大约是78也能放下不过精度会损失不少所以位宽别太小。3. Verilog实现从接口到流水线逐级拆解3.1 模块接口定义与整体架构我先给一个典型的旋转模式CORDIC模块接口输入是相位角输出是cos和sin。这个模块可以直接当数控振荡器NCO用也可以作为底层IP集成到更大的系统里。module cordic_rot #( parameter PHASE_WIDTH 18, // 输入角度位宽1符号2整数15小数 parameter DATA_WIDTH 20, // 内部数据位宽 parameter ITER_NUM 16 // 迭代级数 )( input wire clk, input wire rst_n, input wire en, // 输入有效 input wire signed [PHASE_WIDTH-1:0] phase_in, // 弧度定点数 output reg signed [DATA_WIDTH-1:0] cos_out, output reg signed [DATA_WIDTH-1:0] sin_out, output reg valid // 输出有效 );内部采用16级流水线结构每级有两组寄存器一组是旋转后的x、y坐标一组是剩余角度z。第i级的旋转方向由z_i的符号决定。输入使能信号en在第一级打入然后跟随流水线逐级传递同时用移位寄存器产生valid信号这样就知道输出数据哪一拍是有效的。这里有个很重要的设计点valid信号的延迟必须和数据处理路径完全对齐。我见过不少工程因为valid延迟没对齐导致后续模块采到无效数据整个链路输出乱跳。最简单的方法就是写一个和流水线深度相同的移位寄存器en信号进valid出。3.2 角度常量表的生成与量化角度表是CORDIC的“灵魂”每个角度值对应一次迭代的旋转步长。16级迭代需要16个角度arctan(2^0)、arctan(2^-1)、arctan(2^-2)……一直到arctan(2^-15)。这些角度是固定常量可以在Verilog里用localparam硬编码也可以用case语句生成。我习惯用Python先把角度值算好再生成Verilog代码片段避免手算出错。量化位数和数据位宽挂钩比如PHASE_WIDTH18位去掉1位符号和2位整数小数位F15角度值就乘以2^15四舍五入取整。import math F 15 # 小数位宽 ITER 16 for i in range(ITER): angle math.atan(2 ** (-i)) q int(round(angle * (2 ** F))) print(f15d{q}, // i{i}, atan{angle:.6f})生成的定点角度表可以直接放到一个case语句里或者放进一个只读数组。注意这些值是正数不需要符号位。每次迭代根据z的符号决定加上还是减去对应角度。在Verilog里case语句的形式如下function [PHASE_WIDTH-2:0] get_atan_table; input [3:0] idx; begin case (idx) 4d0: get_atan_table 15d25736; // atan(1) 4d1: get_atan_table 15d15193; // atan(0.5) 4d2: get_atan_table 15d8026; // atan(0.25) // ... 其他迭代级 default: get_atan_table 15d0; endcase end endfunction具体量化值会因为F取15略有不同实际工程里用Python重新生成一份最稳妥。3.3 流水线迭代级代码详解下面是核心的单级流水线逻辑。以第i级为例输入是x_in、y_in、z_in输出是x_out、y_out、z_out。旋转方向由z_in的最高位决定z_in为正符号位为0逆时针旋转z_in为负符号位为1顺时针旋转。wire z_neg z_in[PHASE_WIDTH-1]; // 符号位1表示负数 wire [DATA_WIDTH-1:0] x_shift x_in i; // 算术右移保留符号 wire [DATA_WIDTH-1:0] y_shift y_in i; always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_out 0; y_out 0; z_out 0; end else begin if (z_neg) begin // 顺时针旋转x加y的移位项y减x的移位项z加角度 x_out x_in y_shift; y_out y_in - x_shift; z_out z_in angle_i; end else begin // 逆时针旋转x减y的移位项y加x的移位项z减角度 x_out x_in - y_shift; y_out y_in x_shift; z_out z_in - angle_i; end end end有两点必须强调。第一右移必须用算术右移因为x_in和y_in是有符号数如果用逻辑右移负数的高位会补0整个数据会变成巨大的正数结果彻底错误。第二每级移位量i是不同的第0级右移0位第1级右移1位第2级右移2位。流水线展开写时可以直接用generate-for循环参数化避免手工复制16遍代码。使用generate-for的完整流水线结构大致是这样genvar i; generate for (i 0; i ITER_NUM; i i 1) begin: cordic_stage cordic_stage #( .DATA_WIDTH(DATA_WIDTH), .PHASE_WIDTH(PHASE_WIDTH), .SHIFT(i), .ANGLE(get_atan_table(i)) ) u_stage ( .clk(clk), .rst_n(rst_n), .x_in(i 0 ? x0 : stage_x[i-1]), .y_in(i 0 ? y0 : stage_y[i-1]), .z_in(i 0 ? phase_abs : stage_z[i-1]), .x_out(stage_x[i]), .y_out(stage_y[i]), .z_out(stage_z[i]) ); end endgenerate实际上get_atan_table用function在generate里调用会有综合限制更稳妥的做法是把角度表定义成常量数组然后按循环变量i索引localparam [PHASE_WIDTH-2:0] atan_table [0:ITER_NUM-1] { 18d25736, 18d15193, // 必须和实际位宽对齐 // ... }; // generate块内 wire signed [PHASE_WIDTH-1:0] angle_i {{1{1b0}}, atan_table[i]};总之流水线展开的关键是每级用不同的移位量和不同的角度值其他逻辑完全一样。综合工具会复制16份相同的逻辑资源消耗大约是串行结构的16倍但换来的是每个时钟周期都能出结果。3.4 象限映射预处理与初始值注入之前说CORDIC迭代的收敛范围大约在±99.88度超过这个范围算法不收敛。但实际输入角度可能覆盖±π甚至更大所以必须在进入迭代之前做象限映射。常用的做法是把所有角度折叠到第一象限0到π/2用符号标志记录折叠方式迭代完成后根据标志恢复cos/sin的符号。象限映射逻辑不难但要仔细。我写过一版输入phase_in先取绝对值判断它落在哪个象限然后映射到第一象限的角度。映射规则是第二象限π-x映射成x第三象限x-π映射成x第四象限-x映射成x。cos和sin的最终符号根据原始象限决定。初始值注入方面旋转模式把x0预置为1/K的定点值y0为0。1/K的量化和数据位宽相关比如DATA_WIDTH20取2整数17小数1/K×2^177962320位有符号数最大是524287没问题。如果采用不预置而事后补偿的方案初始x0就直接赋1最后再加一个乘法器对结果乘0.60725。两种方案各有优劣预置法省逻辑但限制初始值事后补偿法灵活但要消耗乘法器。4. 仿真验证与实测误差到底能压到多低4.1 Testbench怎么写才能真正发现问题很多人写CORDIC的testbench就是给几个特定角度跑一遍看波形对不对。这样做不是不行但对精度评估来说信息量太少了。我建议至少做两类测试定点对比测试和全范围扫描测试。定点对比测试是拿CORDIC输出和Matlab/Python计算的浮点参考值对比。把输入角度和输出数据都导出成文本然后在Python里统一做误差统计。可以算最大绝对误差、均方根误差还可以看误差的分布是随机噪声还是系统性偏差。系统性偏差通常来自增益补偿不精确或者角度表量化误差随机噪声主要来自迭代截断。区分这两类误差对定位问题很有帮助。全范围扫描测试是把输入角度从-π到π按一定步长扫描把所有输出结果记录下来。重点检查两个东西一是角度在边界0、±π/2、±π附近时输出是否正确这些地方最容易出现符号翻转错误二是检查整个范围内误差是否均匀如果某个区间误差突然变大很可能象限映射逻辑在那个边界出问题了。Testbench里还可以加一个随机角度激励用$urandom生成随机相位持续跑几千个周期统计valid输出和采样时刻是否稳定。这种随机测试特别适合找时序对齐的问题。4.2 实测误差分析迭代次数和位宽对精度的影响我实际做过一组对比实验固定输入输出位宽18位分别用12级、16级、20级迭代扫描全范围角度测得最大绝对误差大致如下迭代级数最大绝对误差cos最大绝对误差sin备注12约3.2e-3约3.5e-3迭代误差主导16约2.1e-4约2.3e-4逼近理论精度20约1.9e-4约2.0e-4受限于数据位宽不再明显改善从数据能看到两个规律迭代从12级加到16级误差降低了一个数量级但从16级加到20级误差几乎没有变化。原因就是内部数据位宽20位截断误差已经盖过了迭代误差。想进一步提升精度光加迭代级数没用得同时加数据位宽。另一个观察是误差在角度接近±π/2时略微偏大。这是因为cos在0附近比较平坦对角度误差不敏感sin在π/2附近导数接近0同样不敏感但CORDIC在迭代过程中z的残差总会留一点表现在输出上就是靠近某些角度时会出现小的“台阶”。这个现象在要求高精度信号合成的场景需要注意如果完全不能接受可以后级加一个小的补偿修正或者提高定点数位宽。4.3 资源占用与性能实测在Cyclone V上一块中等规模的FPGA上16级流水线、数据位宽20位的CORDIC模块实测资源占用大约是不到700个Logic ElementLE基本不用DSP和BRAM。最高时钟频率可以跑到200MHz以上每周期输出一个结果。这个资源占用和性能指标对绝大多数中低端FPGA都是很轻松的负担。作为对比如果是串行迭代结构资源占用会降到200个LE以内但计算一次需要16个时钟周期吞吐率只有流水线的1/16。所以架构选择完全看项目需求如果正弦波发生器后面接的DAC采样率只有几十kHz串行结构绰绰有余如果要用在高速解调链路里数据率是几十兆甚至上百兆流水线结构才是正解。5. 常见问题速查与避坑经验5.1 精度上不去先查这几个地方遇到CORDIC精度不如预期我总结了一套排查顺序。第一查数据位宽。内部位宽如果和输出位宽一样没有余量那么每个中间级的截断误差一路累积到最后精度通常会比理论值差很多。把内部位宽扩到输出位宽加3到4位精度立刻改善。第二查角度表量化。角度表如果直接截断而不是四舍五入每个角度值会损失约0.5个LSB累积起来是系统性的负偏。把角度表全部换成四舍五入量化能明显减少系统性误差。第三查增益补偿。如果算出来的cos/sin波形形状对但幅值整体偏大约1.647倍那就是增益没补偿。用预置初始值1/K的方法最省事但注意1/K量化后要检查是否超出位宽范围。第四查象限映射。如果误差只在特定角度区间突然恶化大概率是象限映射的边界情况处理错了。比如角度正好等于π/2时符号位的判断和绝对值处理要特别小心。5.2 时序收敛和资源优化经验CORDIC流水线每一级的组合逻辑很浅就是一次移位和一次加减所以时序收敛一般都挺顺利。但如果你把16级全部写成一个大的组合逻辑块中间不加寄存器组合逻辑链会非常长时序必挂。记得每一级之间要有一级寄存器采用真正意义上的流水线结构。如果资源很紧张可以考虑半流水设计把16级迭代拆成两组8级的反馈结构每个时钟周期算半级4个时钟周期算完一级迭代总面积大约只有全流水的1/4吞吐率降低4倍。这个方案在很多中低速率应用里性价比很高。另一个优化技巧是角度表的实现。16个宽度18位的小常量用LUT就能放下不需要用BRAM。但如果你把它写成一个ROM有些综合工具会为了优化把它放进BRAM反而浪费了宝贵的存储资源。可以在综合属性里把角度表声明为luts来避免这种情况。5.3 流水线深度带来的对齐问题使用流水线CORDIC时数据从输入到输出需要固定N个时钟周期的延迟。如果你的系统里还有其他并行路径必须确保这些路径的延迟和CORDIC一致否则后续模块会采到时间上没对齐的数据。这个问题的隐蔽性很高因为单看CORDIC本身功能正确但整个链路一配合就出错。解决办法是给每个路径加上匹配延迟。用SRL16或者移位寄存器做延迟对齐是最常见的方式。工程上还要注意延迟匹配要按周期精确对齐哪怕差一个周期输出的星座图或者波形都会有明显的相位跳变。6. 项目实战两个可以直接套用的应用场景6.1 用CORDIC做高精度数控振荡器NCO数控振荡器NCO在通信系统里几乎是标配。传统NCO用相位累加器加查找表实现查找表的深度直接决定频率分辨率和杂散性能。用CORDIC替代查找表后一个很大的优势是可以把ROM省掉而且输出正弦/余弦是逐点实时计算出来的不存在查找表的相位截断杂散。我在一个数字解调项目里用过这个结构。相位累加器输出32位相位取高18位送给CORDICCORDIC输出cos和sin给后级混频器。实测无杂散动态范围比同参数的查找表方案高了不少关键是查表方案为了压低杂散需要存两万个点的ROM而CORDIC方案几乎不占BRAM。把省下来的BRAM留给后面的滤波器整体链路性能提升明显。6.2 用向量模式做电机控制里的角度与速度估算永磁同步电机的无感FOC控制里经常需要根据αβ轴的反电动势估算转子角度和速度。这个本质上就是算atan(eβ/eα)。用CORDIC向量模式把eα和eβ作为输入x0、y0迭代完成后z输出就是转子角度x输出近似幅值一举两得。实际工程里还有个细节CORDIC向量模式在x0为负时需要特殊处理否则求出来的角度范围会不对。我的做法是先判断x0的符号如果为负就把输入向量旋转180度得到一个正x0的等效向量迭代完成后再把角度加回π。这样整个角度范围可以达到-π到π满足FOC控制的全部需求。类似的处理在很多通信和伺服应用里都会遇到可以作为一个通用技巧记住。写到这里关于CORDIC算法的原理、定点化设计、Verilog实现和工程实践经验基本都讲到了。我个人实际做下来的体会是CORDIC这种“用移位和加减法做数学运算”的思路在现代FPGA资源日益丰富的今天依然有着不可替代的价值。它不追求算法本身的花哨胜在结构规整、时序友好、资源开销可控。如果你手头正好有角度计算、三角函数或者坐标变换的需求不妨按这篇文章的思路搭一个出来试试跑一遍仿真看看误差曲线再放到板子上测一测实际效果整个流程走通之后你会对FPGA里的数值运算有一种完全不同的掌控感。