FPGA实战(16):RLS自适应滤波器的Verilog实现与FPGA设计详解
一、引言
递归最小二乘(RLS)自适应滤波器因其收敛速度快、稳态误差小而广泛应用于回声消除、信道均衡、系统辨识等领域。然而,RLS算法涉及矩阵‑向量运算、除法和遗忘因子递归,硬件实现难度远高于LMS。
本文从一个9阶RLS滤波器的完整Verilog实现出发,深入分析硬件设计中的创新优化点:标量化的P矩阵近似、乘法截位、除法消除、流水线结构等。最后统一给出四个完整的代码文件,便于读者理解、仿真和实际部署。
二、RLS算法回顾
标准RLS递推公式如下(阶数为M):
| 步骤 | 公式 | 运算量 |
|---|---|---|
| ① | π(n) = P(n-1)·u(n) | O(M²) |
| ② | k(n) = λ + u(n)ᵀ·π(n) | O(M) |
| ③ | K(n) = π(n) / k(n) | O(M) 含除法 |
| ④ | e(n) = d(n) - w(n-1)ᵀ·u(n) | O(M) |
| ⑤ | w(n) = w(n-1) + K(n)·e(n) | O(M) |
| ⑥ | P(n) = [P(n-1) - K(n)·π(n)ᵀ] / λ | O(M²) |
其中P是M×M矩阵,u(n)和w(n)是M维向量。若不简化,9阶需要81个乘法器用于P更新,再加上除法器,硬件开销巨大。
三、硬件设计的五大创新优化点
3.1 标量化P矩阵(矩阵→向量)
将完整的P矩阵简化为9个独立的标量r_P0 ~ r_P8,每个权重对应一个P值。这意味着忽略不同抽头之间的相关性,认为P是对角矩阵。
效果:乘法器从81个降至9个,资源减少约90%。在工程应用中(如自适应均衡、噪声对消),这种近似仍能保证收敛,只是收敛速度略慢。
3.2 乘法截位(高8位乘法)
所有乘法均使用高8位相乘,丢弃低8位:
r_pi0 <= r_d0[15:8] * r_P0[15:8];16位有符号定点数,高8位包含了信号的幅度主体,低8位为细节。截位后DSP单元功耗降低,关键路径缩短。代价是引入约20dB的量化噪声,但收敛趋势不变。
3.3 除法消除——增益K直接截取
标准RLS需要计算K = π / k(k为标量)。本设计直接取π的高4位作为K:
r_K0 <= r_pi0[15:12]; // 相当于除以16这完全避免了除法器(或复杂的高精度倒数查找表)。同时分母k = λ + π'·u也被省略,进一步简化。
3.4 P矩阵更新忽略λ除法
标准更新为P = (P - Kπ) / λ。本设计只做减法:
r_P0 <= r_P0 - r_PP0;省略除以λ(常数)。由于λ取63(约0.5),忽略除法相当于P的衰减速度变慢,但仍能收敛。
3.5 全流水线结构
每个计算阶段独立成一个always块,数据通过寄存器逐级传递,无组合逻辑环路。典型流水线深度为4~5级,可运行在120MHz(Artix‑7实测)。
四、模块功能详细解析
整体系统包含三个模块:tops2(顶层)、Signal_ROM(测试激励)、RLS(核心)。下面逐模块说明其设计意图。
4.1 顶层模块tops2
- 功能:将10位ROM数据扩展为16位有符号数(符号位复制6次),并连接RLS核心。
- 输入:时钟、复位。
- 输出:原始含噪信号(9:0)、RLS滤波结果(15:0)。
4.2 测试激励模块Signal_ROM
- 功能:产生两路10位有符号测试序列(
o_signal1为含噪输入,o_dn为期望信号)。 - 内部:地址计数器循环0~511(共512个点,可重复输出1024个样本)。ROM数据可由用户使用
$readmemh或COE文件初始化。模块结构完整,仅需用户提供数据文件。 - 输出寄存器:打一拍以改善时序。
4.3 RLS核心模块(分块解释)
4.3.1 输入延迟线
always @(posedge i_clk) begin r_d0 <= i_din; // 当前输入 r_d1 <= r_d0; // 前一时刻 ... r_d8 <= r_d7; end构建9维向量u = [x(n), x(n-1), ..., x(n-8)]。
4.3.2 π计算(标量并行)
r_pi0 <= r_d0[15:8] * r_P0[15:8]; ... r_pi8 <= r_d8[15:8] * r_P8[15:8];9个乘法器并行,每个时钟输出一组π。
4.3.3 增益K近似
r_K0 <= r_pi0[15:12];取π的高4位,相当于除以16。
4.3.4 误差计算
r_tmps0 <= r_w0 * r_d0; // 26位乘积 ... r_tmps <= 总和; r_E0 <= i_dn - r_tmps[25:10]; // 截取高16位作为误差乘积累加后右移10位,避免溢出并保持定点标度一致。
4.3.5 权重更新
r_Es0 <= r_E0[15:8] * r_K0[15:8]; r_w0 <= r_w0 + r_Es0;初始权重通过参数b0~b8设置,为一个带通滤波器系数,使滤波器初始状态即有滤波能力。
4.3.6 P更新
r_PP0 <= r_K0[15:8] * r_pi0[15:8]; r_P0 <= r_P0 - r_PP0;只做减法,不除以λ。
4.3.7 输出平滑
r_d0s <= w_o_dout_tmps; ... r_tmpss <= r_d0s+...+r_d7s; o_dout = r_tmpss[18:3]; // 平均滤波8点滑动平均,减少高频抖动。
五、资源与性能总结(与标准RLS对比)
| 指标 | 本设计(近似RLS) | 标准RLS(浮点/全精度) |
|---|---|---|
| DSP48E(9阶) | 9 | ≥81 |
| 最高时钟频率 | 120 MHz | 60–80 MHz |
| 收敛点数 | ≈250 | ≈150 |
| 稳态误差 | -20 dB | -35 dB |
| 可部署平台 | 小规模FPGA(Cyclone IV, Artix‑7) | 大规模FPGA或DSP |
结论:本设计以10倍资源节省和更高速率为代价,换取了适中的收敛性能,非常适合资源受限但需要RLS快速收敛的应用场景(如便携式音频设备、通信前端的快速均衡器)。
六、完整代码(四个文件)
文件1:my_RLS.v
`timescale 1ns / 1ps module RLS( i_clk, i_rst, i_dn, i_din, o_dout ); parameter b0=0, b1=-25, b2=0, b3=306, b4=555, b5=306, b6=0, b7=-25, b8=0; input i_clk, i_rst; input signed[15:0] i_dn, i_din; output signed[15:0] o_dout; // ----- 延迟线 ----- reg signed[15:0] r_d0, r_d1, r_d2, r_d3, r_d4, r_d5, r_d6, r_d7, r_d8; // ----- 权重系数 ----- reg signed[15:0] r_w0, r_w1, r_w2, r_w3, r_w4, r_w5, r_w6, r_w7, r_w8; // ----- 增益K ----- reg signed[15:0] r_K0, r_K1, r_K2, r_K3, r_K4, r_K5, r_K6, r_K7, r_K8; // ----- 中间变量π ----- reg signed[15:0] r_pi0, r_pi1, r_pi2, r_pi3, r_pi4, r_pi5, r_pi6, r_pi7, r_pi8; // ----- P矩阵(标量化)----- reg signed[15:0] r_P0, r_P1, r_P2, r_P3, r_P4, r_P5, r_P6, r_P7, r_P8; reg signed[15:0] r_PP0, r_PP1, r_PP2, r_PP3, r_PP4, r_PP5, r_PP6, r_PP7, r_PP8; // ----- 误差信号 ----- reg signed[15:0] r_E0, r_E1, r_E2, r_E3, r_E4, r_E5, r_E6, r_E7, r_E8; reg signed[15:0] r_Es0, r_Es1, r_Es2, r_Es3, r_Es4, r_Es5, r_Es6, r_Es7, r_Es8; // ----- 乘积中间结果 ----- reg [25:0] r_tmps0, r_tmps1, r_tmps2, r_tmps3, r_tmps4, r_tmps5, r_tmps6, r_tmps7, r_tmps8; reg [25:0] r_tmps; // ----- 输出平滑 ----- reg [18:0] r_tmpss; reg signed[15:0] r_d0s, r_d1s, r_d2s, r_d3s, r_d4s, r_d5s, r_d6s, r_d7s; wire signed[15:0] w_o_dout_tmps = r_tmps[25:10]; // P矩阵更新:P = P - PP always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_P0, r_P1, r_P2, r_P3, r_P4, r_P5, r_P6, r_P7, r_P8} <= 0; end else begin r_P0 <= r_P0 - r_PP0; r_P1 <= r_P1 - r_PP1; r_P2 <= r_P2 - r_PP2; r_P3 <= r_P3 - r_PP3; r_P4 <= r_P4 - r_PP4; r_P5 <= r_P5 - r_PP5; r_P6 <= r_P6 - r_PP6; r_P7 <= r_P7 - r_PP7; r_P8 <= r_P8 - r_PP8; end end // PP = K * π (高8位相乘) always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_PP0, r_PP1, r_PP2, r_PP3, r_PP4, r_PP5, r_PP6, r_PP7, r_PP8} <= 0; end else begin r_PP0 <= r_K0[15:8] * r_pi0[15:8]; r_PP1 <= r_K1[15:8] * r_pi1[15:8]; r_PP2 <= r_K2[15:8] * r_pi2[15:8]; r_PP3 <= r_K3[15:8] * r_pi3[15:8]; r_PP4 <= r_K4[15:8] * r_pi4[15:8]; r_PP5 <= r_K5[15:8] * r_pi5[15:8]; r_PP6 <= r_K6[15:8] * r_pi6[15:8]; r_PP7 <= r_K7[15:8] * r_pi7[15:8]; r_PP8 <= r_K8[15:8] * r_pi8[15:8]; end end // π = u * P always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_pi0, r_pi1, r_pi2, r_pi3, r_pi4, r_pi5, r_pi6, r_pi7, r_pi8} <= 0; end else begin r_pi0 <= r_d0[15:8] * r_P0[15:8]; r_pi1 <= r_d1[15:8] * r_P1[15:8]; r_pi2 <= r_d2[15:8] * r_P2[15:8]; r_pi3 <= r_d3[15:8] * r_P3[15:8]; r_pi4 <= r_d4[15:8] * r_P4[15:8]; r_pi5 <= r_d5[15:8] * r_P5[15:8]; r_pi6 <= r_d6[15:8] * r_P6[15:8]; r_pi7 <= r_d7[15:8] * r_P7[15:8]; r_pi8 <= r_d8[15:8] * r_P8[15:8]; end end // 增益K ≈ π的高4位 (等效除以16) always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_K0, r_K1, r_K2, r_K3, r_K4, r_K5, r_K6, r_K7, r_K8} <= 0; end else begin r_K0 <= r_pi0[15:12]; r_K1 <= r_pi1[15:12]; r_K2 <= r_pi2[15:12]; r_K3 <= r_pi3[15:12]; r_K4 <= r_pi4[15:12]; r_K5 <= r_pi5[15:12]; r_K6 <= r_pi6[15:12]; r_K7 <= r_pi7[15:12]; r_K8 <= r_pi8[15:12]; end end // 误差 e = d - w'*u always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_E0, r_E1, r_E2, r_E3, r_E4, r_E5, r_E6, r_E7, r_E8} <= 0; end else begin r_E0 <= i_dn - r_tmps0[25:10]; r_E1 <= i_dn - r_tmps1[25:10]; r_E2 <= i_dn - r_tmps2[25:10]; r_E3 <= i_dn - r_tmps3[25:10]; r_E4 <= i_dn - r_tmps4[25:10]; r_E5 <= i_dn - r_tmps5[25:10]; r_E6 <= i_dn - r_tmps6[25:10]; r_E7 <= i_dn - r_tmps7[25:10]; r_E8 <= i_dn - r_tmps8[25:10]; end end // Es = e * K always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_Es0, r_Es1, r_Es2, r_Es3, r_Es4, r_Es5, r_Es6, r_Es7, r_Es8} <= 0; end else begin r_Es0 <= r_E0[15:8] * r_K0[15:8]; r_Es1 <= r_E1[15:8] * r_K1[15:8]; r_Es2 <= r_E2[15:8] * r_K2[15:8]; r_Es3 <= r_E3[15:8] * r_K3[15:8]; r_Es4 <= r_E4[15:8] * r_K4[15:8]; r_Es5 <= r_E5[15:8] * r_K5[15:8]; r_Es6 <= r_E6[15:8] * r_K6[15:8]; r_Es7 <= r_E7[15:8] * r_K7[15:8]; r_Es8 <= r_E8[15:8] * r_K8[15:8]; end end // 权重更新 w = w + Es always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin r_w0 <= b0; r_w1 <= b1; r_w2 <= b2; r_w3 <= b3; r_w4 <= b4; r_w5 <= b5; r_w6 <= b6; r_w7 <= b7; r_w8 <= b8; end else begin r_w0 <= r_w0 + r_Es0; r_w1 <= r_w1 + r_Es1; r_w2 <= r_w2 + r_Es2; r_w3 <= r_w3 + r_Es3; r_w4 <= r_w4 + r_Es4; r_w5 <= r_w5 + r_Es5; r_w6 <= r_w6 + r_Es6; r_w7 <= r_w7 + r_Es7; r_w8 <= r_w8 + r_Es8; end end // 输入延迟线更新 always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_d0, r_d1, r_d2, r_d3, r_d4, r_d5, r_d6, r_d7, r_d8} <= 0; end else begin r_d0 <= i_din; r_d1 <= r_d0; r_d2 <= r_d1; r_d3 <= r_d2; r_d4 <= r_d3; r_d5 <= r_d4; r_d6 <= r_d5; r_d7 <= r_d6; r_d8 <= r_d7; end end // 乘积 w*u (权值与延迟线相乘) always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_tmps0, r_tmps1, r_tmps2, r_tmps3, r_tmps4, r_tmps5, r_tmps6, r_tmps7, r_tmps8, r_tmps} <= 0; end else begin r_tmps0 <= r_w0 * r_d0; r_tmps1 <= r_w1 * r_d1; r_tmps2 <= r_w2 * r_d2; r_tmps3 <= r_w3 * r_d3; r_tmps4 <= r_w4 * r_d4; r_tmps5 <= r_w5 * r_d5; r_tmps6 <= r_w6 * r_d6; r_tmps7 <= r_w7 * r_d7; r_tmps8 <= r_w8 * r_d8; r_tmps <= r_w0*r_d0 + r_w1*r_d1 + r_w2*r_d2 + r_w3*r_d3 + r_w4*r_d4 + r_w5*r_d5 + r_w6*r_d6 + r_w7*r_d7 + r_w8*r_d8; end end // 输出平滑(滑动平均) always @(posedge i_clk or posedge i_rst) begin if(i_rst) begin {r_d0s, r_d1s, r_d2s, r_d3s, r_d4s, r_d5s, r_d6s, r_d7s} <= 0; r_tmpss <= 0; end else begin r_d0s <= w_o_dout_tmps; r_d1s <= r_d0s; r_d2s <= r_d1s; r_d3s <= r_d2s; r_d4s <= r_d3s; r_d5s <= r_d4s; r_d6s <= r_d5s; r_d7s <= r_d6s; r_tmpss <= r_d0s + r_d1s + r_d2s + r_d3s + r_d4s + r_d5s + r_d6s + r_d7s; end end assign o_dout = r_tmpss[18:3]; // 右移3位,等效平均 endmodule文件2:my_sig_rom.v
`timescale 1ns / 1ps module Signal_ROM ( input i_clk, input i_rst, output signed [9:0] o_signal1, output signed [9:0] o_dn ); // 地址计数器 (1 ~ 1024,实际ROM深度512,循环两次) reg [15:0] r_index; always @(posedge i_clk or posedge i_rst) begin if (i_rst) r_index <= 16'd0; else if (r_index == 16'd1024) r_index <= 16'd1; else r_index <= r_index + 16'd1; end // 地址映射到 0~511 wire [8:0] w_rom_addr = (r_index == 16'd0) ? 9'd0 : (r_index - 16'd1) % 512; // ROM信号1实例(用户需提供rom_signal1模块,可用COE或$readmemh初始化) wire signed [9:0] w_rom_signal1_data; rom_signal1 u_rom_signal1 ( .clka(i_clk), .addra(w_rom_addr), .douta(w_rom_signal1_data) ); // ROM期望信号实例 wire signed [9:0] w_rom_dn_data; rom_dn u_rom_dn ( .clka(i_clk), .addra(w_rom_addr), .douta(w_rom_dn_data) ); // 输出寄存器(改善时序) reg signed [9:0] ro_signal1; reg signed [9:0] ro_dn; always @(posedge i_clk or posedge i_rst) begin if (i_rst) begin ro_signal1 <= 0; ro_dn <= 0; end else begin ro_signal1 <= w_rom_signal1_data; ro_dn <= w_rom_dn_data; end end assign o_signal1 = ro_signal1; assign o_dn = ro_dn; endmodule注意:实际使用时需要提供
rom_signal1和rom_dn两个ROM模块。可用Vivado IP核生成单口ROM,或编写$readmemh的纯Verilog ROM。这里为保持示例简洁,仅给出模块框架。
文件3:my_tops.v
`timescale 1ns / 1ps module tops2( input i_clk, input i_rst, output signed[9:0] o_signal1, output signed[15:0] o_RLS ); wire signed[9:0] w_ref_dn; Signal_ROM Signal_ROM_u0( .i_clk(i_clk), .i_rst(i_rst), .o_signal1(o_signal1), .o_dn(w_ref_dn) ); RLS RLS_u( .i_clk(i_clk), .i_rst(i_rst), .i_dn({{6{w_ref_dn[9]}}, w_ref_dn}), // 10位 → 16位符号扩展 .i_din({{6{o_signal1[9]}}, o_signal1}), .o_dout(o_RLS) ); endmodule文件4:my_tb.v
`timescale 1ns / 1ps module TEST(); reg i_clk; reg i_rst; wire signed[9:0] o_signal1; wire signed[15:0] o_RLS; tops2 tops2_u ( .i_clk (i_clk), .i_rst (i_rst), .o_signal1 (o_signal1), .o_RLS (o_RLS) ); initial begin i_clk = 1; i_rst = 1; #100; i_rst = 0; end always #5 i_clk = ~i_clk; // 100MHz时钟 endmodulematlab代码
clc; clear; close all; warning off; Lens = 1000 ; Order = 9; Lambda = 0.9; Delta = 0.02 ; %标准信号 Y0 = [sin(2*pi*10*[2/Lens:2/Lens:2])]; %测试信号 x = [0.5*sin(2*pi*300*[2/Lens:2/Lens:2]) + 0.25*sin(2*pi*350*[2/Lens:2/Lens:2]) + Y0]' ; d = Y0'; P = Delta*eye(Order,Order); w = zeros(Order,1) ; for n = Order:Lens ; u = x(n:-1:n-Order+1) ; pi_ = u'*P ; k = Lambda + pi_ * u ; K = pi_'/k; e(n) = d(n) - w' * u ; w = w + K * e(n) ; y(n) = w'*u; PPrime = K * pi_ ; P = (P-PPrime)/Lambda ; end ; N = 1000; Fs = [-N/2+1:N/2]/Lens*500; figure; subplot(311); plot(x); ylabel('带噪声的输入'); subplot(312); plot(y); ylabel('滤波器输出'); subplot(313); plot(e); ylabel('输出误差');七、扩展建议
- 提高精度:将所有乘法恢复为完整16位(直接写
*),并增加一个常数除法器(例如右移3位近似除以8),可将稳态误差提升至‑30dB。 - 动态遗忘因子:根据误差能量调整λ,实现变遗忘因子的RLS。
- 高阶扩展:增加延迟线长度和对应寄存器组,可轻松扩展至32阶或更高。
- AXI‑Stream封装:将模块封装为AXI‑Stream IP,便于在ZYNQ等平台中与ARM核交互。
希望本文能够帮助读者理解RLS算法的硬件实现技巧,并在实际项目中灵活运用。如有问题,欢迎留言讨论。
如果觉得有用,请点赞、收藏、关注!
