自适应陷波器的设计

1. 自适应陷波器原理

陷波器也是一种滤波器，其主要特点在于对输入信号的某些已知频率信号进行滤除。这里所说的已知频率信号，是指已知干扰信号的频率，不知信号的相位和幅度信息。

例如，由于我国采用的是50Hz的交流电，所以在需要对信号进行采集处理和分析时，经常会存在50Hz工频信号，对有用信号的处理造成很大的干扰，因此很有必要设计50Hz的陷波器，对这个已知频率信号进行滤除。所谓自适应陷波器，就是可以根据干扰信号频率的相位和幅度，自适应地调整滤波器系数，跟踪其参数变化，保持对干扰信号的有效滤除。

陷波器的设计方法较多，除了本文所讲之外，还可以直接设计一个带阻滤波器，使干扰信号的频率无法通过。显然，由于硬件电路的滤波器无法做到很窄的过渡带，因此会对有用信号造成损失。

本文所讲的自适应陷波器，其目的是采用LMS算法，在保证其他有用信号不损失的情况下，有效抑制输入信号中某一频率的干扰（如果需要滤除多个频率，只需增加相应的运算资源）。下图为同时滤除两个干扰频率（${\omega }_1、{\omega }_2$）信号的自适应陷波器原理图。

图中，$x(t)$是叠加有两个干扰频率的输入信号，自适应陷波器就是要滤出${\omega }_1、{\omega }_2$这两个频率的干扰信号，$s(t)$是需要保留的有用信号，可以是任意形式的信号。需要说明的是，在输入信号中，我们仅知道干扰信号的频率是${\omega }_1、{\omega }_2$，不知道其具体幅度$A_1、A_2$及相位值${\theta }_1、{\theta }_2$，自适应算法（LMS）的目的就是要估计两个未知的幅度值及相位值。为了估计出单个频率的幅度值和相位值，我们需要用两路相互正交的单频信号，如$cos({\omega }_{1}t)、sin({\omega }_{2}t)$，通过调整其权值$w_1、w_2$（均为实数），就可以合成与干扰信号$A_{1}cos({\omega }_{1}t+{\theta }_1)$完全相同的信号；由于需要有效滤除两个单频干扰信号，因此需要输入4路参考信号，通过调整相应权值实现单频干扰信号的有效滤除。

经LMS算法调整后的权值，与4路参考信号乘加输出的值，即为两个单频干扰的估计值，即$y(t)=\widehat{A_1}cos({\omega }_{1}t+\widehat{{\theta }_1})+\widehat{A_2}cos({\omega }_{2}t+\widehat{{\theta }_2})$。LMS算法的误差信号$e(t)$就是滤除了干扰的有用信号。

经过上面的讨论，显然，如果干扰信号只是一个单频信号，如50Hz的工频干扰，则参考信号只需要$cos100\pi t、sin100\pi t$两路信号；如果需要滤除3个单频干扰，则参考信号需要6路。

数字LMS算法中，权值均是对输入信号进行调制，期望信号$d(n)$均是确知的信号，不需要权值调整。在进行自适应陷波器的LMS算法仿真及实现时，需要将输入信号$x(t)$当做LMS算法的期望信号$d(n)$；将4路参考信号$r(t)$当做LMS算法的输入信号$x(n)$；将权值调整后的乘加值$y(t)$当做LMS算法的估计值$y(n)$；将误差信号$e(t)$当做LMS算法的误差信号$e(n)$，也是自适应陷波器的最终输出结果。

2. 自适应陷波器的MATLAB仿真

MATLAB仿真程序中，采样频率为4KHz，输入信号中有两个单频干扰信号$(f_1=50Hz、f_2=10Hz)$，要求对这两个单频干扰信号进行滤除。仿真有用信号为两种形式信号（400Hz的单频信号、随机信号）情况下的陷波器算法。LMS算法的步长因子设置为1/128。
源码

len=4000; %数据长度
fs=4000;  %采样频率

t=1:len;t=t/fs;

%两个单频干扰信号的频率/参考信号的频率
f1=50;f2=10;
%生成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;x3;x4];

%生成干扰信号
J1=2*cos(2*pi*f1.*t+pi/3);
J2=2*sin(2*pi*f2.*t+pi/6);

%生成有用信号
s=cos(2*pi*400.*t); %频率为400Hz的单频信号
% s=rand(1,len);     %或者为：随机信号

%生成混有干扰信号的输入信号
d=J1+J2+s;

%LMS算法中间变量初始化
w=zeros(4,len+1);
w(:,1)=ones(4,1)/2;
e=zeros(1,len);
aw=zeros(4,len);
%LMS算法
for k=1:len
    y(k)=w(:,k)'*x(:,k);
    e(k)=d(k)-y(k);
    %aw(:,k)=2*(1/128)*x(:,k)*conj(e(k));     %LMS算法
    aw(:,k)=2*(1/128)*sign(x(:,k))*conj(e(k));%符号LMS算法
    w(:,k+1)=w(:,k)+aw(:,k);
end

%绘图
disp_len=300;  %显示300个数据点
% disp_len=2000; %显示2000个数据点
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('滤除干扰后的信号');

仿真结果
有用信号为单频信号的陷波器仿真波形（显示300个数据）

有用信号为单频信号的陷波器仿真波形（显示2000个数据）

有用信号为随机信号的陷波器仿真波形（显示300个数据）

有用信号为随机信号的陷波器仿真波形（显示2000个数据）

从图中可以看出，陷波器仍然有效滤除了两个单频干扰，输出的信号与有用信号完全一致。