FPGA圈子待久了你会发现CORDIC这个名字隔三差五就会跳出来——做数字解调要算相位做电机控制要算sin/cos做频谱分析要从FFT结果里抠幅度和角度甚至做图像旋转都躲不开它。我最早接触CORDIC是在一个QPSK解调项目里那时候图省事直接拖了IP核后来换平台发现IP核版本不一致、位宽还要重新配折腾了几次之后下定决心自己用Verilog写一个。这篇文章就是把整个思路摊开聊CORDIC到底在算什么、参数怎么定、流水线怎么写、误差哪里来、上板有哪些坑。适合已经写过一些FPGA模块、想摆脱只会调IP核状态的工程师也适合准备面试时被问CORDIC原理的同学。1. CORDIC的数学底子如何把三角函数变成移位和加法1.1 FPGA上为什么不能直接算sin在MCU里一句sin(x)背后可能是查表加插值也可能是多项式逼近反正CPU帮你扛了。FPGA没有这种运行时函数库想在一个时钟周期内拿到sin值最朴素的想法是搞一张大查找表输入角度做地址输出直接是结果。表越大精度越高但16位精度就需要64K深度的ROM资源占用很可观。另一个思路是用泰勒展开或者多项式逼近问题是FPGA里做乘法的代价不低阶数一高组合逻辑路径变长时序立马崩给你看。CORDICCoordinate Rotation Digital Computer坐标旋转数字计算机换了个思路不直接算角度对应的三角函数值而是让一个已知向量通过一系列微旋转逐步逼近目标角度。每次旋转的角度预先定好让这个角度的正切值恰好是2的负整数次幂于是旋转操作就退化成了移位和加法。FPGA里移位是布线就干完的事加法器是基本功所以整个算法天然适合硬件实现。1.2 伪旋转先丢掉模长再干掉乘法常规的二维旋转公式是x x·cosθ - y·sinθ y x·sinθ y·cosθ把cosθ提出来可以改写成x cosθ·(x - y·tanθ) y cosθ·(x·tanθ y)如果我们暂时不管向量长度被缩放了多少只关心角度转得对不对那么括号里的部分就够用了x x - y·tanθ y x·tanθ y这一步在CORDIC里叫伪旋转pseudo-rotation。问题是tanθ还是不好算。CORDIC的关键跳跃在于把目标角度θ拆成一串特殊角θ_i arctan(2^(-i))的和。因为tanθ_i 2^(-i)伪旋转公式变成x x - y · 2^(-i) y x · 2^(-i) y右移i位就是乘2^(-i)。在硬件上右移不占任何逻辑资源于是旋转一个任意角度被拆成了若干次右移加减法。1.3 角度累加器越转越接近目标每次旋转到底是加还是减取决于当前角度累加器z的符号。初始时z等于你输入的目标角度如果z还大于0说明转得还不够这一级就往正方向转一个θ_i如果z已经小于0说明转过了就往回转一个θ_iz_{i1} z_i - d_i·θ_i其中d_i 1当z_i≥0否则d_i -1这个过程本质上是一个二分式的角度逼近。经过N级迭代后z会收敛到接近0而x和y则收敛到目标角度对应的cos和sin值。角度分辨率大约等于2^(-N)弧度迭代到第16级时分辨率已经到0.00087°量级对付绝大多数工控和通信场景绰绰有余。1.4 被放大的模长1.64676这个常数的来历伪旋转虽然角度对了但每转一次向量长度都会被放大。观察伪旋转矩阵[[1, -2^(-i)], [2^(-i), 1]]它的行列式是1 2^(-2i)所以每一步模长会乘sqrt(1 2^(-2i))。把所有级连起来Gain ∏ sqrt(1 2^(-2i))i从0到∞收敛到约1.6467602也就是说经过足够多级迭代x和y都被放大了1.64676倍。要在最后拿到真实的cos和sin值就得乘上0.6072529这个倒数。我习惯把补偿放在输入端初始向量x直接给0.607而不是1.0这样整个流水线内部的数值始终不会超过1后面截位、符号处理都省心。如果你想在调试时直接看到理论上的cos/sin中间值也可以把补偿放在输出端代价是内部数据位宽要额外留2bit因为中间最大会涨到1.647。2. 动手前的关键决策迭代级数、角度格式与象限折叠2.1 迭代级数选16级背后是精度和资源的平衡迭代级数和精度之间有一个非常直观的关系N级迭代之后剩余的角度误差不超过最后一个未参与迭代的角度近似等于2^(-N)弧度。换算成角度2^(-16)弧度大约是0.00087°这个精度对应大约14.8位有效二进制位。而如果输入角度用的是16位有符号数整个角度域的量化精度是π/32768 ≈ 0.0055°也就是说16级迭代的角度分辨率已经远小于输入角度的量化间隔了。我把不同迭代级数和精度的关系整理了一下迭代级数理论角度误差上限16位角度下是否够用资源消耗8 级约0.22°不够误差明显低12 级约0.014°勉强接近量化极限中16 级约0.00087°足够低于1个LSB中高20 级约0.000055°过剩精度受位宽限制高还有一个工程上的细节当迭代超过数据位宽对应的精度极限后角度查找表里的值已经小于1个LSB继续迭代只是在空转消耗加法器和寄存器却带不来可感知的精度提升。所以16位数据宽度配14到16级迭代是甜点区我最终选了16级。2.2 角度格式1.0代表π不是代表2π这是新手最容易翻车的地方。CORDIC的标准数学描述里角度范围是[-π, π)但是FPGA里只有二进制数没有π这个单位。大家习惯用定点数表示角度但不统一的是满量程到底对应多少。我见过不少实现把16位有符号数的满量程对应成2π也就是0x8000代表-π、0x7FFF代表接近π。这样每个LSB对应的角度是2π/65536看起来相位累加器用起来方便但和CORDIC迭代角度表配对时特别容易出错。我的做法是让1.0对应π输入范围[-1, 1)映射到[-π, π)0.5对应π/2-0.5对应-π/2。这样做的直接好处是角度累加器的最高两位天然就是象限编号做象限折叠时零额外开销。对应的arctan查找表的每一项也要按角度/π来存。比如第0级角度是arctan(1) π/4在16位定点下存的就是0.25 × 32768 8192。这一点务必和整个系统的角度格式保持一致否则仿真到一半会发现输出完全对不上。2.3 象限折叠让CORDIC只工作在有效收敛域内标准CORDIC的收敛域大约在±99.7°以内直接输入120°或者-120°迭代过程会发散输出完全不对。好在三角函数有对称性我们可以把任意角度先折回第一或第四象限算完后再根据原始象限恢复符号。具体做法是看输入角度的高两位输入角度范围angle_in[15:14]折叠后角度输出sin输出cos0 ~ π/200angle_insin_corecos_coreπ/2 ~ π010x8000 - angle_insin_core-cos_core-π ~ -π/210angle_in 0x8000-sin_core-cos_core-π/2 ~ 011-angle_in-sin_corecos_core这里0x8000在16位有符号数里是-1按补码加法折回时恰好起到加π的效果。折叠之后的mapped_angle一定落在第一象限范围内送到CORDIC核心去迭代就万无一失。象限信息本身要打一拍寄存下来和最后的sin_core、cos_core对齐再去做符号恢复。3. 流水线Verilog核心实现级间信号流与可复用代码3.1 为什么用全流水线而不是状态机复用CORDIC的每一次迭代都依赖上一次的结果天然是一个串行链。用状态机加共享ALU的方式实现面积省但每个角度要算16拍吞吐率只有流水线的1/16。在FPGA这种并行资源多的平台我更倾向于全流水16级每一级都放一组独立的加法器和移位逻辑输入连续送入虽然延迟固定16拍但每个时钟周期都能吐一个结果。内部数据位宽我选择在输入16位基础上加宽2bit用18位S3.15格式做中间运算。原因在1.4里已经提到伪旋转会让模长最大增长到1.64676倍如果中间也用16位迭代到第10级左右就会溢出饱和输出波形会出现一段平头。加宽2bit后中间数据最大到1.647也能装下最后输出前再截回16位。3.2 参数化模块代码下面是一个旋转模式rotation mode的CORDIC核心代码输入是16位归一化角度输出是Q1.15格式的sin和cos。为了保持文章篇幅可控我略去了最终的符号恢复寄存器组但核心级联逻辑是完整的。module cordic_rotation #( parameter WIDTH 16, parameter ITER 16 )( input wire clk, input wire rst_n, input wire valid_in, input wire signed [WIDTH-1:0] angle_in, // 归一化角度, 1.0代表pi output reg signed [WIDTH-1:0] sin_out, output reg signed [WIDTH-1:0] cos_out, output reg valid_out ); localparam IW WIDTH 2; // 内部18位数据宽度 localparam signed [WIDTH-1:0] INIT_X 16d19899; // 0.60725的Q15表示 // 角度查找表: atan(2^(-i)) / pi * 32768 function [WIDTH-1:0] get_angle; input integer idx; begin case (idx) 0: get_angle 16d8192; 1: get_angle 16d4836; 2: get_angle 16d2555; 3: get_angle 16d1297; 4: get_angle 16d651; 5: get_angle 16d326; 6: get_angle 16d163; 7: get_angle 16d81; 8: get_angle 16d41; 9: get_angle 16d20; 10: get_angle 16d10; 11: get_angle 16d5; 12: get_angle 16d3; 13: get_angle 16d1; default: get_angle 16d0; endcase end endfunction // 流水线寄存器组 reg signed [IW-1:0] x_reg [0:ITER]; reg signed [IW-1:0] y_reg [0:ITER]; reg signed [WIDTH-1:0] z_reg [0:ITER]; reg [1:0] quad_reg [0:ITER]; reg valid_reg [0:ITER]; wire [1:0] quadrant angle_in[WIDTH-1:WIDTH-2]; // 象限折叠 reg signed [WIDTH-1:0] mapped_angle; always (*) begin case (quadrant) 2b00: mapped_angle angle_in; 2b01: mapped_angle 16sh8000 - angle_in; 2b10: mapped_angle angle_in 16sh8000; default: mapped_angle -angle_in; endcase end // 第0级: 输入初始化 always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_reg[0] {2b00, INIT_X}; y_reg[0] {IW{1b0}}; z_reg[0] {2b00, mapped_angle}; quad_reg[0] 2b00; valid_reg[0] 1b0; end else begin x_reg[0] {2b00, INIT_X}; y_reg[0] {IW{1b0}}; z_reg[0] {{2{1b0}}, mapped_angle}; quad_reg[0] quadrant; valid_reg[0] valid_in; end end // 16级流水线 genvar i; generate for (i 0; i ITER; i i 1) begin : pipe_stage wire signed [IW-1:0] x_shift x_reg[i] i; wire signed [IW-1:0] y_shift y_reg[i] i; wire signed [WIDTH-1:0] angle_step get_angle(i); always (posedge clk or negedge rst_n) begin if (!rst_n) begin x_reg[i1] {IW{1b0}}; y_reg[i1] {IW{1b0}}; z_reg[i1] {WIDTH{1b0}}; quad_reg[i1] 2b00; valid_reg[i1] 1b0; end else if (valid_reg[i]) begin if (z_reg[i] 0) begin x_reg[i1] x_reg[i] - y_shift; y_reg[i1] y_reg[i] x_shift; z_reg[i1] z_reg[i] - angle_step; end else begin x_reg[i1] x_reg[i] y_shift; y_reg[i1] y_reg[i] - x_shift; z_reg[i1] z_reg[i] angle_step; end quad_reg[i1] quad_reg[i]; valid_reg[i1] 1b1; end else begin valid_reg[i1] 1b0; end end end endgenerate // 输出级: 补偿增益并截回16位 wire signed [IW-1:0] x_core x_reg[ITER]; wire signed [IW-1:0] y_core y_reg[ITER]; always (posedge clk or negedge rst_n) begin if (!rst_n) begin sin_out {WIDTH{1b0}}; cos_out {WIDTH{1b0}}; valid_out 1b0; end else begin // 0.60725 * 32768 ≈ 19900, 右移15位去掉增益 sin_out (y_core * 16d19900) 15; cos_out (x_core * 16d19900) 15; valid_out valid_reg[ITER]; end end endmodule这段代码里有几个地方要特别说明。x_reg、y_reg用 i做算术右移这样对有符号数是带符号扩展的右移如果用不同仿真器对signed reg的行为可能不一致上板综合也可能出现你意想不到的结果所以这里统一用算术右移。角度查找表虽然写成了function但综合时所有调用点的输入都是常数i所以会被综合成16组常量不会真的生成一块ROM。象限信息quad_reg和valid_reg要跟着流水线一级一级打拍因为最后恢复符号时必须用和x_core、y_core同一拍的象限。如果图省事直接用最开始的quadrant输出会错位一拍表现是波形整体对不上而且越到后面越难排查。3.3 有效标志与复位流水线对齐和亚稳态处理整个流水线在复位释放后的第一个有效输入进入经过16级迭代、1级补偿输出valid_out会滞后valid_in正好17拍。valid_reg每级打一拍就是为了生成一个和输出数据严格对齐的有效标志下游模块直接采这个标志就行不用关心具体延迟。复位方面FPGA里常见的坑是异步复位信号本身有亚稳态风险。我的习惯是异步复位、同步释放也就是rst_n先经过两级触发器同步再作为整个模块的复位信号。CORDIC内部没有状态机不需要特殊的上电初始化序列复位之后所有寄存器清零流水线自然进入空转状态有valid才往下走。4. 仿真与误差实测偏差到底从哪里来4.1 测试平台怎么搭CORDIC的testbench比一般模块要稍微讲究一点因为输出有固定延迟还要和理想值做差。我的做法是写一个for循环从-16384到16383遍历输入角度相当于从-π到π扫一遍每个点都和$cos、$sin的理想值比较统计最大绝对误差和均方根误差。注意参考模型的输入角度也要归一化也就是角度值除以π这样才和模块内部格式一致。reg signed [15:0] angle_tb; reg valid_tb; wire signed [15:0] sin_tb, cos_tb; reg signed [15:0] sin_ref, cos_ref; real angle_real; real err_abs; initial begin for (integer i -16384; i 16384; i i 1) begin angle_tb i; valid_tb 1; angle_real $itor(i) / 16384.0 * 3.1415926535; sin_ref $rtoi($sin(angle_real) * 32768.0); cos_ref $rtoi($cos(angle_real) * 32768.0); repeat(20) (posedge clk); end valid_tb 0; repeat(40) (posedge clk); $finish; end扫全范围的价值在于能直观看到误差在哪些点最大。实测下来我那个16位实现的最大绝对误差基本在1到2个LSB之间也就是0.005°到0.011°左右。这个精度水平对于大部分解调、测角、电机控制应用都够用了。4.2 三个误差源分别是什么把误差拆开看主要来自三方面。量化误差是最底的输入角度只有16位本来就只有有限个可表示的角度这个误差无法消除。迭代截断误差来自只迭代16级剩余角度理论上还有2^(-16)弧度左右但前面算过这个值小于角度量化LSB淹没在量化误差里。最后是补偿截位误差乘0.60725后右移15位低位被丢掉的1bit就是±1LSB的量级。三个误差叠加最大误差落在1到2个LSB是非常正常的。误差源数值量级是否可以消除输入角度量化误差1 LSB ≈ 0.0055°不能受位宽限制迭代截断误差~0.00087°增加级数可降低增益补偿截位误差±1 LSB增加输出位宽可降低4.3 仿真环境与常见报错如果你用的是Intel版ModelSim或Questa注意免费的Starter Edition对代码行数有限制CORDIC这种参数化模块展开后行数可能超标。遇到failure to obtain a Verilog simulation license这类报错优先查环境变量LM_LICENSE_FILE有没有指向正确的license文件或者干脆换开源的Icarus Verilog加GTKWave做功能仿真CORDIC不依赖厂商库开源工具完全跑得动。仿真时还有一个隐蔽的问题有符号数和无符号数的混用。比如象限折叠里用了16sh8000如果你写的是16h8000在一些工具里会被当成无符号数处理相减的结果就完全不对了。所以所有角度、数据通路的信号声明里signed关键字一定不能漏。5. 上板前后踩过的坑从仿真通过到实际可用的距离5.1 增益补偿位置不对中间级溢出我第一次写这个模块时输入x直接给了1.0想着最后乘0.60725拉回来就行。仿真单点看没错但扫全范围时发现0°和90°附近的输出不对仔细看波形x_reg在第9级就因为溢出饱和了。根源就是伪旋转的模长放大效应中间值最大到1.64716位Q1.15表示不下。后来把输入x改成0.60725并且把内部位宽加宽到18位这个问题才彻底解决。5.2 象限符号恢复时的边界尖刺象限恢复的表格看起来简单但有一个边界情况很坑当cos_m正好等于-32768时取负运算-cos_m在16位补码里会溢出回绕成-32768本身导致输出突然跳变。虽然这个点只在cos-1、sin0附近出现但如果你做频谱分析时恰好有信号落在那里会观察到一个奇怪的毛刺。我在实际项目里把输出限幅在[-32767, 32767]解决代价是损失了1个LSB的动态范围换来了输出单调性。5.3 用组合逻辑直通16级时序收敛不了有人为了省寄存器把16级迭代全部用组合逻辑直通。仿真毫无问题但综合布局布线之后发现Fmax只有几十MHz根本跑不到需求。CORDIC的每一级是一组加法器和移位器16级级联的组合逻辑路径相当长。如果对资源特别敏感可以做迭代式结构状态机控制一个ALU循环16次面积能缩到几分之一但吞吐率也会降到1/16。我的经验是在现在的FPGA资源条件下全流水线通常是最划算的——16级流水只多占16组寄存器却能保证每个周期出一个结果。5.4 什么时候别自己写和IP核的取舍自己写CORDIC最大的优势是透明位宽、迭代级数、延迟、中间数据格式全在自己手里想怎么改怎么改不依赖厂商IP核的版本和授权。IP核的优势在于经过了充分验证支持浮点、高精度、资源优化等高级选项项目周期紧的时候直接拖一个最稳。我的建议是如果你需要20位以上的精度、或者要处理IEEE754浮点输入、或者时间真的来不及用IP核如果你想彻底搞清楚算法细节、要裁剪成任意位宽、要在不同厂商平台间无缝迁移自己写一定值得。Xilinx那个老应用笔记XAPP523是迭代式结构原理讲得很清楚代码风格偏老建议当原理参考不要直接往新项目里搬。往后面扩展的话把旋转模式里d_i的符号源从z的符号换成y的符号同一个模块就变成了向量模式vector mode可以直接算arctan(y/x)和sqrt(x²y²)。我后来做QPSK解调时就是靠这一个CORDIC模块两种模式既算了载波相位又算了幅值。最后说一句实操心得CORDIC写完之后不管仿真多漂亮上板前先在实机上扫一个完整周期的点把误差曲线拉出来看一眼这比任何理论分析都更能让你放心。