1. 为什么要在FPGA里用CORDIC算三角函数很多人第一次接触在FPGA上算sin和cos第一反应是去找IP核或者干脆用查找表。查找表确实简单但精度和资源是一对矛盾——你要16位精度ROM就得占一大块Block RAM你要省资源精度就惨不忍睹。而CORDIC这个算法有意思的地方在于它压根不用乘法器只靠移位和加减就能把三角函数算出来。这在FPGA里简直是天作之合因为移位在硬件里就是连线加法器是FPGA最不缺的东西。CORDIC全称是Coordinate Rotation Digital Computer坐标旋转数字计算机。名字听着唬人核心思想其实很朴素把一个角度不断分解成一系列越来越小的固定角度每次判断当前角度是正还是负决定往哪个方向转。转着转着角度就逼近目标了同时坐标也跟着转到了目标位置。旋转模式下如果你把初始向量放在x轴上x1, y0那么旋转到目标角度θ之后新的x坐标就是cos(θ)y坐标就是sin(θ)。这就是旋转模式算三角函数的本质。我选择EGo1这块板卡来做验证原因很实际。EGo1是Xilinx Artix-7系列的入门级开发板资源不算富裕但足够跑CORDIC板载的时钟、按键、LED和数码管刚好能构成一个完整的验证闭环。你不需要额外的示波器或者逻辑分析仪用板载资源就能看到结果对不对。对于学习CORDIC来说这种自闭环的验证方式比仿真波形直观得多。这篇文章适合谁看如果你已经写过Verilog或者VHDL知道什么是时序逻辑和组合逻辑但还没真正把CORDIC跑通过那这篇就是写给你的。我会从迭代公式的推导讲起把每一级流水线为什么这么设计说清楚然后给出完整的Verilog代码最后讲上板验证时怎么用数码管把结果展示出来。整个过程不需要任何IP核纯手写RTL。提示CORDIC的迭代次数和精度是直接挂钩的。N次迭代大约能得到N位的精度但角度覆盖范围有限后面会详细说这个问题怎么处理。2. CORDIC旋转模式的迭代公式到底怎么来的2.1 从一次旋转说起假设平面上有一个点(x, y)我们要把它绕原点旋转角度θ。根据旋转矩阵x x·cos(θ) - y·sin(θ) y x·sin(θ) y·cos(θ)这个公式没问题但问题在于cos(θ)和sin(θ)本身就是我们要算的东西这就成了循环依赖。CORDIC的巧妙之处在于它把旋转角θ拆成一堆预先定好的小角度每个小角度都有一个特点它的正切值是2的负整数次幂。具体来说第i次旋转的角度是atan(2^(-i))。为什么选这个角度因为tan(atan(2^(-i))) 2^(-i)而乘以2^(-i)在硬件里就是右移i位。这样一来旋转矩阵里的cos和sin就被消掉了只剩下移位和加法。把cos(θ)从旋转矩阵里提出来x cos(θ)·(x - y·tan(θ)) y cos(θ)·(y x·tan(θ))当tan(θ) ±2^(-i)时括号里的运算就变成了x cos(θ)·(x - y·d·2^(-i)) y cos(θ)·(y x·d·2^(-i))其中d是方向因子1表示逆时针转-1表示顺时针转。2.2 那个cos(θ)乘积去哪了每次迭代都有一个cos(atan(2^(-i)))的乘积因子。N次迭代下来总的乘积是K ∏ cos(atan(2^(-i))) (i从0到N-1)这个K是一个常数当N足够大时K约等于0.607252935。也就是说如果初始向量是(1, 0)经过N次迭代后得到的向量长度是K而不是1。所以如果你想要真正的cos和sin值初始x应该设为1/K也就是约1.646760258。在硬件实现里通常直接把初始x设为1/K的定点数表示。比如16位定点数1位符号1位整数14位小数1/K ≈ 1.64676对应的定点值就是1.64676 × 2^14 ≈ 26980十六进制是0x6964。2.3 角度累加器的作用每次迭代还需要一个角度累加器z初始值是目标角度θ。每次迭代根据z的符号决定旋转方向如果z 0说明当前角度还不够需要逆时针转d 1如果z 0说明转多了需要顺时针转d -1然后z更新为z - d·atan(2^(-i))。迭代足够多次后z趋近于0此时(x, y)就是(cos(θ), sin(θ))。这里有个细节atan(2^(-i))的值需要预先算好存起来。对于16位精度通常迭代16次需要16个角度常数。这些常数用定点数表示存在一个case语句或者ROM里。迭代次数iatan(2^(-i))角度值定点表示Q2.14045.000°0x2000126.565°0x12E4214.036°0x09FB37.125°0x051143.576°0x028B51.789°0x014660.895°0x00A370.448°0x005180.224°0x002990.112°0x0014100.056°0x000A110.028°0x0005120.014°0x0003130.007°0x0001140.003°0x0001150.002°0x0000注意角度范围问题CORDIC旋转模式能覆盖的角度范围是-99.9°到99.9°左右因为所有atan(2^(-i))的和约等于99.9°。如果你要算的角度超过这个范围需要先做象限预处理。比如算120°的sin可以先算60°的sin然后根据象限关系转换。这个后面在代码里会处理。3. Verilog实现从单周期迭代到全流水线3.1 迭代式实现与流水线式的取舍最直观的写法是用一个状态机每个时钟周期做一次迭代16次迭代就是16个周期出结果。这种写法资源省但吞吐率低。另一种写法是全流水线16级流水线排开每个时钟周期都能吐出一个新结果代价是寄存器用量大。在EGo1这种入门板上两种写法都能跑。我建议先用迭代式把功能调通因为迭代式调试起来简单你可以把中间每一级的x、y、z都接到ILA上观察。等确认算法没问题了再改成流水线提升吞吐率。迭代式的核心代码结构大概是这样module cordic_iterative ( input wire clk, input wire rst_n, input wire start, input wire [15:0] angle_in, // Q2.14格式范围-180°到180° output reg [15:0] cos_out, output reg [15:0] sin_out, output reg done ); // 角度常数表Q2.14格式 function [15:0] atan_table; input [3:0] idx; case (idx) 4d0: atan_table 16h2000; // 45.000° 4d1: atan_table 16h12E4; // 26.565° 4d2: atan_table 16h09FB; // 14.036° 4d3: atan_table 16h0511; // 7.125° 4d4: atan_table 16h028B; // 3.576° 4d5: atan_table 16h0146; // 1.789° 4d6: atan_table 16h00A3; // 0.895° 4d7: atan_table 16h0051; // 0.448° 4d8: atan_table 16h0029; // 0.224° 4d9: atan_table 16h0014; // 0.112° 4d10: atan_table 16h000A; // 0.056° 4d11: atan_table 16h0005; // 0.028° 4d12: atan_table 16h0003; // 0.014° 4d13: atan_table 16h0001; // 0.007° 4d14: atan_table 16h0001; // 0.003° 4d15: atan_table 16h0000; // 0.002° default: atan_table 16h0000; endcase endfunction reg [4:0] iter_cnt; reg [15:0] x, y, z; reg busy; // 象限预处理 reg [15:0] angle_pre; reg sign_flip; always (posedge clk or negedge rst_n) begin if (!rst_n) begin iter_cnt 5d0; busy 1b0; done 1b0; x 16d0; y 16d0; z 16d0; end else begin done 1b0; if (start !busy) begin busy 1b1; iter_cnt 5d0; // 初始值x 1/K ≈ 1.64676Q2.14格式 x 16h6964; y 16d0; // 角度预处理如果角度绝对值大于90°做象限转换 if (angle_in[15]) begin // 负角度取反加一得到绝对值 z ~angle_in 1b1; sign_flip 1b1; end else begin z angle_in; sign_flip 1b0; end end else if (busy) begin if (iter_cnt 5d15) begin busy 1b0; done 1b1; // 根据符号决定输出 if (sign_flip) begin cos_out x; sin_out ~y 1b1; // sin(-θ) -sin(θ) end else begin cos_out x; sin_out y; end end else begin iter_cnt iter_cnt 1b1; if (z[15] 1b0) begin // z 0逆时针转 x x - (y iter_cnt); y y (x iter_cnt); z z - atan_table(iter_cnt[3:0]); end else begin // z 0顺时针转 x x (y iter_cnt); y y - (x iter_cnt); z z atan_table(iter_cnt[3:0]); end end end end end endmodule这段代码有几个地方需要特别注意。第一是算术右移对于有符号数才能保持符号位。第二x和y的更新必须用同一时刻的旧值不能一个用新值一个用旧值否则迭代就错了。第三角度预处理只处理了负角度的情况对于大于90°的正角度需要额外判断。3.2 象限预处理到底怎么处理CORDIC旋转模式的收敛范围是±99.9°但实际应用中角度可能是0到360°。处理方法是利用三角函数的对称性如果角度在90°到180°之间sin(θ) sin(180°-θ)cos(θ) -cos(180°-θ)如果角度在180°到270°之间sin(θ) -sin(θ-180°)cos(θ) -cos(θ-180°)如果角度在270°到360°之间sin(θ) -sin(360°-θ)cos(θ) cos(360°-θ)在定点数表示里Q2.14格式下180°对应0x8000也就是-3276890°对应0x400016384。判断象限就是看角度值的高两位。更简洁的做法是先把角度归一化到-90°到90°之间记录象限信息迭代完成后再根据象限修正符号。这样CORDIC核心只需要处理-90°到90°的范围完全在收敛范围内。3.3 流水线版本的关键改动流水线版本把16次迭代展开成16级每级用独立的寄存器。这样每个时钟周期都能接收新角度并吐出新结果。代价是寄存器数量大约是迭代式的16倍但在Artix-7上这点资源完全不是问题。流水线版本的核心改动是每一级都有自己独立的x、y、z寄存器级与级之间用寄存器打拍。角度常数表变成组合逻辑每级直接查表。方向判断在每一级独立进行因为每级的z值不同。// 流水线单级示例 module cordic_stage ( input wire clk, input wire [15:0] x_in, input wire [15:0] y_in, input wire [15:0] z_in, input wire [3:0] stage, output reg [15:0] x_out, output reg [15:0] y_out, output reg [15:0] z_out ); wire [15:0] atan_val atan_table(stage); wire dir ~z_in[15]; // z 0 时逆时针 always (posedge clk) begin if (dir) begin x_out x_in - (y_in stage); y_out y_in (x_in stage); z_out z_in - atan_val; end else begin x_out x_in (y_in stage); y_out y_in - (x_in stage); z_out z_in atan_val; end end endmodule流水线版本上板后你可以用板载的50MHz时钟每个周期算一个角度理论上每秒能算5000万个三角函数值。当然实际输出到数码管或者DAC的时候瓶颈在输出接口上。4. EGo1上板验证从代码到看得见的结果4.1 板卡资源分配与约束文件EGo1的时钟是100MHz板载晶振但CORDIC不需要跑那么快。我一般用分频后的1Hz到1kHz来驱动角度输入这样数码管刷新和肉眼观察都方便。约束文件里需要绑定几个关键信号时钟输入EGo1的100MHz时钟在E3引脚具体看板卡原理图复位按键通常用BTN0或者SW0数码管段选和位选EGo1用的是共阳极数码管段选低电平有效LED用来指示done信号或者溢出# EGo1约束示例部分 set_property PACKAGE_PIN E3 [get_ports clk_100m] set_property IOSTANDARD LVCMOS33 [get_ports clk_100m] set_property PACKAGE_PIN D9 [get_ports rst_n] set_property IOSTANDARD LVCMOS33 [get_ports rst_n] # 数码管段选 set_property PACKAGE_PIN B4 [get_ports {seg[0]}] set_property PACKAGE_PIN A4 [get_ports {seg[1]}] # ... 其余段选引脚4.2 用数码管显示sin和cos值数码管显示定点数需要做二进制到BCD的转换。16位定点数Q2.14实际值范围是-2到2但sin和cos的范围是-1到1。显示的时候我通常把结果乘以10000然后显示整数部分。比如cos_out 0x6964这是1.64676的定点表示。但实际cos值应该在-1到1之间。等等这里有个容易搞混的地方CORDIC迭代完成后输出的x就是cos值但它的定点格式和输入角度不一样。输入角度是Q2.14范围-180到180输出cos/sin是Q1.14范围-1到1。所以0x6964对应的实际值是0x6964 / 2^14 1.64676这不对cos不可能大于1。问题出在哪初始x设的是1/K 1.64676但迭代完成后x会缩小K倍最终结果才是cos值。所以如果你初始x设1/K迭代完直接输出x就是cos。但如果你初始x设1.0迭代完输出的是K·cos需要再乘1/K修正。我建议初始x直接设1/K的定点值这样输出就是真正的cos和sin不需要额外修正。验证一下当角度为0时cos应该是1.0对应定点0x4000Q1.14格式下1.0 16384 0x4000。你可以用这个作为测试用例。4.3 实测中遇到的三个坑第一个坑右移的符号问题。Verilog里是逻辑右移是算术右移。对于有符号数必须用否则负数右移会变成正数。我一开始用了结果角度在负半轴时输出完全乱套。第二个坑迭代次数和精度的关系。我一开始只迭代了8次想着省点资源结果cos(45°)算出来误差有0.01左右。后来加到12次误差降到0.001以内。16次迭代的误差在0.0001量级对于数码管显示来说完全够用。但如果你要做DDS或者通信调制可能需要更多迭代次数。第三个坑角度输入的格式。EGo1上的拨码开关或者按键输入的是整数你需要把它转换成Q2.14格式的角度值。比如输入45°对应的定点值是45/180 × 32768 8192 0x2000。这个转换在Testbench里容易搞错上板前一定要在仿真里验证几个关键角度。角度Q2.14定点值理论cos理论sin0°0x00001.00000.000030°0x0AAA0.86600.500045°0x20000.70710.707160°0x2AAA0.50000.866090°0x40000.00001.00005. 精度、资源与速度的三角平衡5.1 迭代次数对精度的影响到底有多大CORDIC的精度主要受两个因素影响迭代次数和定点数位宽。理论上N次迭代能得到大约N位的精度但实际上因为角度常数的量化误差和累加误差有效精度会略低。我做过一组实测用16位定点数分别迭代8、12、16次然后跟MATLAB的double精度结果对比。8次迭代的最大误差约0.00512次约0.000816次约0.0001。对于数码管显示4位有效数字12次迭代就够了。对于音频信号生成16位DAC16次迭代是底线。如果你需要更高精度有两个方向增加迭代次数或者增加位宽。增加迭代次数的边际效益递减因为后面的atan值越来越小对结果的修正越来越微弱。增加位宽更有效但资源消耗线性增长。5.2 资源占用实测数据在Artix-7 XC7A35T上综合后的资源占用配置LUTFFDSPBlock RAM迭代式16次31219800流水线16级1847112000流水线12级142386400可以看到CORDIC完全不消耗DSP和Block RAM这是它最大的优势。流水线版本的LUT和FF用量大约是迭代式的6倍但对于XC7A35T约20800个LUT来说即使流水线版本也只用了不到10%的资源。5.3 什么时候该用CORDIC什么时候该用查找表这个问题没有标准答案但有一个经验法则如果你的精度要求不超过10位查找表更简单直接如果精度要求12位以上CORDIC的资源效率更高。查找表的资源消耗随精度指数增长。10位精度需要1024个条目每个条目16位就是16Kb的ROM。12位精度就是64Kb16位精度就是1Mb。而CORDIC的资源消耗随精度线性增长16位精度和32位精度的资源差距只有一倍左右。另一个考虑因素是灵活性。查找表只能输出预先存好的值如果你想动态改变频率或者相位查找表需要重新加载。CORDIC是纯计算输入什么角度就算什么天然支持任意角度。6. 从单点计算到连续信号生成6.1 相位累加器与CORDIC的配合单独算一个角度的sin和cos意义不大实际应用中通常是连续生成正弦波。这时候需要一个相位累加器每个时钟周期累加一个频率控制字累加器的输出作为CORDIC的角度输入。相位累加器的位宽决定了频率分辨率。比如32位累加器时钟100MHz频率分辨率是100MHz / 2^32 ≈ 0.023Hz。频率控制字K对应的输出频率是K × 100MHz / 2^32。但CORDIC的角度输入是16位所以需要把32位相位截断到16位。截断会引入相位噪声这是DDS的固有缺陷。解决办法是用相位抖动或者更高位宽的CORDIC但那是另一个话题了。6.2 输出到DAC的接口设计EGo1板载没有高速DAC但你可以用PWM加低通滤波器的方式输出模拟信号。把CORDIC算出的sin值Q1.14格式取高8位作为PWM的占空比PWM频率设到100kHz以上经过RC低通滤波就能得到还算干净的正弦波。PWM的位宽决定了输出信号的SNR。8位PWM的理论SNR约50dB对于示波器观察来说足够了。如果你想要更好的效果可以用外接的并行DAC比如TLC7528或者AD9708。// PWM生成示例 reg [7:0] pwm_cnt; reg pwm_out; always (posedge clk) begin pwm_cnt pwm_cnt 1b1; pwm_out (pwm_cnt sin_out[14:7]) ? 1b1 : 1b0; end6.3 实测波形与误差分析我用EGo1的PWM输出接了一个简单的RC滤波器1kΩ 100nF截止频率约1.6kHz。生成1kHz正弦波时示波器上看波形很干净THD大约在-40dB左右。主要失真来源是PWM的量化噪声和RC滤波器的非线性。如果你把CORDIC的输出直接接到逻辑分析仪或者ILA上可以看到sin和cos的数值序列。在1kHz输出时每个周期有1000个采样点假设采样率1MHz波形非常平滑。把ILA的数据导出到MATLAB做FFT可以看到除了基波之外只有很小的谐波分量。注意PWM输出时sin_out是有符号数需要先加上偏置变成无符号数才能作为占空比。具体做法是sin_out 0x4000然后取高位。7. 一些容易忽略但很重要的细节7.1 复位后的第一个结果不可信CORDIC迭代式实现里复位后x、y、z的初始值需要几个周期才能稳定。如果你在start信号拉高后立刻读结果读到的可能是中间状态。我的做法是等done信号拉高后再读done信号在最后一次迭代完成后的下一个周期拉高。流水线版本也有类似问题。流水线填满需要16个周期前16个输出是无效的。你需要一个valid信号在流水线填满后拉高表示输出有效。7.2 角度常数的定点表示要统一atan_table里的常数是Q2.14格式但x和y是Q1.14格式。z的更新是z ± atan_val两者都是Q2.14没问题。但x和y的更新涉及移位移位不改变定点格式。所以整个迭代过程中x和y始终保持Q1.14z始终保持Q2.14。输出的时候x和y直接就是Q1.14的cos和sin。如果你在代码里混用了不同格式的定点数结果会差一个2的幂次。这种错误在仿真里很容易发现因为输出值会大得离谱或者小得离谱。7.3 Testbench怎么写才能覆盖所有情况一个好的Testbench应该覆盖0°、90°、180°、270°这些边界角度正负45°这些典型角度以及接近收敛边界的角度比如89°和91°。每个角度都要跟MATLAB或者Python的计算结果对比。// Testbench片段 initial begin // 测试0° angle_in 16h0000; start 1b1; #20 start 1b0; wait(done); #10; $display(0°: cos%d, sin%d, cos_out, sin_out); // 测试45° angle_in 16h2000; start 1b1; #20 start 1b0; wait(done); #10; $display(45°: cos%d, sin%d, cos_out, sin_out); // 测试-45° angle_in 16hE000; // -45°的补码 start 1b1; #20 start 1b0; wait(done); #10; $display(-45°: cos%d, sin%d, cos_out, sin_out); end预期结果0°时cos0x40001.0sin0x000045°时cossin0x2D410.7071-45°时cos0x2D41sin0xD2BF-0.7071。7.4 上板调试时ILA的用法Vivado的ILA是调试CORDIC的利器。把x、y、z、iter_cnt、done都接到ILA上触发条件设为start的上升沿。这样你可以看到整个迭代过程每一级的x、y、z变化一目了然。我一般会观察z的收敛情况如果z在最后几次迭代时还在大幅变化说明迭代次数不够如果z很快就趋近于0然后保持不变说明迭代次数有富余。理想情况下z应该在最后一次迭代后接近0但又不完全为0。ILA的采样深度要设够至少能覆盖一次完整的迭代过程16个周期加上前后的空闲周期。采样深度1024对于迭代式足够了流水线版本需要2048以上。8. 从EGo1到实际项目的移植建议EGo1验证通过之后如果你要把CORDIC移植到实际项目里有几个地方需要调整。首先是时钟频率实际项目可能跑在200MHz甚至更高这时候组合逻辑的延迟会成为瓶颈。解决办法是在流水线中间插入额外的寄存器或者降低单级逻辑的复杂度。其次是位宽实际项目可能需要18位、24位甚至32位精度。位宽增加后角度常数表也要相应扩展。32位精度的CORDIC在Artix-7上大约消耗5000个LUT仍然可以接受。最后是接口实际项目可能用AXI-Stream或者自定义的并行接口。CORDIC核心本身是流式的加一个简单的握手信号就能适配大多数接口。我个人在多个项目里用过CORDIC从简单的信号发生器到复杂的通信调制解调它从来没让我失望过。唯一需要注意的是CORDIC的精度和迭代次数是硬绑定的你不能指望16次迭代给出20位的精度。在项目初期就把精度需求定清楚后面会省很多事。 FPGA实战：手写Verilog实现CORDIC算法计算三角函数
1. 为什么要在FPGA里用CORDIC算三角函数很多人第一次接触在FPGA上算sin和cos第一反应是去找IP核或者干脆用查找表。查找表确实简单但精度和资源是一对矛盾——你要16位精度ROM就得占一大块Block RAM你要省资源精度就惨不忍睹。而CORDIC这个算法有意思的地方在于它压根不用乘法器只靠移位和加减就能把三角函数算出来。这在FPGA里简直是天作之合因为移位在硬件里就是连线加法器是FPGA最不缺的东西。CORDIC全称是Coordinate Rotation Digital Computer坐标旋转数字计算机。名字听着唬人核心思想其实很朴素把一个角度不断分解成一系列越来越小的固定角度每次判断当前角度是正还是负决定往哪个方向转。转着转着角度就逼近目标了同时坐标也跟着转到了目标位置。旋转模式下如果你把初始向量放在x轴上x1, y0那么旋转到目标角度θ之后新的x坐标就是cos(θ)y坐标就是sin(θ)。这就是旋转模式算三角函数的本质。我选择EGo1这块板卡来做验证原因很实际。EGo1是Xilinx Artix-7系列的入门级开发板资源不算富裕但足够跑CORDIC板载的时钟、按键、LED和数码管刚好能构成一个完整的验证闭环。你不需要额外的示波器或者逻辑分析仪用板载资源就能看到结果对不对。对于学习CORDIC来说这种自闭环的验证方式比仿真波形直观得多。这篇文章适合谁看如果你已经写过Verilog或者VHDL知道什么是时序逻辑和组合逻辑但还没真正把CORDIC跑通过那这篇就是写给你的。我会从迭代公式的推导讲起把每一级流水线为什么这么设计说清楚然后给出完整的Verilog代码最后讲上板验证时怎么用数码管把结果展示出来。整个过程不需要任何IP核纯手写RTL。提示CORDIC的迭代次数和精度是直接挂钩的。N次迭代大约能得到N位的精度但角度覆盖范围有限后面会详细说这个问题怎么处理。2. CORDIC旋转模式的迭代公式到底怎么来的2.1 从一次旋转说起假设平面上有一个点(x, y)我们要把它绕原点旋转角度θ。根据旋转矩阵x x·cos(θ) - y·sin(θ) y x·sin(θ) y·cos(θ)这个公式没问题但问题在于cos(θ)和sin(θ)本身就是我们要算的东西这就成了循环依赖。CORDIC的巧妙之处在于它把旋转角θ拆成一堆预先定好的小角度每个小角度都有一个特点它的正切值是2的负整数次幂。具体来说第i次旋转的角度是atan(2^(-i))。为什么选这个角度因为tan(atan(2^(-i))) 2^(-i)而乘以2^(-i)在硬件里就是右移i位。这样一来旋转矩阵里的cos和sin就被消掉了只剩下移位和加法。把cos(θ)从旋转矩阵里提出来x cos(θ)·(x - y·tan(θ)) y cos(θ)·(y x·tan(θ))当tan(θ) ±2^(-i)时括号里的运算就变成了x cos(θ)·(x - y·d·2^(-i)) y cos(θ)·(y x·d·2^(-i))其中d是方向因子1表示逆时针转-1表示顺时针转。2.2 那个cos(θ)乘积去哪了每次迭代都有一个cos(atan(2^(-i)))的乘积因子。N次迭代下来总的乘积是K ∏ cos(atan(2^(-i))) (i从0到N-1)这个K是一个常数当N足够大时K约等于0.607252935。也就是说如果初始向量是(1, 0)经过N次迭代后得到的向量长度是K而不是1。所以如果你想要真正的cos和sin值初始x应该设为1/K也就是约1.646760258。在硬件实现里通常直接把初始x设为1/K的定点数表示。比如16位定点数1位符号1位整数14位小数1/K ≈ 1.64676对应的定点值就是1.64676 × 2^14 ≈ 26980十六进制是0x6964。2.3 角度累加器的作用每次迭代还需要一个角度累加器z初始值是目标角度θ。每次迭代根据z的符号决定旋转方向如果z 0说明当前角度还不够需要逆时针转d 1如果z 0说明转多了需要顺时针转d -1然后z更新为z - d·atan(2^(-i))。迭代足够多次后z趋近于0此时(x, y)就是(cos(θ), sin(θ))。这里有个细节atan(2^(-i))的值需要预先算好存起来。对于16位精度通常迭代16次需要16个角度常数。这些常数用定点数表示存在一个case语句或者ROM里。迭代次数iatan(2^(-i))角度值定点表示Q2.14045.000°0x2000126.565°0x12E4214.036°0x09FB37.125°0x051143.576°0x028B51.789°0x014660.895°0x00A370.448°0x005180.224°0x002990.112°0x0014100.056°0x000A110.028°0x0005120.014°0x0003130.007°0x0001140.003°0x0001150.002°0x0000注意角度范围问题CORDIC旋转模式能覆盖的角度范围是-99.9°到99.9°左右因为所有atan(2^(-i))的和约等于99.9°。如果你要算的角度超过这个范围需要先做象限预处理。比如算120°的sin可以先算60°的sin然后根据象限关系转换。这个后面在代码里会处理。3. Verilog实现从单周期迭代到全流水线3.1 迭代式实现与流水线式的取舍最直观的写法是用一个状态机每个时钟周期做一次迭代16次迭代就是16个周期出结果。这种写法资源省但吞吐率低。另一种写法是全流水线16级流水线排开每个时钟周期都能吐出一个新结果代价是寄存器用量大。在EGo1这种入门板上两种写法都能跑。我建议先用迭代式把功能调通因为迭代式调试起来简单你可以把中间每一级的x、y、z都接到ILA上观察。等确认算法没问题了再改成流水线提升吞吐率。迭代式的核心代码结构大概是这样module cordic_iterative ( input wire clk, input wire rst_n, input wire start, input wire [15:0] angle_in, // Q2.14格式范围-180°到180° output reg [15:0] cos_out, output reg [15:0] sin_out, output reg done ); // 角度常数表Q2.14格式 function [15:0] atan_table; input [3:0] idx; case (idx) 4d0: atan_table 16h2000; // 45.000° 4d1: atan_table 16h12E4; // 26.565° 4d2: atan_table 16h09FB; // 14.036° 4d3: atan_table 16h0511; // 7.125° 4d4: atan_table 16h028B; // 3.576° 4d5: atan_table 16h0146; // 1.789° 4d6: atan_table 16h00A3; // 0.895° 4d7: atan_table 16h0051; // 0.448° 4d8: atan_table 16h0029; // 0.224° 4d9: atan_table 16h0014; // 0.112° 4d10: atan_table 16h000A; // 0.056° 4d11: atan_table 16h0005; // 0.028° 4d12: atan_table 16h0003; // 0.014° 4d13: atan_table 16h0001; // 0.007° 4d14: atan_table 16h0001; // 0.003° 4d15: atan_table 16h0000; // 0.002° default: atan_table 16h0000; endcase endfunction reg [4:0] iter_cnt; reg [15:0] x, y, z; reg busy; // 象限预处理 reg [15:0] angle_pre; reg sign_flip; always (posedge clk or negedge rst_n) begin if (!rst_n) begin iter_cnt 5d0; busy 1b0; done 1b0; x 16d0; y 16d0; z 16d0; end else begin done 1b0; if (start !busy) begin busy 1b1; iter_cnt 5d0; // 初始值x 1/K ≈ 1.64676Q2.14格式 x 16h6964; y 16d0; // 角度预处理如果角度绝对值大于90°做象限转换 if (angle_in[15]) begin // 负角度取反加一得到绝对值 z ~angle_in 1b1; sign_flip 1b1; end else begin z angle_in; sign_flip 1b0; end end else if (busy) begin if (iter_cnt 5d15) begin busy 1b0; done 1b1; // 根据符号决定输出 if (sign_flip) begin cos_out x; sin_out ~y 1b1; // sin(-θ) -sin(θ) end else begin cos_out x; sin_out y; end end else begin iter_cnt iter_cnt 1b1; if (z[15] 1b0) begin // z 0逆时针转 x x - (y iter_cnt); y y (x iter_cnt); z z - atan_table(iter_cnt[3:0]); end else begin // z 0顺时针转 x x (y iter_cnt); y y - (x iter_cnt); z z atan_table(iter_cnt[3:0]); end end end end end endmodule这段代码有几个地方需要特别注意。第一是算术右移对于有符号数才能保持符号位。第二x和y的更新必须用同一时刻的旧值不能一个用新值一个用旧值否则迭代就错了。第三角度预处理只处理了负角度的情况对于大于90°的正角度需要额外判断。3.2 象限预处理到底怎么处理CORDIC旋转模式的收敛范围是±99.9°但实际应用中角度可能是0到360°。处理方法是利用三角函数的对称性如果角度在90°到180°之间sin(θ) sin(180°-θ)cos(θ) -cos(180°-θ)如果角度在180°到270°之间sin(θ) -sin(θ-180°)cos(θ) -cos(θ-180°)如果角度在270°到360°之间sin(θ) -sin(360°-θ)cos(θ) cos(360°-θ)在定点数表示里Q2.14格式下180°对应0x8000也就是-3276890°对应0x400016384。判断象限就是看角度值的高两位。更简洁的做法是先把角度归一化到-90°到90°之间记录象限信息迭代完成后再根据象限修正符号。这样CORDIC核心只需要处理-90°到90°的范围完全在收敛范围内。3.3 流水线版本的关键改动流水线版本把16次迭代展开成16级每级用独立的寄存器。这样每个时钟周期都能接收新角度并吐出新结果。代价是寄存器数量大约是迭代式的16倍但在Artix-7上这点资源完全不是问题。流水线版本的核心改动是每一级都有自己独立的x、y、z寄存器级与级之间用寄存器打拍。角度常数表变成组合逻辑每级直接查表。方向判断在每一级独立进行因为每级的z值不同。// 流水线单级示例 module cordic_stage ( input wire clk, input wire [15:0] x_in, input wire [15:0] y_in, input wire [15:0] z_in, input wire [3:0] stage, output reg [15:0] x_out, output reg [15:0] y_out, output reg [15:0] z_out ); wire [15:0] atan_val atan_table(stage); wire dir ~z_in[15]; // z 0 时逆时针 always (posedge clk) begin if (dir) begin x_out x_in - (y_in stage); y_out y_in (x_in stage); z_out z_in - atan_val; end else begin x_out x_in (y_in stage); y_out y_in - (x_in stage); z_out z_in atan_val; end end endmodule流水线版本上板后你可以用板载的50MHz时钟每个周期算一个角度理论上每秒能算5000万个三角函数值。当然实际输出到数码管或者DAC的时候瓶颈在输出接口上。4. EGo1上板验证从代码到看得见的结果4.1 板卡资源分配与约束文件EGo1的时钟是100MHz板载晶振但CORDIC不需要跑那么快。我一般用分频后的1Hz到1kHz来驱动角度输入这样数码管刷新和肉眼观察都方便。约束文件里需要绑定几个关键信号时钟输入EGo1的100MHz时钟在E3引脚具体看板卡原理图复位按键通常用BTN0或者SW0数码管段选和位选EGo1用的是共阳极数码管段选低电平有效LED用来指示done信号或者溢出# EGo1约束示例部分 set_property PACKAGE_PIN E3 [get_ports clk_100m] set_property IOSTANDARD LVCMOS33 [get_ports clk_100m] set_property PACKAGE_PIN D9 [get_ports rst_n] set_property IOSTANDARD LVCMOS33 [get_ports rst_n] # 数码管段选 set_property PACKAGE_PIN B4 [get_ports {seg[0]}] set_property PACKAGE_PIN A4 [get_ports {seg[1]}] # ... 其余段选引脚4.2 用数码管显示sin和cos值数码管显示定点数需要做二进制到BCD的转换。16位定点数Q2.14实际值范围是-2到2但sin和cos的范围是-1到1。显示的时候我通常把结果乘以10000然后显示整数部分。比如cos_out 0x6964这是1.64676的定点表示。但实际cos值应该在-1到1之间。等等这里有个容易搞混的地方CORDIC迭代完成后输出的x就是cos值但它的定点格式和输入角度不一样。输入角度是Q2.14范围-180到180输出cos/sin是Q1.14范围-1到1。所以0x6964对应的实际值是0x6964 / 2^14 1.64676这不对cos不可能大于1。问题出在哪初始x设的是1/K 1.64676但迭代完成后x会缩小K倍最终结果才是cos值。所以如果你初始x设1/K迭代完直接输出x就是cos。但如果你初始x设1.0迭代完输出的是K·cos需要再乘1/K修正。我建议初始x直接设1/K的定点值这样输出就是真正的cos和sin不需要额外修正。验证一下当角度为0时cos应该是1.0对应定点0x4000Q1.14格式下1.0 16384 0x4000。你可以用这个作为测试用例。4.3 实测中遇到的三个坑第一个坑右移的符号问题。Verilog里是逻辑右移是算术右移。对于有符号数必须用否则负数右移会变成正数。我一开始用了结果角度在负半轴时输出完全乱套。第二个坑迭代次数和精度的关系。我一开始只迭代了8次想着省点资源结果cos(45°)算出来误差有0.01左右。后来加到12次误差降到0.001以内。16次迭代的误差在0.0001量级对于数码管显示来说完全够用。但如果你要做DDS或者通信调制可能需要更多迭代次数。第三个坑角度输入的格式。EGo1上的拨码开关或者按键输入的是整数你需要把它转换成Q2.14格式的角度值。比如输入45°对应的定点值是45/180 × 32768 8192 0x2000。这个转换在Testbench里容易搞错上板前一定要在仿真里验证几个关键角度。角度Q2.14定点值理论cos理论sin0°0x00001.00000.000030°0x0AAA0.86600.500045°0x20000.70710.707160°0x2AAA0.50000.866090°0x40000.00001.00005. 精度、资源与速度的三角平衡5.1 迭代次数对精度的影响到底有多大CORDIC的精度主要受两个因素影响迭代次数和定点数位宽。理论上N次迭代能得到大约N位的精度但实际上因为角度常数的量化误差和累加误差有效精度会略低。我做过一组实测用16位定点数分别迭代8、12、16次然后跟MATLAB的double精度结果对比。8次迭代的最大误差约0.00512次约0.000816次约0.0001。对于数码管显示4位有效数字12次迭代就够了。对于音频信号生成16位DAC16次迭代是底线。如果你需要更高精度有两个方向增加迭代次数或者增加位宽。增加迭代次数的边际效益递减因为后面的atan值越来越小对结果的修正越来越微弱。增加位宽更有效但资源消耗线性增长。5.2 资源占用实测数据在Artix-7 XC7A35T上综合后的资源占用配置LUTFFDSPBlock RAM迭代式16次31219800流水线16级1847112000流水线12级142386400可以看到CORDIC完全不消耗DSP和Block RAM这是它最大的优势。流水线版本的LUT和FF用量大约是迭代式的6倍但对于XC7A35T约20800个LUT来说即使流水线版本也只用了不到10%的资源。5.3 什么时候该用CORDIC什么时候该用查找表这个问题没有标准答案但有一个经验法则如果你的精度要求不超过10位查找表更简单直接如果精度要求12位以上CORDIC的资源效率更高。查找表的资源消耗随精度指数增长。10位精度需要1024个条目每个条目16位就是16Kb的ROM。12位精度就是64Kb16位精度就是1Mb。而CORDIC的资源消耗随精度线性增长16位精度和32位精度的资源差距只有一倍左右。另一个考虑因素是灵活性。查找表只能输出预先存好的值如果你想动态改变频率或者相位查找表需要重新加载。CORDIC是纯计算输入什么角度就算什么天然支持任意角度。6. 从单点计算到连续信号生成6.1 相位累加器与CORDIC的配合单独算一个角度的sin和cos意义不大实际应用中通常是连续生成正弦波。这时候需要一个相位累加器每个时钟周期累加一个频率控制字累加器的输出作为CORDIC的角度输入。相位累加器的位宽决定了频率分辨率。比如32位累加器时钟100MHz频率分辨率是100MHz / 2^32 ≈ 0.023Hz。频率控制字K对应的输出频率是K × 100MHz / 2^32。但CORDIC的角度输入是16位所以需要把32位相位截断到16位。截断会引入相位噪声这是DDS的固有缺陷。解决办法是用相位抖动或者更高位宽的CORDIC但那是另一个话题了。6.2 输出到DAC的接口设计EGo1板载没有高速DAC但你可以用PWM加低通滤波器的方式输出模拟信号。把CORDIC算出的sin值Q1.14格式取高8位作为PWM的占空比PWM频率设到100kHz以上经过RC低通滤波就能得到还算干净的正弦波。PWM的位宽决定了输出信号的SNR。8位PWM的理论SNR约50dB对于示波器观察来说足够了。如果你想要更好的效果可以用外接的并行DAC比如TLC7528或者AD9708。// PWM生成示例 reg [7:0] pwm_cnt; reg pwm_out; always (posedge clk) begin pwm_cnt pwm_cnt 1b1; pwm_out (pwm_cnt sin_out[14:7]) ? 1b1 : 1b0; end6.3 实测波形与误差分析我用EGo1的PWM输出接了一个简单的RC滤波器1kΩ 100nF截止频率约1.6kHz。生成1kHz正弦波时示波器上看波形很干净THD大约在-40dB左右。主要失真来源是PWM的量化噪声和RC滤波器的非线性。如果你把CORDIC的输出直接接到逻辑分析仪或者ILA上可以看到sin和cos的数值序列。在1kHz输出时每个周期有1000个采样点假设采样率1MHz波形非常平滑。把ILA的数据导出到MATLAB做FFT可以看到除了基波之外只有很小的谐波分量。注意PWM输出时sin_out是有符号数需要先加上偏置变成无符号数才能作为占空比。具体做法是sin_out 0x4000然后取高位。7. 一些容易忽略但很重要的细节7.1 复位后的第一个结果不可信CORDIC迭代式实现里复位后x、y、z的初始值需要几个周期才能稳定。如果你在start信号拉高后立刻读结果读到的可能是中间状态。我的做法是等done信号拉高后再读done信号在最后一次迭代完成后的下一个周期拉高。流水线版本也有类似问题。流水线填满需要16个周期前16个输出是无效的。你需要一个valid信号在流水线填满后拉高表示输出有效。7.2 角度常数的定点表示要统一atan_table里的常数是Q2.14格式但x和y是Q1.14格式。z的更新是z ± atan_val两者都是Q2.14没问题。但x和y的更新涉及移位移位不改变定点格式。所以整个迭代过程中x和y始终保持Q1.14z始终保持Q2.14。输出的时候x和y直接就是Q1.14的cos和sin。如果你在代码里混用了不同格式的定点数结果会差一个2的幂次。这种错误在仿真里很容易发现因为输出值会大得离谱或者小得离谱。7.3 Testbench怎么写才能覆盖所有情况一个好的Testbench应该覆盖0°、90°、180°、270°这些边界角度正负45°这些典型角度以及接近收敛边界的角度比如89°和91°。每个角度都要跟MATLAB或者Python的计算结果对比。// Testbench片段 initial begin // 测试0° angle_in 16h0000; start 1b1; #20 start 1b0; wait(done); #10; $display(0°: cos%d, sin%d, cos_out, sin_out); // 测试45° angle_in 16h2000; start 1b1; #20 start 1b0; wait(done); #10; $display(45°: cos%d, sin%d, cos_out, sin_out); // 测试-45° angle_in 16hE000; // -45°的补码 start 1b1; #20 start 1b0; wait(done); #10; $display(-45°: cos%d, sin%d, cos_out, sin_out); end预期结果0°时cos0x40001.0sin0x000045°时cossin0x2D410.7071-45°时cos0x2D41sin0xD2BF-0.7071。7.4 上板调试时ILA的用法Vivado的ILA是调试CORDIC的利器。把x、y、z、iter_cnt、done都接到ILA上触发条件设为start的上升沿。这样你可以看到整个迭代过程每一级的x、y、z变化一目了然。我一般会观察z的收敛情况如果z在最后几次迭代时还在大幅变化说明迭代次数不够如果z很快就趋近于0然后保持不变说明迭代次数有富余。理想情况下z应该在最后一次迭代后接近0但又不完全为0。ILA的采样深度要设够至少能覆盖一次完整的迭代过程16个周期加上前后的空闲周期。采样深度1024对于迭代式足够了流水线版本需要2048以上。8. 从EGo1到实际项目的移植建议EGo1验证通过之后如果你要把CORDIC移植到实际项目里有几个地方需要调整。首先是时钟频率实际项目可能跑在200MHz甚至更高这时候组合逻辑的延迟会成为瓶颈。解决办法是在流水线中间插入额外的寄存器或者降低单级逻辑的复杂度。其次是位宽实际项目可能需要18位、24位甚至32位精度。位宽增加后角度常数表也要相应扩展。32位精度的CORDIC在Artix-7上大约消耗5000个LUT仍然可以接受。最后是接口实际项目可能用AXI-Stream或者自定义的并行接口。CORDIC核心本身是流式的加一个简单的握手信号就能适配大多数接口。我个人在多个项目里用过CORDIC从简单的信号发生器到复杂的通信调制解调它从来没让我失望过。唯一需要注意的是CORDIC的精度和迭代次数是硬绑定的你不能指望16次迭代给出20位的精度。在项目初期就把精度需求定清楚后面会省很多事。