本篇博客源于b站上的一个视频，视频里up主不仅讲述了FPGA上实现LMS算法，还讲述了在实现一个算法时的整体的思路，很值得学习，想了解这种思路的可以去看一下视频，这里主要将算法与算法实现过程总结为笔记，视频链接如下：
verilog实现LMS自适应滤波算法（FPGA）
另外视频中用到的文献如下：

刘雄飞,高金定,齐海兵. LMS自适应滤波器FPGA实现的新方法[J]. 压电与声光,2007,29(1):87-89. DOI:10.3969/j.issn.1004-2474.2007.01.029.

LMS算法

LMS:Least Mean Square,最小均方
自适应滤波器的模型如下图所示：

输入信号矢量：$ \pmb{x}(n)=[x(n) x(n-1) ... x(n-L)]^T $

输出y为：$ y(n)=\sum_{k=0}^{L}w_k(n)x(n-k) $

权系数可以组成权系数矢量：$ \pmb{w}(n)=[w_0(n) w_1(n) ... w_L(n)]^T $

因此输出可以表示为：$ y(n)= \pmb{w}^T (n)\pmb{x}(n) $

误差信号：$ e(n)=d(n)-\pmb{w}^T (n)\pmb{x}(n) $

定义均方误差：$ J(n)=E[ e^2 (n)]=E[ d^2 (n)]+\pmb{w}^T (n) E[ \pmb{x}(n)\pmb{x}^T(n) ]\pmb{w}(n)-2E[ d(n)\pmb{x}^T(n) ]\pmb{w}(n) $

因此定义输入信号的自相关矩阵为：$ \pmb{R}=E[ \pmb{x}(n)\pmb{x}^T(n) ] $

定义期望信号与输入信号的互相关矩阵为：$ \pmb{P}=E[ d(n)\pmb{x}^T(n) ] $

因此均方误差可表示为：$ J(n)=E[ d^2 (n)]+\pmb{w}^T \pmb{R}\pmb{w}-2\pmb{P}^T \pmb{w} $

梯度：$ \nabla=\frac{\partial J}{\partial w}=2\pmb{R}\pmb{w}-2\pmb{P} $

令梯度等于零，即可求得最小均方误差对应的维纳解：$ \pmb{w}^* =\pmb{R}^{-1} \pmb{P} $

LMS算法直接利用瞬态均方误差对瞬时抽头向量求梯度：

因此得到了LMS算法的迭代式：
$ y(n)= \pmb{w}^T (n)\pmb{x}(n) $

$ e(n)=d(n)-\pmb{w}^T (n)\pmb{x}(n) $

$ \pmb{w}(n+1)=\pmb{w}(n)+2\mu e(n)\pmb{x}(n) $
其中μ为迭代步长

而在FPGA实现时，可以直接使用FIR结构，框图如下所示：

另外还可以在matlab进行一下简单仿真，看一下LMS算法的效果：

close all;clear all;clc
N=2000;
n=8;
UU=randn(n,N);
VV=randn(1,N)*0.1+1;
wo=randn(n,1);
DD=wo'*UU+VV;

%% LMS
mu=0.01;
w_lms=randn(n,1);
for ii=1:N
   Err_LMS(ii)=norm(wo-w_lms);
   di=DD(ii);
   ui=UU(:,ii);
   ei=di-w_lms'*ui;
   w_lms=w_lms+2*mu*ei*ui;
   
end
plot(Err_LMS);

Verilog实现

顶层模块：

module FPGA_LMS(
    input sys_clk,
    input sys_rst_n,
    input [15:0] x_in,
    input [15:0] d_in,
    output [15:0] y_out,
    output [15:0] e_out
    );

    wire [127:0] coeff;

    assign e_out = d_in - y_out;

    filter filter1(
        .sys_clk(sys_clk),
        .sys_rst_n(sys_rst_n),
        .x_in(x_in),
        .coeff(coeff),
        .y_out(y_out)

    );

    coeff coeff1(
        .sys_clk(sys_clk),
        .sys_rst_n(sys_rst_n),
        .x_in(x_in),
        .e_in(e_out),
        .w_out(coeff)
    );


endmodule

输出更新模块：

module filter(
    input sys_clk,
    input sys_rst_n,
    input [15:0] x_in,
    input [127:0] coeff,
    output [15:0] y_out
    );

    // y(k)=w(k)*x(k)
    reg [127:0] x_in_r;
    wire [34:0] m0,m1,m2,m3,m4,m5,m6,m7;

    assign m0 = x_in_r[15:0] * coeff[15:0];
    assign m1 = x_in_r[31:16] * coeff[31:16] + m0;
    assign m2 = x_in_r[47:32] * coeff[47:32] + m1;
    assign m3 = x_in_r[63:48] * coeff[63:48] + m2;
    assign m4 = x_in_r[79:64] * coeff[79:64] + m3;
    assign m5 = x_in_r[95:80] * coeff[95:80] + m4;
    assign m6 = x_in_r[111:96] * coeff[111:96] + m5;
    assign m7 = x_in_r[127:112] * coeff[127:112] + m6;
    assign y_out = m7[15:0];

    always @(posedge sys_clk or negedge sys_rst_n) begin
        if(!sys_rst_n)
        begin
            x_in_r <= 128'd0;
        end
        else
        begin
            x_in_r <= {x_in_r[111:0],x_in};
        end
    end


endmodule

系数更新模块：

module coeff(
    input sys_clk,
    input sys_rst_n,
    input [15:0] x_in,
    input [15:0] e_in,
    output reg [127:0] w_out
    );
    // w(k+1)=w(k)+μ*e(k)*x(k)

    reg [127:0] x_in_r;
    wire [47:0] m0,m1,m2,m3,m4,m5,m6,m7;

    always @(posedge sys_clk or negedge sys_rst_n) begin
        if(!sys_rst_n)
        begin
            x_in_r <= 128'd0;
        end
        else
        begin
            x_in_r <= {x_in_r[111:0],x_in};
        end
    end

    always @(posedge sys_clk or negedge sys_rst_n) begin
        if(!sys_rst_n)
        begin
            w_out <= 128'd0;
        end
        else
        begin
            w_out[15:0] <= m0[15:0];
            w_out[31:16] <= m1[15:0];
            w_out[47:32] <= m2[15:0];
            w_out[63:48] <= m3[15:0];
            w_out[79:64] <= m4[15:0];
            w_out[95:80] <= m5[15:0];
            w_out[111:96] <= m6[15:0];
            w_out[127:112] <= m7[15:0];

        end

    end

assign m0 = x_in_r[15:0] * e_in + {{32{w_out[15]}},w_out[15:0]};
assign m1 = x_in_r[31:16] * e_in + {{32{w_out[31]}},w_out[31:16]};
assign m2 = x_in_r[47:32] * e_in + {{32{w_out[32]}},w_out[47:32]};
assign m3 = x_in_r[63:48] * e_in + {{32{w_out[63]}},w_out[63:48]};
assign m4 = x_in_r[79:64] * e_in + {{32{w_out[79]}},w_out[79:64]};
assign m5 = x_in_r[95:80] * e_in + {{32{w_out[95]}},w_out[95:80]};
assign m6 = x_in_r[111:96] * e_in + {{32{w_out[111]}},w_out[111:96]};
assign m7 = x_in_r[127:112] * e_in + {{32{w_out[127]}},w_out[127:112]};

endmodule

在Verilog代码里的乘法可以由别的方法实现，比如LUT、IP核、原语等
接下来就要验证代码是否能够实现功能，视频里采用的是正弦波去噪的案例，不过由于UP主没有提供测试的波形文件，因此这里我们自己制作一下，使用matlab即可：