浅谈自适应滤波器---（自适应陷波器）

  原文地址：http://blog.csdn.net/hj199404182515/article/details/71527011

  陷波器顾名思义就是对特定频率的信号有着很强的衰减的滤波器，也即阻带带宽极窄的带阻滤波器。在传统的数字陷波器设计中，为了能使某一频率信号得到足够大的衰减，通常的做法就是把阶数选的足够高来达到很大的衰减；但同时计算量也变得更大了。而且设计的过程复杂，不利于动态的调整。为了解决上述存在的问题自适应陷波器孕育而生。博主在一次心电信号(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的瞬时梯度估计算法的仿真代码：

[plain] view plain copy

 print?

1. clc;  
2. clear all;  
3. close all;  
4. %%  
5. %************************生成仿真信号**************************************  
6. Fs = 500;                                                     %设置采样频率  
7. t = 0:1/Fs:3;    
8. t = t';  
9. Size_t = size(t,1);  
10. F1 = 7;  
11. F2 = 13;  
12. F3 = 23;  
13. F4 = 50;  
14. SNR = -100;            %信干比 Unit:dB  
15.   
16. Signal = 10^(SNR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号  
17. noise = 0.95*sin(2*pi*F4*t+pi/2);  
18. Signal_noise = Signal + noise;                               %加入50Hz工频干扰  
19.         
20. %%  
21. %*************************************************************************  
22. M = 2;                                    %定义FIR滤波器阶数  
23. Signal_Len = Size_t;                      %定义信号数据的个数  
24. niu = 1;                                  %算法调节步长控制因子  
25. y_out = zeros(Signal_Len,1);              %滤波器输出  
26. error_out = zeros(Signal_Len,1);          %误差输出              
27. w_out = zeros(Signal_Len,M);              %系数输出  
28. for i=1:Signal_Len  
29.     %数据输入  
30.     if i == 1           %如果是第一次进入  
31.         w = zeros(M,1); %初始化滤波器抽头系数  
32.         x = zeros(M,1); %初始化信号向量  
33.     end  
34.     
35.     d = Signal_noise(i);                     %输入新的期望信号  
36.     x = [sin(2*pi*F4*(i-1)/Fs)  
37.          cos(2*pi*F4*(i-1)/Fs)];                 %输入新的信号矢量  
38.       
39.     %算法正体  
40.     y = x' * w;                              %计算滤波器输出  
41.     error = d - y;                           %计算误差  
42.     w_forward = w + niu * error * x;         %计算滤波器系数向量  
43.     %变量更替  
44.     w = w_forward;  
45.     %滤波结果存储  
46.     y_out(i) = y;  
47.     error_out(i) = error;  
48.     w_out(i,:) = w';  
49. end  
50. %%  
51. subplot(2,1,1);  
52. plot(t,Signal);  
53. title('原始信号');  
54. xlabel('时间t/s');  
55. subplot(2,1,2);  
56. plot(t,Signal_noise);  
57. title('加入干扰噪声的信号');  
58. xlabel('时间t/s');  
59.   
60. figure;  
61. subplot(2,1,1);  
62. plot(t,y_out);  
63. title('滤波器输出');  
64. xlabel('时间t/s');  
65. subplot(2,1,2);  
66. plot(t,error_out,'b');  
67. title('输出误差');  
68. xlabel('时间t/s');  
69. axis([0,3,-3*10^-5,3*10^-5]);  
70.   
71. figure;  
72. plot(t(1:Signal_Len),w_out(1),'r',t(1:Signal_Len),w_out(2),'b');  
73. title('自适应滤波器系数');  
74. xlabel('时间t/s');  

运行一下得到如下的结果：

可以看到在信干比-100dB的情况下，有用的信号得到了不错的恢复。

接下来给出的是基于RLS的仿真代码：

[plain] view plain copy

 print?

1. clc;  
2. clear all;  
3. close all;  
4. %%  
5. %************************生成仿真信号**************************************  
6. Fs = 500;                                                     %设置采样频率  
7. t = 0:1/Fs:3;    
8. t = t';  
9. Size_t = size(t,1);  
10. F1 = 7;  
11. F2 = 13;  
12. F3 = 23;  
13. F4 = 50;  
14. SIR = -100;            %信干比 Unit:dB  
15.   
16. Signal = 10^(SIR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号  
17. noise = 0.95*sin(2*pi*F4*t + pi/4);  
18. Signal_noise = Signal + noise;                               %加入50Hz工频干扰  
19.        
20. %%  
21. %*************************************************************************  
22. M = 2;                         %定义FIR滤波器阶数  
23. lamda = 0.90;                  %定义遗忘因子       
24. Signal_Len = Size_t;           %定义信号数据的个数  
25. I = eye(M);                    %生成对应的单位矩阵  
26. c = 0.01;                      %小正数 保证矩阵P非奇异  
27. y_out = zeros(Signal_Len,1);  
28. Eta_out = zeros(Signal_Len,1);  
29. w_out = zeros(Signal_Len,M);  
30. for i=1:Signal_Len  
31.     %输入数据  
32.     if i == 1                 %如果是第一次进入  
33.         P_last = I/c;  
34.         w_last = zeros(M,1);    
35.     end  
36.     d = Signal_noise(i);                         %输入新的期望信号  
37.     x = [sin(2*pi*F4*(i-1)/Fs)  
38.          cos(2*pi*F4*(i-1)/Fs)];                 %输入新的信号矢量  
39.     %算法正体  
40.     K = (P_last * x)/(lamda + x'* P_last * x);   %计算增益矢量  
41.     y = x'* w_last;                          %计算FIR滤波器输出  
42.     Eta = d - y;                             %计算估计的误差  
43.     w = w_last + K * Eta;                    %计算滤波器系数矢量  
44.     P = (I - K * x')* P_last/lamda;          %计算误差相关矩阵  
45.     %变量更替  
46.     P_last = P;  
47.     w_last = w;  
48.     %滤波结果存储  
49.     y_out(i) = y;  
50.     Eta_out(i) = Eta;  
51.     w_out(i,:) = w';  
52. end  
53. %%  
54. figure;  
55. subplot(3,1,1);  
56. plot(t,Signal);  
57. title('原始信号');  
58. xlabel('时间t/s');  
59. subplot(3,1,2);  
60. plot(t,Signal_noise);  
61. string = ['加入50Hz工频干扰的信号,信干比SIR = ',num2str(SIR),'dB'];  
62. title(string);  
63. xlabel('时间t/s');  
64. subplot(3,1,3);  
65. plot(t,Eta_out,'b',t,Signal,'r');  
66. title('输出误差');  
67. xlabel('时间t/s');  
68. 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，也就是低通滤波器，这时的输出就相当平滑了，但是还有一个值比较大的近乎直流的浮动的电压，这个电压主要是人体的静电引起的，因此还需要估计一下这个浮动的电压，然后将其减掉从而最终得到了干净的心电信号。

        之前学习自适应一直是做的仿真，这次终于把它运用到了实际的项目中，效果非常显著，更有说服力了。如果大家有什么疑问的可以联系我，我们一起探讨。