自适应陷波器的原理与实现

前言：陷波器可以将某个频率成分滤除而保留其余信号，当然也可以通过设计带阻滤波器完成相同的功能，但带阻滤波器会对有用信号造成损失，且设计一个阻带很窄的滤波器消耗的资源也会很多。自适应陷波器则可以根据干扰频率信号的幅度和相位，自适应地调整滤波器的系数，跟踪其参数变化，有效滤除干扰信号，而不会损失有用信号。

目录

一、自适应陷波器的原理

二、MATLAB仿真

三、自适应滤波器的FPGA实现

四、实际测试

一、自适应陷波器的原理

采用LMS算法，在保证不损失其他有用信号的情况下，有效地抑制输入信号中某一频率的干扰，如需滤除多个频率，只需要增加相应的输入信号。

图1 可同时滤除两个频率干扰信号的自适应陷波器原理图

x(t)是叠加有两个频率干扰信号的输入信号，s(t)是需要保留的有用信号，在输入信号中，我们只知道频率w1和w2，但不知道信号的幅度A和相位θ，自适应陷波器可以估计两个未知的幅度和相位值，这需要我们提供两路相互正交的单频信号，如cos(w1t)、sin(w1t)，通过调整其权值，从而合成出与干扰信号完全相同的信号A1cos(w1t+θ1)，从而实现干扰信号的滤除。

通过上面原理图，并结合LMS算法可知，参考信号为r(t)，期望信号为x(t)，输出信号y(t)即为干扰信号，而误差则是有用信号s(t)，也就是自适应陷波器最终的输出结果。

二、MATLAB仿真

基于符号LMS算法对自适应陷波器进行仿真，干扰信号频率分别为50Hz、10Hz，有用信号为400Hz，抽样频率4KHz。

Matlab代码如下；

len = 4000;    
fs = 4000;     
t = 1:len;t = t/fs;
u = 1/128;   

f1 = 50;
f2 = 10;
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;x3;x4]; 

j1 = cos(2*pi*f1.*t + pi/3);
j2 = sin(2*pi*f2.*t + pi/6);

f = 400;
s = 2*cos(2*pi*f.*t);
d = j1 + j2 + s;

w = zeros(4,len+1);    
w(:,1) = ones(4,1)/2; 
e = zeros(1,len);      
dw = zeros(4,len);

for k=1:len
    y(k) = w(:,k)'*x(:,k);
    e(k) = d(k) - y(k);
    dw(:,k) =  u*sign(x(:,k))*e(k);
    w(:,k+1) = w(:,k) + dw(:,k);
end

disp_len = 500;
ax = 1:disp_len+1;
subplot(511);
plot(ax,s(len-disp_len:len));legend('有用信号');
subplot(512);
plot(ax,j1(len-disp_len:len));legend('50Hz信号');
subplot(513);
plot(ax,j2(len-disp_len:len));legend('10Hz信号');
subplot(514);
plot(ax,d(len-disp_len:len));legend('混合信号');
subplot(515);
plot(ax,e(len-disp_len:len));legend('滤波后的信号');

图2 仿真结果

可以看出滤除了混合信号中的两个干扰信号，输出信号与有用信号一致。

三、自适应滤波器的FPGA实现

由于数据经过了16bit量化，如上节介绍的那样，必须注意定点数小数点的位置，防止计算结果出错。

Verilog代码如下：

module NotchFilter
(
      input                  clk       ,
      input                  rst       ,
      input signed[15:0]     din       ,

      output signed[15:0]    dout
);

wire signed[15:0] Xin_Reg[3:0];
nco nco_inst1
(
      .phi_inc_i    (32'd8947849),
      .clk          (clk),
      .reset_n      (rst),
      .clken        (1'b1),
      .fsin_o       (Xin_Reg[1]),
      .fcos_o       (Xin_Reg[0]),
      .out_valid    ()
);

nco nco_inst2
(
      .phi_inc_i    (32'd1789570),      //10Hz
      .clk          (clk),
      .reset_n      (rst),
      .clken        (1'b1),
      .fsin_o       (Xin_Reg[3]),
      .fcos_o       (Xin_Reg[2]),
      .out_valid    ()
);

reg [2:0] count;
reg signed[15:0] Rin;
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
           count <= 3'd0;
      else if(count == 3'd5)
           begin
                 count <= 3'd0;
                 Rin <= din;
           end
      else
           count <= count + 1'b1;

reg signed[15:0] W_Reg[3:0];
reg signed[15:0] DW_Reg[3:0];
reg [2:0] k;
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
           begin
                 for(k=0;k<4;k=k+1)
                      W_Reg[k] <= 16'd1;        
           end
      else if(count == 3'd5)
           begin
                 for(k=0;k<4;k=k+1)
                      W_Reg[k] <= W_Reg[k] + DW_Reg[k];
           end

wire signed[31:0] Y_Reg[3:0];
genvar v;
generate
      for(v=0;v<=3;v=v+1)
           begin:u
                 mult u
                 (
                      .aclr ( !rst ),
                      .clock ( clk ),
                      .dataa ( Xin_Reg[v] ),
                      .datab ( W_Reg[v] ),
                      .result ( Y_Reg[v] )
                 );
           end
endgenerate

reg signed[34:0] Y_out;
reg signed[34:0] E_out;
always@(posedge clk or negedge rst)
      if(rst == 1'b0)
           begin
                 Y_out <= 35'd0;
                 E_out <= 35'd0;
           end
      else
           begin
                 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]};
                 E_out <= {{4{Rin[15]}},Rin,15'd0} - Y_out;
           end
assign dout = E_out[30:15];

always@(posedge clk or negedge rst)
      if(rst == 1'b0)
           begin
                 for(k=0;k<4;k=k+1)
                      DW_Reg[k] <= 16'd0;
           end
      else    
           begin
                 for(k=0;k<4;k=k+1)
                      if(E_out[34] == 1'b1)
                            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
endmodule

两路干扰信号由NCO核产生，数据输入速率为4KHz，系统时钟为24KHz，合理分配时钟周期以提高系统整体的运算速率，编写仿真文件，输入数据由外部文件读入。

图2 modelsim仿真结果

四、实际测试

编写板载测试文件，由于DA模块是8bit无符号型，这里注意有符号数据转换为无符号型，并截取高位。同样采用NCO核产生200Hz、50Hz和10Hz的叠加信号作为测试输入，观察经过陷波器滤波后的输出结果，预测只会保留200Hz信号。

测试信号波形：

图3 测试信号

滤波输出：

图4 输出信号

示波器输出如图所示，看出经自适应陷波器处理后，仅剩下200Hz的单频信号，而50Hz与10Hz信号已经自动被滤除。

至此，完成了自适应陷波器的设计与验证，用MATLAB对LMS算法进行了仿真验证，在使用FPGA实现该算法时要特别注意中间变量的数据范围，由此确定数据字长和小数点的位置，防止结果出错；同时还要根据各运算步骤的运算量合理分配时钟，以提高系统整体运算速率。

参考文献：

[1] 杜勇. 数字滤波器的MATLAB与FPGA实现——Altera/Verilog版（第2版）. 北京：电子工业出版社，2019.