陷波器顾名思义就是对特定频率的信号有着很强的衰减的滤波器，也即阻带带宽极窄的带阻滤波器。在传统的数字陷波器设计中，为了能使某一频率信号得到足够大的衰减，通常的做法就是把阶数选的足够高来达到很大的衰减；但同时计算量也变得更大了。而且设计的过程复杂，不利于动态的调整。为了解决上述存在的问题自适应陷波器孕育而生。博主在一次心电信号(ECG)工频干扰滤除的项目中就使用了自适应陷波器，其效果要远远好于传统的陷波器，它能在信干比-100dB的时候依然能够很好的将有用的信号提取出来。
当我们知道原始信号里的干扰信号频率是多少时（例如最常见的50Hz工频干扰），这时我们只需要知道这个干扰信号的相位和幅度，然后就可以完全的“再现”这个干扰信号，然后我们就可以直接的从原始信号中将其减去，从而就得到了我们想要的信号成分。这一过程实际上就是自适应陷波器的基本工作原理。关于自适应实现算法的具体原理可以参考博主之前写过的博客文章《浅谈自适应滤波器》，下面给出的是今天要介绍的自适应陷波器的结构图：
其中x(k)是带有特定频率干扰的信号，也即输入的原始信号，可见它是从自适应滤波器的期望信号端输入的；而sin(2π f0/fs k)和cos(2π f0/fs k)，是我们已知的频率为f0的干扰信号（其中fs是采样率），将它们分别乘以W1和W2进行适当的线性组合，就可以使其输出y(k)接近实际的干扰，最后输出的误差e(k)就是我们感兴趣的信号。如何进行W1和W2的选择这其实是自适应算法需要完成的工作，具体来说就是采用某种准则来构建一个关于误差e(k)的函数，通过是误差e(k)最小化来求得W1和W2，最常用的准则就是使误差的均方和最小。这也是本次博主要和大家分享的具体的自适应陷波器。
Matlab仿真
首先给出的是基于LMS的瞬时梯度估计算法的仿真代码：
clc;
clear all;
close all;
%%
%************************生成仿真信号**************************************
Fs = 500; %设置采样频率
t = 0:1/Fs:3;
t = t';
Size_t = size(t,1);
F1 = 7;
F2 = 13;
F3 = 23;
F4 = 50;
SNR = -100; %信干比 Unit:dBSignal = 10^(SNR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号
noise = 0.95*sin(2*pi*F4*t+pi/2);
Signal_noise = Signal + noise; %加入50Hz工频干扰%%
%*************************************************************************
M = 2; %定义FIR滤波器阶数
Signal_Len = Size_t; %定义信号数据的个数
niu = 1; %算法调节步长控制因子
y_out = zeros(Signal_Len,1); %滤波器输出
error_out = zeros(Signal_Len,1); %误差输出
w_out = zeros(Signal_Len,M); %系数输出
for i=1:Signal_Len%数据输入if i == 1 %如果是第一次进入w = zeros(M,1); %初始化滤波器抽头系数x = zeros(M,1); %初始化信号向量endd = Signal_noise(i); %输入新的期望信号x = [sin(2*pi*F4*(i-1)/Fs)cos(2*pi*F4*(i-1)/Fs)]; %输入新的信号矢量%算法正体y = x' * w; %计算滤波器输出error = d - y; %计算误差w_forward = w + niu * error * x; %计算滤波器系数向量%变量更替w = w_forward;%滤波结果存储y_out(i) = y;error_out(i) = error;w_out(i,:) = w';
end
%%
subplot(2,1,1);
plot(t,Signal);
title('原始信号');
xlabel('时间t/s');
subplot(2,1,2);
plot(t,Signal_noise);
title('加入干扰噪声的信号');
xlabel('时间t/s');figure;
subplot(2,1,1);
plot(t,y_out);
title('滤波器输出');
xlabel('时间t/s');
subplot(2,1,2);
plot(t,error_out,'b');
title('输出误差');
xlabel('时间t/s');
axis([0,3,-3*10^-5,3*10^-5]);figure;
plot(t(1:Signal_Len),w_out(1),'r',t(1:Signal_Len),w_out(2),'b');
title('自适应滤波器系数');
xlabel('时间t/s');
可以看到在信干比-100dB的情况下，有用的信号得到了不错的恢复。
接下来给出的是基于RLS的仿真代码：
clc;
clear all;
close all;
%%
%************************生成仿真信号**************************************
Fs = 500; %设置采样频率
t = 0:1/Fs:3;
t = t';
Size_t = size(t,1);
F1 = 7;
F2 = 13;
F3 = 23;
F4 = 50;
SIR = -100; %信干比 Unit:dBSignal = 10^(SIR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号
noise = 0.95*sin(2*pi*F4*t + pi/4);
Signal_noise = Signal + noise; %加入50Hz工频干扰%%
%*************************************************************************
M = 2; %定义FIR滤波器阶数
lamda = 0.90; %定义遗忘因子
Signal_Len = Size_t; %定义信号数据的个数
I = eye(M); %生成对应的单位矩阵
c = 0.01; %小正数 保证矩阵P非奇异
y_out = zeros(Signal_Len,1);
Eta_out = zeros(Signal_Len,1);
w_out = zeros(Signal_Len,M);
for i=1:Signal_Len%输入数据if i == 1 %如果是第一次进入P_last = I/c;w_last = zeros(M,1); endd = Signal_noise(i); %输入新的期望信号x = [sin(2*pi*F4*(i-1)/Fs)cos(2*pi*F4*(i-1)/Fs)]; %输入新的信号矢量%算法正体K = (P_last * x)/(lamda + x'* P_last * x); %计算增益矢量y = x'* w_last; %计算FIR滤波器输出Eta = d - y; %计算估计的误差w = w_last + K * Eta; %计算滤波器系数矢量P = (I - K * x')* P_last/lamda; %计算误差相关矩阵%变量更替P_last = P;w_last = w;%滤波结果存储y_out(i) = y;Eta_out(i) = Eta;w_out(i,:) = w';
end
%%
figure;
subplot(3,1,1);
plot(t,Signal);
title('原始信号');
xlabel('时间t/s');
subplot(3,1,2);
plot(t,Signal_noise);
string = ['加入50Hz工频干扰的信号,信干比SIR = ',num2str(SIR),'dB'];
title(string);
xlabel('时间t/s');
subplot(3,1,3);
plot(t,Eta_out,'b',t,Signal,'r');
title('输出误差');
xlabel('时间t/s');
axis([0,3,-2*10^-5,2*10^-5]);%
可以看到仿真的结果相比于LMS的要好一些。
实际项目运用情况简介
心电信号(ECG 10mV以下)以及脑电信号(EEG 100uV以下)都是非常微弱的信号，其频带在0.1Hz到80Hz之间，在采集的过程中周围的各种50Hz工频干扰及其倍频三倍频等干扰会通过大地，空气等介质传播到人体上，进而被采集进来形成干扰。其强度一般是远大于心电信号的，所以要进行后一步详细的分析，干扰必须要被滤除掉。这时利用自适应陷波器是最为有效的方法。
当然这时用的自适应陷波器和之前介绍的又有所区别，主要区别是现在要去除的干扰频率有三个，所以需要对上述的滤波器进行改进，一种方法是采用串行的方式滤除，即首先滤除50Hz的干扰，然后滤除100Hz的，最后滤除150Hz的。另一种是采用并行的方式同时滤除这三种干扰，下面我给出的是这种方式的结构图：
下面给出的是采集的一段心电信号处理的结果：
可以看到原始信号中的50Hz干扰非常强，如果做频谱分析的话会发现还有100Hz和150Hz的干扰因此我设计的是滤除50Hz、100Hz和150Hz的自适应陷波器，这时出来的结果已经是把50Hz、100Hz和150Hz的干扰滤除干净了但是波形视乎还是不够平滑，这主要是大于100Hz的噪声引起的，因此自适应陷波后还要将其通过一个截止频率为80Hz的LPF，也就是低通滤波器，这时的输出就相当平滑了，但是还有一个值比较大的近乎直流的浮动的电压，这个电压主要是人体的静电引起的，因此还需要估计一下这个浮动的电压，然后将其减掉从而最终得到了干净的心电信号。
之前学习自适应一直是做的仿真，这次终于把它运用到了实际的项目中，效果非常显著，更有说服力了。如果大家有什么疑问的可以联系我，我们一起探讨。
