LMS自适应滤波算法（MATLAB仿真+fpga实现）
·
MATLAB代码仿真
1
g=100; %蒙特-卡罗仿真统计次数
N=1024; % 输入信号序列长度
k=128; %fir滤波器长度
pp=zeros(g,N-k); %每次循环仿真的误差信号存于矩阵pp中 zero用于生成矩阵
u=1/256; %步长因子
snr=[3,-3]; %输入信号信噪比参数 一个3db，一个-3db
t=1:N;
s=sin(0.1*pi*t); %生成正弦波信号
xn=zeros(1,N); %存放输入信号
y=zeros(1,N); %存放输出信号
w=zeros(1,k); %存放权值信号
e=zeros(1,N); %存放误差信号
for type=1:4
for q=1:g
noise=rand(1,length(s)); %随机噪声
if type==1
SNR=snr(1);d=s; % 场景1：期望信号为正弦波，SNR=3dB
elseif type==2
SNR=snr(1);d=sqrt(10^(-SNR/10))*noise; % 场景2：期望信号为噪声，SNR=3dB
elseif type==3
SNR=snr(2);d=s; % 场景3：期望信号为正弦波，SNR=-3dB
else
SNR=snr(2);d=sqrt(10^(-SNR/10))*noise; % 场景4：期望信号为噪声，SNR=-3dB
end
xn=sqrt(10^(-SNR/10))*noise+s;
y(1:k)=xn(1:k); % 初始化滤波器输出
%LMS
for i=(k+1):N
XN=xn((i-k+1):(i));
y(i)=w*XN'; %注意'是矩阵行列转置
e(i)=d(i)-y(i); %误差
w=w+u*e(i)'*XN; % 更新权值
end
pp(q,:)=(e(k+1:N)).^2; %求每次仿真后误差信号的平方值
end
figure(1);
subplot(311);
plot(s(300:450)); %截取一段信号绘图
title('信号s时域波形');
if type==1
subplot(312);plot(xn(300:450));
title('信号s加噪声后（snr=3db）');
elseif type==3
subplot(313);plot(xn(300:450));
title('信号s加噪声后（snr=-3db）');
end
%求取各次循环仿真的误差统计均值
for b=1:N-k
bi(b)=sum(pp(:,b))/g;
end
%绘制自适应滤波后的输出信号
figure(2);
if type==1
subplot(311);
plot(y(300:450));title('自适应滤波输出，期望为sin，snr=3db');
elseif type==3
subplot(312);
plot(y(300:450));title('自适应滤波输出，期望为sin，snr=-3db');
elseif type==4
subplot(313);y=xn-y; %期望信号为噪声，系统相当于干扰抵消系统
plot(y(300:450));title('自适应滤波输出，期望为噪声，snr=3db');
end
%绘制误差信号
figure(3);
if type==1
subplot(311);
plot(bi(1:100));title('误差均方信号，snr=3db，期望为sin');
elseif type==3
subplot(312);
plot(bi(1:100));title('误差均方信号，snr=-3db，期望为sin');
elseif type==4
subplot(313);
plot(bi(1:100));title('误差均方信号，snr=3db，期望为噪声');
end
end
仿真结果
从图中可以看出：期望信号选择噪声信号或者正弦波信号均可得到较好结果
2
将输入的信号改一下，改成s=sin(0.1*pi*t)+sin(0.3*pi*t);，参考信号改成s0=sin(0.1*pi*t);
g=100; %蒙特-卡罗仿真统计次数
N=1024; % 输入信号序列长度
k=128; %fir滤波器长度
pp=zeros(g,N-k); %每次循环仿真的误差信号存于矩阵pp中 zero用于生成矩阵
u=1/256; %步长因子
snr=[3,-3]; %输入信号信噪比参数 一个3db，一个-3db
t=1:N;
s0=sin(0.1*pi*t); %参考信号
s=sin(0.1*pi*t)+sin(0.3*pi*t); %生成正弦波信号
xn=zeros(1,N); %存放输入信号
y=zeros(1,N); %存放输出信号
w=zeros(1,k); %存放权值信号
e=zeros(1,N); %存放误差信号
for type=1:4
for q=1:g
noise=rand(1,length(s)); %随机噪声
if type==1
SNR=snr(1);d=s0; % 场景1：期望信号为正弦波，SNR=3dB
elseif type==2
SNR=snr(1);d=sqrt(10^(-SNR/10))*noise; % 场景2：期望信号为噪声，SNR=3dB
elseif type==3
SNR=snr(2);d=s0; % 场景3：期望信号为正弦波，SNR=-3dB
else
SNR=snr(2);d=sqrt(10^(-SNR/10))*noise; % 场景4：期望信号为噪声，SNR=-3dB
end
xn=sqrt(10^(-SNR/10))*noise+s;
y(1:k)=xn(1:k); % 初始化滤波器输出
%LMS
for i=(k+1):N
XN=xn((i-k+1):(i));
y(i)=w*XN'; %注意'是矩阵行列转置
e(i)=d(i)-y(i); %误差
w=w+u*e(i)'*XN; % 更新权值
end
pp(q,:)=(e(k+1:N)).^2; %求每次仿真后误差信号的平方值
end
figure(1);
subplot(311);
plot(s(300:450)); %截取一段信号绘图
title('信号s时域波形');
if type==1
subplot(312);plot(xn(300:450));
title('信号s加噪声后（snr=3db）');
elseif type==3
subplot(313);plot(xn(300:450));
title('信号s加噪声后（snr=-3db）');
end
%求取各次循环仿真的误差统计均值
for b=1:N-k
bi(b)=sum(pp(:,b))/g;
end
%绘制自适应滤波后的输出信号
figure(2);
if type==1
subplot(311);
plot(y(300:450));title('自适应滤波输出，期望为sin，snr=3db');
elseif type==3
subplot(312);
plot(y(300:450));title('自适应滤波输出，期望为sin，snr=-3db');
elseif type==4
subplot(313);y=xn-y; %期望信号为噪声，系统相当于干扰抵消系统
plot(y(300:450));title('自适应滤波输出，期望为噪声，snr=3db');
end
%绘制误差信号
figure(3);
if type==1
subplot(311);
plot(bi(1:100));title('误差均方信号，snr=3db，期望为sin');
elseif type==3
subplot(312);
plot(bi(1:100));title('误差均方信号，snr=-3db，期望为sin');
elseif type==4
subplot(313);
plot(bi(1:100));title('误差均方信号，snr=3db，期望为噪声');
end
end
仿真结果
3自适应陷滤波器
clear all;
len=4000; %数据长度
fs=12.5*10^6; %采样频率
u=1/128; %步长因子
t=1:len;
t=t/fs;
f0=500*10^3; %500khz
%
f1=200*10^3; %200khz
f2=50*10^3; %50khz
%生成4路参考信号 两两正交
x1=cos(2*pi*f1.*t);
x2=sin(2*pi*f1.*t);
%x3=cos(2*pi*f2.*t);
%x4=sin(2*pi*f2.*t);
x=[x1;x2];
%生成干扰信号
j1=2*cos(2*pi*f1.*t+pi/3);
j2=2*sin(2*pi*f2.*t+pi/6);
%生成有用信号
s=cos(2*pi*f0.*t); %500khz
%将干扰信号混入有用信号
d=j1+s;
%LMS算法中间变量初始化
w=zeros(2,len+1);
w(:,1)=ones(2,1)/2;
e=zeros(1,len);
aw=zeros(2,len);
%LMS
for k=1:len
y(k)=w(:,k)'*x(:,k);
e(k)=d(k)-y(k);
%aw(:,k)=2*u*x(:,k)*conj(e(k)); %LMS算法
aw(:,k)=2*u*sign(x(:,k))*conj(e(k)); %符号LMS算法
w(:,k+1)=w(:,k)+aw(:,k);
end
%绘图
disp_len=1000; %显示1000个数据点
ax=1:disp_len+1;
subplot(611);
plot(ax,s(len-disp_len:len));legend('有用信号');
subplot(612);
plot(ax,j1(len-disp_len:len));legend('200khz干扰信号');
subplot(613);
plot(ax,j2(len-disp_len:len));legend('50khz干扰信号');
subplot(614);
plot(ax,d(len-disp_len:len));legend('叠加干扰信号');
subplot(615);
plot(ax,e(len-disp_len:len));legend('滤除干扰后的信号eout');
subplot(616);
plot(ax,y(len-disp_len:len));legend('滤除干扰后的信号yout');
仿真结果
选择eout可能会有问题，特定情境下
lfsr模块，用于仿真时的噪声模拟
module lfsr(
input clk,
input en,
input rst_n,
output [15:0] lfsr_out
);
//Create a 16-bit linear feedback shift register with
//maximal polynomial x^16 + x^14 + x^13 + x^11 + 1.
reg [15:0]lfsr = 16'd1;
wire feedback;
assign feedback = ((lfsr[15] ^ lfsr[13]) ^ lfsr[12]) ^ lfsr[10];
assign lfsr_out = lfsr;
//Update the linear feedback shift register.
always @(posedge clk) begin
if(!rst_n)
lfsr <= 16'd1;
else if(en)
lfsr <= {lfsr[14:0], feedback};
end
endmodule
verilog代码
此处是一个示例，输入一个方波或是三角波，输出基波的正弦波
notch_filter自适应滤波器
//////////////////////////////////////////////////////////////////////////////////
// Company: HDU
// Engineer: Aqua
// Create Date: 2025/04/24 21:38:31
// Design Name:
// Module Name: notch_filter 自适应陷滤波器
//////////////////////////////////////////////////////////////////////////////////
module notch_filter(
input clk, // fpga系统时钟信号，速率为数据速率的6倍 75Mhz
input clk_50, // ila
input rst_n,
input [23:0] amp, // 分频系数
input signed [15:0] din,
input [33:0] freq, // 想要滤出的信号频率
output reg signed [34:0] e_out, // 误差信号e_out
output reg signed [34:0] y_out, // 滤波器输出yout
output reg signed [15:0] dout // 输出滤波器输出
);
// ---------------------基准信号生成----------------------
// 例化dds ip ，产生2路基准信号
wire signed [15:0] xin_reg[1:0];
wire [31:0] sin200k,sin50k;
// f = data*fclk/2^B
// f为dds出波频率，fclk为dds驱动时钟，B为相位累加字位宽
// 例如此处，fclk = 75Mhz，B = 16，f = 200k 计算出data为174.7取175
// 期望信号
dds_0 dds_0_inst(
.aclk(clk_50),
.s_axis_config_tvalid(1'b1),
.s_axis_config_tdata(('d65536*freq)/50_000_000),
.m_axis_data_tdata({xin_reg[1],xin_reg[0]})
);
// // 50kHz的正弦波
// dds_0 dds_0_inst_1(
// .aclk(clk),
// .s_axis_config_tvalid(1'b1),
// .s_axis_config_tdata(16'd44),
// .m_axis_data_tdata({xin_reg[3],xin_reg[2]})
// );
// 3位计数器，计数周期为6
reg [2:0] cnt;
reg signed [15:0] rin;
always@(posedge clk or negedge rst_n) begin
if(!rst_n)
cnt <= 3'd0;
else begin
if(cnt == 3'd5) begin
cnt <= 3'd0;
rin <= din; // 每6个时钟周期取一次采集的波形数据
end
else
cnt <= cnt + 1'b1;
end
end
// 权值数据
reg signed [15:0] w_reg[1:0]; // 这是一个数组，包含4个16位的有符号数
reg signed [15:0] dw_reg[1:0];
reg [2:0] k;
reg [2:0] k1;
always@(posedge clk or negedge rst_n) begin
if(!rst_n) begin
//初始化位移寄存器的值为0
k1 <= 3'd0;
for(k1=0;k1<2;k1=k1+1) begin
// 初始化权值为1
w_reg[k1] <= 16'b0011111111111111111;
end
end
else begin
if(cnt == 3'd5) begin
k1 <= 3'd0;
for(k1=0;k1<2;k1=k1+1) begin
w_reg[k1] <= w_reg[k1] + dw_reg[k1];
end
end
end
end
// 4个乘法器，2级流水线，并行完成 权值与基准信号的乘法运算
wire signed [31:0] y_reg[1:0];
mult mult_0(
.SCLR(!rst_n),
.CLK(clk),
.A(xin_reg[0]),
.B(w_reg[0]),
.P(y_reg[0])
);
mult mult_1(
.SCLR(!rst_n),
.CLK(clk),
.A(xin_reg[1]),
.B(w_reg[1]),
.P(y_reg[1])
);
// mult mult_2(
// .SCLR(!rst_n),
// .CLK(clk),
// .A(xin_reg[2]),
// .B(w_reg[2]),
// .P(y_reg[2])
// );
// mult mult_3(
// .SCLR(!rst_n),
// .CLK(clk),
// .A(xin_reg[3]),
// .B(w_reg[3]),
// .P(y_reg[3])
// );
always@(posedge clk or negedge rst_n) begin
if(!rst_n) begin
// 初始化移位寄存器的值为0
y_out <= 35'd0;
e_out <= 35'd0;
end
else begin
// yout 在一个时钟周期内完成更新
y_out <= {{3{y_reg[0][31]}},y_reg[0]} + {{3{y_reg[1][31]}},y_reg[1]} ; //+ {{3{y_reg[2][31]}},y_reg[2]} + {{3{y_reg[3][31]}},y_reg[3]}
// eout 在二个时钟周期内完成更新
e_out <= {{4{rin[15]}},rin,15'd0} - y_out; // 误差
end
end
always@(posedge clk) begin
if(cnt == 3'd4)
dout <= y_out[30:15]; // 取高位输出
end
// 根据误差信号e_out的符号，求dw的值，延时一个时钟周期
always@(posedge clk or negedge rst_n) begin
if(!rst_n) begin
// 初始化移位寄存器的值为0
k <= 3'd0;
for(k=0;k<2;k=k+1) begin
dw_reg[k] <= 16'd0;
end
end
else begin
k <= 3'd0;
for(k=0;k<2;k=k+1) begin
if(e_out[34])
dw_reg[k] <= -{{7{xin_reg[k][15]}},xin_reg[k][15:7]};
else
dw_reg[k] <= {{7{xin_reg[k][15]}},xin_reg[k][15:7]};
end
end
end
ila_0 ila_0_inst(
.clk(clk_50),
.probe0(din), // 16bit
.probe1(e_out), // 34bit
.probe2(y_out), // 34bit
.probe3(dout), // 16bit
.probe4(xin_reg[1]), // 16bit
.probe5(freq) // 34bit
);
endmodule
div
时钟分频模块，奇偶分频占空比50%，此处主要是用于产生PLL产生不了的一些比较低的时钟频率，并且比较方便去修改notch_filter的工作频率，因为这个的工作的频率和adc采样频率都得根据需要处理信号的频率来一点一点尝试
// 对于奇数倍和偶数倍分频，且保证占空比为50%
module div
(
input clk,
input rst_n,
input [23:0]amp, //分频系数
output clk_div
);
// 奇数分频
reg [15:0] cnt;
reg clk_pos,clk_neg;
wire clk_div_j;
assign clk_div_j = clk_pos ^ clk_neg;
// 偶数分频
reg clk_div_o;
// 输出
assign clk_div = (amp%2 == 1)?clk_div_j:clk_div_o; //奇数偶数分频选择输出
// 奇数分频
always @(posedge clk or negedge rst_n)begin
if(!rst_n)
cnt <= 0;
else if(cnt==amp-1)
cnt <= 0;
else
cnt <= cnt + 1;
end
always @(posedge clk or negedge rst_n)begin
if(!rst_n)
clk_pos <= 0;
else if(cnt==amp-1)
clk_pos <= ~clk_pos;
end
always @(negedge clk or negedge rst_n)begin
if(!rst_n)
clk_neg <= 0;
else if(cnt==(amp>>1))
clk_neg <= ~clk_neg;
end
// 偶数分频
reg [23:0] cnt_o;
always@(posedge clk or negedge rst_n) begin
if(!rst_n)
cnt_o <= 24'd0;
// else if(amp == 4 && cnt_o == 1)
else if((amp%2==0) && (cnt_o == amp/2-1))
cnt_o <= 24'd0;
else
cnt_o <= cnt_o + 1'b1;
end
always@(posedge clk or negedge rst_n) begin
if(!rst_n)
clk_div_o <= 1'b0;
// else if(amp == 4 && cnt_o == 1)
// clk_div_o <= ~clk_div_o;
else if(amp ==2) // 2分频比较特殊
clk_div_o <= ~clk_div_o;
else if((amp%2==0) && (cnt_o == amp/2-1)) // 其余偶数分频
clk_div_o <= ~clk_div_o;
else
clk_div_o <= clk_div_o;
end
endmodule
main顶层模块
module main(
input clk,
input rst_n,
input [11:0] adc_data, // adc采集数据
output adc_clk, //时钟为LMS的1/6
output [13:0] dac_data, // dac输出数据
output dac_clk
);
wire clk_75m,clk_12m5,clk_100m,clk_5m,clk_30m;
wire clk_lms;
// assign adc_clk = ~clk;
assign dac_clk = clk;
// 仿真时添加随机噪声
//wire [15:0] lfsr_out; // 随机噪声
wire [33:0] freq;
wire signed [15:0] din;
wire signed [15:0] dout; // 滤波器输出
wire signed [15:0] adc_data1;
wire signed [15:0] adc_data2;
//assign din = adc_data1/2 + adc_data2/2;
// adc数据转化成有符号数
assign din = {adc_data,4'b0}-16'd32768;
// 转化成无符号数给DAC
assign dac_data = dout[15:2] + 14'd8192;
wire signed [15:0] useless_data1;
wire signed [15:0] useless_data2;
// f = data*fclk/2^B
// f为dds出波频率，fclk为dds驱动时钟，B为相位累加字位宽
// 例如此处，fclk = 75Mhz，B = 16，f = 200k 计算出data为174.7取175
// // 200kHz的正弦波
// dds_0 dds_0_inst(
// .aclk(clk_12m5),
// .s_axis_config_tvalid(1'b1),
// .s_axis_config_tdata(16'd1049),
// .m_axis_data_tdata({adc_data1,useless_data1})
// );
// dds_0 dds_0_inst1(
// .aclk(clk_12m5),
// .s_axis_config_tvalid(1'b1),
// .s_axis_config_tdata(16'd264),
// .m_axis_data_tdata({adc_data2,useless_data2})
// );
wire locked;
clock0 clock0_inst(
.clk_75m(clk_75m),
.clk_12m5(clk_12m5),
.clk_100m(clk_100m),
.clk_out4(clk_5m),
.clk_out5(clk_30m),
.resetn(rst_n),
.clk_in1(clk),
.locked(locked)
);
reg [23:0] amp = 24'd5;
div div_inst1(
.clk(clk_30m),
.rst_n(rst_n),
.amp(amp),
.clk_div(clk_lms)
);
div div_inst2(
.clk(clk_5m),
.rst_n(rst_n),
.amp(amp),
.clk_div(adc_clk)
);
notch_filter notch_filter_inst(
.clk(clk_lms),
.clk_50(clk),
.rst_n(rst_n),
.amp(amp),
.freq(freq),
.din(din),
.dout(dout)
);
// 测量基波的频率，以此生成该频率正弦波的基准信号
freq_meter_calc freq_meter_calc_inst(
.sys_clk(clk),
.clk_stand(clk_100m),
.sys_rst_n(rst_n),
.clk_test(din[15]),
.freq(freq)
);
endmodule
效果：
显然这不是一个合格的正弦波，因此还需要过一个低通滤波器，此处的高频分量正好是notch_filter的工作时钟频率的6分频（这个模块处理过程需要6个时钟周期），
更多推荐
4
0- 0
已为社区贡献3条内容
所有评论(0)