ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

FPGA核脉冲梯形成形算法:从MATLAB仿真到Verilog工程实现

FPGA核脉冲梯形成形算法:从MATLAB仿真到Verilog工程实现 1. 这块硬骨头到底在解决什么问题做核仪器、核辐射探测或者能谱测量这块的工程师应该对“核脉冲成形”这几个字不陌生。探测器出来的信号不管是NaI闪烁体加光电倍增管、HPGe还是硅探测器经过前置放大之后基本都是同一个毛病指数衰减尾巴特别长。用示波器看就是一个陡峭上升然后慢慢拖下来的脉冲时间常数在微秒量级。这样的信号直接送去测幅度、做能谱问题非常大。脉冲堆积、弹道亏损、基线漂移随便一个就能把你的能量分辨率搞得很难看。梯形成形就是干这个用的。它能把指数衰减的脉冲通过数字滤波器变换成一个梯形——上升沿、平顶、下降沿都特别规整。为什么要费劲做这一层变换因为梯形脉冲的幅度和平顶采样值能非常好地反映输入脉冲的能量信息同时对弹道亏损不敏感对堆积也有一定的容忍度。而FPGA配合高速ADC做实时数字脉冲处理又是目前核信号处理的主流路线。所以这篇博客就围绕一个很具体的工程问题展开怎么把MATLAB里验证好的梯形成形算法一步一步搬到Verilog上最终在50MHz采样率下稳定跑起来。这套东西适合正在做核脉冲处理、数字能谱仪、辐射监测、PET读出电子学的同学参考也适合刚入门FPGA数字信号处理、想找一个不太复杂又能学到真东西的工程项目的朋友。我尽量把从算法原理到代码实现再到上板调试的整个链路都讲透包括位宽怎么算、延迟线怎么做、流水线怎么排、上板之后会遇到哪些坑。文章里的代码和思路都是实际工程里验证过的不是demo级别的玩具。2. 算法原理和MATLAB验证先把数学模型打通2.1 从指数脉冲到梯形递推公式是怎么来的先说数学模型。设ADC采样后得到的离散输入序列为 x[n]它近似是一个单指数衰减信号也就是 x[n] A · e^(-n/τ) 这种形态。梯形成形的目标是输出一个 y[n]波形形状为上升 k 个点平顶 l 个点再下降 k 个点。这个变换在z域上看非常优雅。梯形成形滤波器的传递函数可以写成H(z) [(1 - z^(-k)) · (1 - z^(-k-l))] / [(1 - z^(-1)) · (1 - z^(-1))]分子上是两个差分器分母上是两个累加器。差分器负责把信号的前后沿“切”出来累加器再把切出来的片段“拼”起来最后拼出一个梯形轮廓。把这个传递函数展开成时域递推表达式就是大家经常看到的经典算法d[n] x[n] - x[n-k] - x[n-k-l] x[n-2k-l]p[n] p[n-1] d[n]y[n] y[n-1] p[n]其中 k 是上升沿对应的采样点数l 是平顶宽度对应的采样点数两个参数设置的单位都是“采样点个数”。我换成人话解释一下。你可以把输入信号想象成一个人从远处跑过来速度很快然后慢悠悠走过去。梯形成形做的事情是给这个人设了三道闸门第一道闸门让前 k 个点进来形成上升沿中间 l 个点保持队伍长度不变形成平顶最后 k 个点再放走形成下降沿。差分器负责计算“进来了多少人”和“走了多少人”的差累加器负责把人头数累计起来最后得到的就是队伍的实时长度——也就是梯形的形状。这个递推式的实现成本特别低。每个输出点只需要做一次4项加减法和两次累加没有任何乘法在FPGA里用纯加法器流水线就能跑得飞快。这也是为什么梯形成形算法在核脉冲处理领域被用到烂——不是因为它新颖而是因为它便宜、稳定、效果好。有一点要提醒上面公式在序列起点附近会有边界问题。比如当 n k 的时候x[n-k] 是不存在的。标准做法是把输入序列前面补零让长度大于 2kl 再开始算否则开头的输出会有一段畸形。MATLAB仿真阶段就得把这个细节处理好不然后面移植到Verilogtestbench对不上号排查起来很痛苦。2.2 MATLAB浮点仿真先把算法跑通在碰Verilog之前我强烈建议先用MATLAB把浮点算法完整跑一遍。原因很简单FPGA上调试信号看不到出了问题定位成本极高。而在MATLAB里你前后左右想怎么画图就怎么画图算法逻辑对不对一目了然。第一步是生成模拟的核脉冲数据。真实场景下核事件到达是随机的近似满足泊松分布。我做仿真时习惯用这样的方式生成测试数据先随机生成一批脉冲的到达时间和幅度每个脉冲按指数衰减模型叠加再加上高斯白噪声。% 生成模拟核脉冲序列 fs 50e6; % 采样率 50MHz T 1/fs; N 20000; % 采样点数 tau 5e-6; % 指数衰减时间常数对应5us decay exp(-T/tau); x zeros(1, N); % 随机脉冲平均每50个采样点来一个脉冲 pulse_pos sort(randperm(N, floor(N/50))); pulse_amp rand(1, length(pulse_pos)) * 30000 5000; for i 1:length(pulse_pos) n0 pulse_pos(i); len min(500, N - n0); % 500点以内的拖尾 n 0:len-1; x(n0:n0len-1) x(n0:n0len-1) pulse_amp(i) * decay.^n; end x x randn(1, N) * 30; % 加噪声然后写一个梯形成形的纯循环版本先把浮点结果跑出来。function y trapezoid_shaping(x, k, l) N length(x); d zeros(1, N); p zeros(1, N); y zeros(1, N); for n 2:N xn x(n); if n - k 1 xn xn - x(n - k); end if n - k - l 1 xn xn - x(n - k - l); end if n - 2*k - l 1 xn xn x(n - 2*k - l); end d(n) xn; p(n) p(n-1) d(n); y(n) y(n-1) p(n); end end跑完之后直接 plot(x) 和 plot(y) 对比你可以清楚看到指数脉冲变成了标准的梯形脉冲。这个阶段建议多试几组 k 和 l 参数。比如 k 取 50l 取 50然后 k 取 150l 取 100观察输出梯形的高矮胖瘦变化。实际工程里k 和 l 的选择直接影响信噪比和堆积容忍能力参数扫描这一步值得做。这里有一个细节因为没有做基线扣除累加器 p 和 y 会把直流偏置累积得特别大导致输出波形整体抬高甚至饱和。核脉冲处理里的基线漂移是个老大难问题。MATLAB仿真阶段可以直接用 detrend 或者减去平均值但写Verilog的时候就必须在电路里显式处理基线了这个后面专门讲。2.3 定点化算法能跑和硬件能跑是两回事MATLAB里用的double浮点,搬进FPGA不可能直接这么搞。FPGA里做浮点不是不行但资源开销大、时序不好收敛对50MHz采样率来说完全是杀鸡用牛刀。更常见的做法是把整个递推式全部转成定点整数运算。位宽规划是这个环节最核心的决策。先算一下数据范围。假设ADC输出是16位无符号数最大满量程是65535实际峰值大约在30000左右。递推式里两次累加p 的数值量和输入幅度、k 值直接相关。一个满幅脉冲经过梯形成形输出幅度大约正比于 A*k如果 A30000k128理论峰值就是3840000需要22位左右。考虑到累加器的历史效应和脉冲堆积直接给32位中间位宽是最稳妥的。我的建议是输入16位中间 d 用 32 位有符号数p 和 y 都用 32 位有符号数最终输出按需右移截取到16位或18位。右移多少位取决于 k/l 参数和上位机期望的幅度范围这个参数在工程里可以做成寄存器上电后再标定调优。为了让Verilog的定点行为在MATLAB里先被验证我一般会额外写一个定点版本仿真函数把所有运算都用int32来模拟再和浮点结果对比看误差。% 定点模拟假设输入已经量化为整数 function yd trapezoid_fixed(x_int, k, l) N length(x_int); acc int32(0); p int32(0); yd zeros(1, N, int32); for n 2:N dval int32(x_int(n)); if n - k 1 dval dval - int32(x_int(n - k)); end if n - k - l 1 dval dval - int32(x_int(n - k - l)); end if n - 2*k - l 1 dval dval int32(x_int(n - 2*k - l)); end p p dval; acc acc p; yd(n) acc; end end把浮点结果和定点结果放在同一张图里对比如果曲线基本重合说明这个位宽方案可行。实测下来32位中间累加器对于16位输入、k 在250以内的情况误差可以控制在0.1%以内完全满足能谱测量的需要。3. Verilog实现从浮点到硬件的关键设计3.1 顶层架构和模块划分MATLAB验证完算法之后开始思考Verilog架构。50MHz采样率不算高但设计得好不好直接决定了资源占用和后续扩展空间。顶层模块的功能很明确每个时钟周期进来一个ADC采样点 x[n]经过三个延迟项和两级积分器输出一个 y[n]同时输出一个 dout_valid 标志告诉后级“数据有效”。此外k 和 l 两个参数通过寄存器接口动态配置方便上板后调参。模块划分我比较习惯这样分输入同步与签转换把ADC的16位无符号数转成有符号数顺便抑制亚稳态延迟线模块提供 x[n-k]、x[n-k-l]、x[n-2k-l] 三个历史数据差分计算单元计算 d[n]两级累加器计算 p[n] 和 y[n]输出对齐模块把 dout 和 dout_valid 打拍对齐为什么要把输出对齐单独拎出来因为延迟线读历史数据可能需要多个时钟周期dout 算出来的时候dout_valid 必须严格同步拉高。一旦valid和数据错位后面的FIFO或者ILA抓波形时就会看到输出波形整体错位一个周期排查起来非常吃亏。顶层接口定义大致是这样module trap_shaper #( parameter DIN_WIDTH 16, parameter ACC_WIDTH 32, parameter MAX_DELAY 4096, parameter ADDR_WIDTH 12 )( input wire clk, input wire rst_n, input wire signed [DIN_WIDTH-1:0] din, input wire din_valid, input wire [ADDR_WIDTH-1:0] k_para, input wire [ADDR_WIDTH-1:0] l_para, output reg signed [ACC_WIDTH-1:0] dout, output reg dout_valid );3.2 延迟线实现三种方案和取舍梯形成形算法的核心资源消耗全在三个历史数据的延迟访问上。 x[n-k]、x[n-k-l]、x[n-2k-l] 这三个点本质上是要把输入数据延迟 k、kl、2kl 拍。k 和 l 的范围通常在几十到几千之间这个深度范围让“延迟线怎么实现”成了一个有讲究的问题。方案一纯移位寄存器链。把数据打到一个深度为 MAX_DELAY 的移位寄存器里需要读哪个延迟点就从对应的寄存器输出引出来。这在参数小的时候最简单比如最大延迟100拍以内直接reg [DIN_WIDTH-1:0] shift_reg [0:MAX_DELAY]就完事了。但问题是如果 MAX_DELAY 做到4096每个通道要4096个16位寄存器那就是65K个FF资源直接爆炸布线也吃紧。方案二用Xilinx的SRL16/SRLC32E原语。SRL可以把一个LUT变成32位移位寄存器资源利用率比FF链高得多。但缺点是只能做固定延迟动态可调的延迟不太方便而且一个SRL只有一个抽头三个延迟点要三路独立延迟链或者是额外的read pointer逻辑复杂度上来了。方案三用Block RAM做环形缓冲。数据一直往RAM里写写地址就是当前时间戳要访问 n-k 时刻的数据只需要计算wr_addr - k作为读地址从同一个RAM里读出来。这种方式最灵活k/l 可以动态配置资源消耗也稳定深度增大只影响地址位宽不影响资源量也就是常用的延迟线RAM方案。我实际工程里推荐方案三。BRAM是FPGA里非常充裕的资源单块36K BRAM可以配置成16K x 16位如果你最大延迟只用到4096那么一块BRAM可以同时服务好几个通道。唯一的问题是BRAM的端口数量有限。一个简单双端口RAM只有一个读端口和一个写端口同一时钟周期只能读一个历史点但我们需要三个延迟点的数据这就触发了下一个关键设计。3.3 流水线计算单元设计三路历史数据不能在同一周期全部读出最直接的办法是把计算拆成三个周期分时读出三个历史值这正好和50MHz的低速率匹配。流水线安排如下周期1读出 x[n-k]计算 temp1 x[n] - x[n-k]周期2读出 x[n-k-l]计算 temp2 temp1 - x[n-k-l]周期3读出 x[n-2k-l]计算 temp3 temp2 x[n-2k-l]同时更新 p 和 y这么安排逻辑链很浅时序非常友好。代价是输出比输入延迟了3个时钟周期对核脉冲处理来说完全可以接受。而且流水线本来就有延迟只要valid信号跟着打拍就行。下面给出核心代码// 环形缓冲写端口 reg [ADDR_WIDTH-1:0] wr_addr; reg [ADDR_WIDTH-1:0] rd_addr; always (posedge clk or negedge rst_n) begin if (!rst_n) begin wr_addr {ADDR_WIDTH{1b0}}; end else if (din_valid) begin wr_addr wr_addr 1b1; mem[wr_addr] din; end end // 三周期流水线读取 always (posedge clk or negedge rst_n) begin if (!rst_n) begin rd_addr {ADDR_WIDTH{1b0}}; stage 2d0; end else if (din_valid) begin case (stage) 2d0: begin rd_addr wr_addr - k_para; stage 2d1; end 2d1: begin rd_addr wr_addr - k_para - l_para; stage 2d2; end 2d2: begin rd_addr wr_addr - {1b0, k_para} - {1b0, k_para} - l_para; stage 2d0; end endcase end end // 中间变量和两级积分 wire signed [DIN_WIDTH-1:0] delay_sample; reg signed [ACC_WIDTH-1:0] t1, t2, dval; reg signed [ACC_WIDTH-1:0] p_acc, y_acc; always (posedge clk or negedge rst_n) begin if (!rst_n) begin t1 {ACC_WIDTH{1b0}}; t2 {ACC_WIDTH{1b0}}; dval {ACC_WIDTH{1b0}}; p_acc {ACC_WIDTH{1b0}}; y_acc {ACC_WIDTH{1b0}}; end else if (din_valid) begin case (stage) 2d0: t1 din - delay_sample; 2d1: t2 t1 - delay_sample; 2d2: begin dval t2 delay_sample; p_acc p_acc dval; y_acc y_acc p_acc; end endcase end end代码里有个细节要注意wr_addr - k_para里的减法在Verilog里如果两边的位宽都是ADDR_WIDTH那么WRAP之后自动取模正好满足环形缓冲的地址回绕需求不需要额外判断。dout 和 dout_valid 的对齐可以从 din_valid 打出3拍延迟。更稳妥的做法是直接用一个2位计数器每次 din_valid 有效时从0数到2数到2时拉高 dout_valid。这样无论中间怎么改只要 din_valid 的节奏稳定输出就不会错位。还有一点k/l参数不要在数据流进行中随意修改。因为环形缓冲里的历史数据是连续流动的中途改变读地址偏移会导致输出波形瞬间错乱这种事我只在调试时干过一次后果就是波形像被鬼畜了一样来回乱跳。正确做法是先把 din_valid 拉低让流水线排空等新参数写入后再重新开始。3.4 仿真验证与MATLAB结果闭环Verilog写完不能直接上板先仿真。仿真最关键的一步是把MATLAB里生成的测试数据喂给testbench再把输出的结果拿回MATLAB对比。这样形成了一个完整的闭环算法来源和验证标准都在MATLAB硬件只是这个算法的忠实搬运工。testbench的核心流程是这样的initial begin $readmemh(input_data.hex, mem); for (i 0; i DATA_NUM; i i 1) begin (posedge clk); din mem[i]; din_valid 1b1; end din_valid 1b0; endMATLAB那边负责把仿真数据写成十六进制文本一行一个数fid fopen(input_data.hex, w); fprintf(fid, %04X\n, int16(x_int)); fclose(fid);仿真完成后同样把 dout 和 dout_valid 写成文件带回MATLAB里和定点仿真结果比较。我当时比对的时候发现两条曲线几乎完全重合差异只有几个最低位说明Verilog对算法的还原度非常高这个信号非常重要可以放心上板了。这里有个仿真技巧为了充分验证流水线testbench里不要每个周期都拉高din_valid要穿插一些空闲周期模拟真实ADC数据流不连续的情况。很多新人写testbench默认每周期都有效结果上板才发现数据流一停输出就乱了对齐逻辑浪费一整天。4. 50MHz采样率部署上板实战与问题排查4.1 资源与时序评估先把读者最关心的资源数据给出来。以Xilinx Artix-7 XC7A35T为例在50MHz采样率、单通道、16位输入、最大延迟4096的配置下整个梯形成形模块的资源消耗大概是这样资源类型消耗量说明Block RAM1 块 36K环形缓冲16K x 16位配置DSP48E10纯加法结构用不到DSPLUT约150个状态机地址计算逻辑FF约200个流水线寄存器与参数缓存时钟频率上限远超150MHz50MHz下预留大量裕量这个资源消耗在一个35T的芯片里只占不到2%。也就是说剩下的大量资源可以用来做多通道扩展、基线估计、堆积修正、能谱统计等更高层的事情。我后来实际做过8通道的版本BRAM共享之后资源增量很小大部分开销都在通道控制逻辑上。时序方面最长的组合逻辑路径通常在地址减法计算上而不是累加器。因为累加器 p_acc 和 y_acc 之间有直接的反馈依赖看起来像是关键路径的主角但你把两级积分拆到一个时钟周期里执行的时候实际上就是从 p_acc 到 y_acc 到 dout 的一串加法链中间必须插一级寄存器把 p_acc 先锁存否则50MHz能过往上超频到100MHz就会开始时序违例。有一个小优化非常管用把dout的两位右移放到输出级做。也就是说FPGA内部保留32位的全精度数据输出给后级时用dout[ACC_WIDTH-1-2 -: 16]截取高16位。这样既保证精度又让后级的FIFO位宽不用做太宽。4.2 接入真实数据流ADC和基线处理上板以后我先把FPGA内部的测试信号发生器打开生成一个正弦波或者可调频率的脉冲序列直接灌进梯形成形模块。确认在纯数字环路里输出波形正常再接真实ADC的数据。ADC数据的接入有几个细节会影响梯形成形效果。第一个是数据同步。如果ADC是并行CMOS接口直接打两拍同步就行如果是LVDS串行接口需要用ISERDES和Bitslip做串并转换这块如果没做过的话建议先用并行接口板子点亮把精力集中在成形算法上别一开始就和LVDS死磕。第二个是基线问题。ADC输入通常有一个直流偏置比如1.25V对应数字码的中间值8000这个偏置在递推式里会被两级积分器不断累加最后导致输出基线严重偏移甚至把有用的脉冲推到负半轴。所以我在FPGA里加了一级基线估计用一个一阶IIR低通滤波器跟踪输入信号的慢变均值然后把每个输入点减去这个估计值。// 基线估计与扣除alpha约等于2^-10 baseline baseline (din - baseline) 10; din_bs din - baseline;这种做法在信号大部分时间是基线、只有偶尔脉冲到达的场景下效果很好。注意alpha系数不能太大否则会把脉冲的能量也一起吃掉太小又跟不上慢漂移一般取2^-8到2^-12之间比较合适。这个参数也可以做成可配置的上板后看波形再调。4.3 常见问题速查表这个环节直接上干货把我实际用ILA抓波形时踩过的坑全部列出来一条一条对照查问题就行。异常现象可能原因解决方法输出梯形顶部有锯齿台阶流水线dout_valid没对齐检查valid打拍是否与数据严格同步梯形上升沿出现反向过冲边界条件没处理好起始段有负序号数据检查延迟线初始化等待一定点数后再输出平顶不平整中间下凹基线扣除过度或信号本身基线漂移调小基线估计alpha系数输出幅度饱和截平累加器位宽不够溢出回绕加宽ACC_WIDTH或增加右移位数大量脉冲堆积时幅度偏小k/l参数太大成形窗覆盖到相邻脉冲减小k/l或加堆积拒绝逻辑修改k/l参数后波形乱跳流水线未排空就切换参数先拉低din_valid等待2kl个周期再改参数长期运行后基线缓慢上移低频噪声或直流偏置累积增加基线估计环路或高通滤波级上板后输出全零din_valid信号未正确拉高用ILA抓din_valid查看检查ADC接口状态这些坑里我觉得最容易被人忽视的是第一个valid对齐。仿真阶段如果输入数据每周期都有效valid对齐错误不容易暴露但一旦数据流有间隙或者你加了FIFO缓冲valid错位的问题就会彻底爆发。调试这个问题的秘诀是用ILA同时抓 din、din_valid、dout、dout_valid 四个信号然后在MATLAB里数一下从din到dout的延迟周期数再检查testbench里预期的是否一致。4.4 参数上板调试与结果验证参数调试阶段我习惯把k、l、右移位数、基线alpha都做成VIOVirtual I/O寄存器这样可以在线修改不用每次改代码重新综合。上板后先用信号发生器注入一个固定幅度的方波脉冲让梯形成形输出一个标准的梯形然后用ILA把波形抓下来人工观察平顶的平坦度。如果手头有实际的探测器信号或者脉冲信号源可以把输出结果通过UART发回上位机在上位机里统计脉冲峰值画能谱直方图。这个环节做完整个系统才算是真正打通了。我测过一组参数k取60l取80频率50MHz信号源发出幅度1V、宽度500ns的脉冲输出梯形上升沿和平顶都非常干净平顶波动在2个LSB以内边缘没有明显过冲。这个结果用来做后续的能谱测量分辨率表现是令人满意的。5. 收尾小经验写到这儿整个项目的技术链路基本完整了。根据我个人的调参经验再说几句可以延伸优化的思路。k和l的参数选择直接影响整个系统的分辨率上限。k可以理解为信噪比调节旋钮k越大累积的采样点越多平均效果越好噪声压低但代价是成形窗变宽高计数率时堆积概率上升。l的主要作用是抵抗弹道亏损和采样时间抖动平顶太窄对相位抖动敏感太宽又牺牲通过率。一般从脉冲的上升沿时间做起点比如50MHz下脉冲上升沿大约2us那就是100个点k选80到120左右l选50到100再根据实测能谱的峰位漂移和FWHM来微调。如果后续要做更高性能的系统可以考虑在梯形成形后面加一级数字基线恢复或者在Zynq平台把成形逻辑放PL、能谱统计放PS利用ARM核跑算法自适应的参数寻优。多通道场景下还可以把各通道的RAM做成共享缓冲池进一步压缩资源。最后分享一个工程习惯FPGA项目里算法仿真通过只是长征第一步真正的调试战场在上板环节而且绝大多数问题都出在时序控制、数据有效标志、位宽边界这类“算法之外”的细节上。所以从第一天开始就把valid、复位、边界这些辅助信号当成核心代码一样打样处理后面能省非常多的时间。这套梯形成形项目做完之后我回头看最值钱的收获不是会了这个算法而是学会了怎么把一个DSP算法从头到尾在FPGA里跑通这套方法论用在其他算法上依然有效。
返回列表