自适应滤波器之LMS算法实现Demo

东方耀

自适应滤波器之LMS算法实现Demo

1. # -*- coding: utf-8 -*-
   
2. __author__ = u'东方耀 微信：dfy_88888'
   
3. __date__ = '2021/4/15 11:27'
   
4. __product__ = 'PyCharm'
   
5. __filename__ = '01_自适应LMS滤波算法'
   
6. 
   
7. import numpy as np
   
8. import matplotlib.pyplot as plt
   
9. 
   
10. # 自适应滤波器之LMS算法实现Demo
    
11. 
    
12. # 定义向量的内积
    
13. # np.convolve()
    
14. # np.dot()
    
15. # def multiVector(A, B):
    
16. #     # 两个向量的内积 对应元素相乘后求和  np.dot()
    
17. #     C = sc.zeros(len(A))
    
18. #
    
19. #     for i in range(len(A)):
    
20. #         C[i] = A[i] * B[i]
    
21. #
    
22. #     return sum(C)
    
23. 
    
24. 
    
25. def inVector(vector, b, a):
    
26.     '''
    
27.     从一个长向量vector中按索引截取并返回倒着的序列
    
28.     :param vector:
    
29.     :param b: 终点索引
    
30.     :param a: 起点索引
    
31.     :return:
    
32.     '''
    
33.     sub_vector = np.zeros(shape=(b-a+1, ))
    
34. 
    
35.     for i in range(b - a + 1):
    
36.         sub_vector[i] = vector[i + a]
    
37. 
    
38.     return sub_vector[::-1]
    
39. 
    
40. 
    
41. # lMS算法的函数
    
42. 
    
43. def LMS(xn, dn, M, learning_rate, itr):
    
44.     '''
    
45.     N阶线性系统出发
    
46.     自适应滤波之LMS算法实现
    
47.     :param xn: 输入原始信号或序列
    
48.     :param dn: 期望的
    
49.     :param M: 滤波器阶数=64
    
50.     :param learning_rate: 学习率
    
51.     :param itr: 采样点数
    
52.     :return:返回的yn与en是 输出和误差
    
53.     '''
    
54.     en = np.zeros(shape=(itr, ))   # 误差 或 损失 信号
    
55.     # 空间复杂度有点高 需要优化
    
56.     W = [[0] * M for i in range(itr)]  # 滤波器的参数 初始化 列表长度=采样点数 每个元素有滤波器阶数个值
    
57.     # print(len(W), len(W[0]))   10000 64
    
58. 
    
59.     for k in range(itr)[M - 1:]:
    
60.         # 为什么要从M-1开始切片：63到最后 为根据前面的信号样点作为期望值
    
61.         # 开始的：k=63  k-M+1=0
    
62.         x = inVector(xn, k, k - M + 1)   # len(x) = 64 滤波器的阶数=滤波器长度 决定了输入信号样点的长度
    
63.         # 那个LMS函数里的d是滑动平均作为期望  问题：如果噪声均值不为0 效果会差吧？
    
64.         # 输入原始信号序列还是xn 但是filter是逐步滑动去滤的
    
65.         d = x.mean()
    
66.         # y是一次滤波操作的输出： W[k - 1]是上一次的滤波器的参数  x呢？不是xn吗？
    
67.         # x是前面的 输入样点的 倒序   因为向量点乘要求：长度一致啊
    
68.         # y = multiVector(W[k - 1], x)  一维滤波
    
69.         y = np.dot(W[k - 1], x)
    
70. 
    
71.         # print("第一次应该为0：", y)
    
72.         # assert 0 == 1, "停"
    
73.         # 误差信号是从k=63开始的
    
74.         en[k] = d - y   # 误差如果很小了 可以不需要迭代了
    
75.         # 循环里 每输入一个序列就更新一次滤波器参数  更新频率很高
    
76.         # 梯度= -2 * learning_rate * en[k] * x
    
77.         # 标准时域LMS算法的更新公式  最终发现：误差越来越小
    
78.         W[k] = np.add(W[k - 1], 2 * learning_rate * en[k] * x)
    
79. 
    
80.     # 求最优时滤波器的输出序列
    
81. 
    
82.     yn = np.inf * np.ones(len(xn))  # inf 是无穷大的意思
    
83.     # 先在原始输入序列信号上 不断改进自适应滤波的权重参数
    
84.     filter_best = W[len(W) - 1]
    
85.     # 最后根据得到的最优滤波器，计算原始输入信号的输出
    
86. 
    
87.     for k in range(itr)[M - 1:]:
    
88.         x = inVector(xn, k, k - M + 1)
    
89.         # x才是真正的输入信号序列，错了，x是滤波操作的输入
    
90. 
    
91.         # yn[k] = multiVector(W[len(W) - 1], x)
    
92.         # 默认：最后的滤波器参数就是最优的
    
93.         yn[k] = np.dot(filter_best, x)
    
94. 
    
95.     return (yn, en)
    
96. 
    
97. 
    
98. if __name__ == "__main__":
    
99.     # 参数初始化
    
100. 
     
101.     itr = 10000  # 采样的点数
     
102.     # sigma = 0.12
     
103.     noise_size = itr
     
104. 
     
105.     X = np.linspace(0, 4 * np.pi, itr, endpoint=True)
     
106.     Y = np.sin(X)
     
107.     signal_array = Y  # [0.0]*noise_size  真实的信号
     
108.     noise_array = np.random.normal(0, 0.3, noise_size)  # 加上高斯白噪声
     
109.     # 高斯白噪声是分析信道加性噪声的理想模型，通信中的主要噪声源——热噪声就属于这类噪声
     
110.     """noise_array = []
     
111.     for x in range(itr):
     
112.         noise_array.append(random.gauss(learning_rate,sigma))"""
     
113.     signal_noise_array = signal_array + noise_array   # 信号+噪声
     
114. 
     
115.     M = 64  # 滤波器的阶数
     
116. 
     
117.     learning_rate = 1e-4
     
118. 
     
119.     xn = signal_noise_array  # 原始输入端的信号为被噪声污染的正弦信号
     
120.     # xn并非输入序列，因为输入序列必须与filter长度一致
     
121. 
     
122.     dn = signal_array  # 对于自适应对消器，用dn作为期望 实际没有用这个
     
123. 
     
124.     # 调用LMS算法
     
125. 
     
126.     (yn, en) = LMS(xn, dn, M, learning_rate, itr)
     
127. 
     
128.     # 画出图形
     
129. 
     
130.     plt.figure(1)
     
131. 
     
132.     plt.plot(xn, label="$signal_noise[        DISCUZ_CODE_0        ]quot;)
     
133. 
     
134.     plt.plot(dn, label="$signal[        DISCUZ_CODE_0        ]quot;)
     
135. 
     
136.     plt.xlabel("Time(s)")
     
137. 
     
138.     plt.ylabel("Volt")
     
139. 
     
140.     plt.title("original signal xn and desired signal dn")
     
141. 
     
142.     plt.legend()
     
143. 
     
144.     plt.figure(2)
     
145.     # plt.plot(xn,label="$xn[        DISCUZ_CODE_0        ]quot;)
     
146.     # plt.plot(xn,label="$xn[        DISCUZ_CODE_0        ]quot;)
     
147.     plt.plot(signal_array, label="$signal_array[        DISCUZ_CODE_0        ]quot;)
     
148.     plt.plot(yn, label="$yn[        DISCUZ_CODE_0        ]quot;)
     
149. 
     
150.     plt.xlabel("Time(s)")
     
151. 
     
152.     plt.ylabel("Volt")
     
153. 
     
154.     plt.title("out signal yn best for filter")
     
155. 
     
156.     plt.legend()
     
157. 
     
158.     plt.figure(3)
     
159. 
     
160.     plt.plot(en, label="$error[        DISCUZ_CODE_0        ]quot;)
     
161. 
     
162.     plt.xlabel("Time(s)")
     
163. 
     
164.     plt.ylabel("Volt")
     
165. 
     
166.     plt.title("error between y and d")
     
167. 
     
168.     plt.legend()
     
169. 
     
170.     plt.show()
     

复制代码

东方耀

LMS算法步长调整原则：滤波器的收敛速度和稳态误差是一对矛盾。当拥有较大的步长因子会加快滤波器的收敛速度，但同时也会增大稳态失调误差。较小的步长因子会减小滤波器的收敛速度但会拥有较小的稳态误差

东方耀

东方耀 发表于 2021-4-16 09:20
LMS算法步长调整原则：滤波器的收敛速度和稳态误差是一对矛盾。当拥有较大的步长因子会加快滤波器的收敛速 ...

为了克服这一矛盾，在现有固定步长算法的基础上我们提出了新的算法—变步长自适应算法。其优点是步长值随误差的改变而改变。在误差较大的时候步长因子也大以获得更快的收敛速度；当收敛过程接近稳态，误差小的时候其步长因子也小以获得较小的稳态误差