守望者--AIR技术交流

标题: Matlab LMS 算法和 RLS 算法实现 [打印本页]

---

作者: 破晓    时间: 2015-8-20 14:24
标题: Matlab LMS 算法和 RLS 算法实现


LMS  算法：
LMS算法的更新迭代公式是：w(n+1) = w(n) + 2ηe * (n) x(n), 算法每次迭代计算仅需要 M + 1 次惩罚和M次加法计算

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %LMS自适应滤波器性能分析
   
3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
4. clear
   
5. clc
   
6. close all
   
7. %%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%%%%%%%%%%%
   
8. w=3;
   
9. % step_len=0.075;%步长
   
10. step_len=0.05;%步长
    
11. 
    
12. variance_v=0.01;%方差
    
13. repeat_times=20;%重复次数
    
14. iteration_times=500;%迭代次数
    
15. filt_len=2;%滤波器长度
    
16. delay=fix(filt_len/2)-1;
    
17. delay1=1;
    
18. 
    
19. 
    
20. 
    
21. %%%%%%%%%%%%%%%%%%%%变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
22. Chan_factor=(1+cos(2*pi*([1:1:3]-2)/w))/2
    
23. Square_err=zeros(repeat_times,iteration_times);
    
24. % Chan_factor1=[-0.195 0.95]
    
25. Chan_factor1=[-0.195 0.95]
    
26. X=zeros(filt_len,1);
    
27. W1=[];
    
28. W2=[];
    
29. Wm1=[repeat_times,iteration_times];
    
30. Wm2=[repeat_times,iteration_times];
    
31. 
    
32. 
    
33. %%%%%%%%%%%%%外循环，重复做repeat_timws次实验，平均值%%%%%%%%%%%%%%%%
    
34. 
    
35. for loop1=1:repeat_times%循环次数
    
36.      
    
37. %%%%%%%%%%%%%%内循环变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%
    
38. Source=zeros(iteration_times,1);
    
39. W=zeros(filt_len,1);
    
40. 
    
41. %%%%%%%%%%%%%内循环，做iteration_times次迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
42.    for i=1:iteration_times
    
43. %%%%%%%%%%%%%%信源，产生等概率分布的正负1%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
44.       if rand<0.5
    
45.           Source(i)=-1;
    
46.       else
    
47.           Source(i)=1;
    
48.       end
    
49.    end
    
50. %%%%%%%%%%%%%%信道，滤波和加噪声%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
51. Revdata=filter(Chan_factor1,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    
52. %  Revdata=filter(Chan_factor,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    
53. %    Revdata=filter(Chan_factor1,1,Source);
    
54.    
    
55. % RT=filter(Chan_factor1,1,Source);
    
56. %%%%%%%%%%%%%%均衡器，基于LMS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    
57.     for loop2=1:iteration_times-filt_len+1
    
58.          
    
59. %         U=Revdata(loop2:loop2+filt_len-1);
    
60. %          y(n)=W'*U=W'*x(n),U=x(n)=输入
    
61. %          d为输出
    
62.         U=Source(loop2:loop2+filt_len-1);
    
63.         d=Revdata(loop2+1);
    
64. %         d=Source(loop2+delay);
    
65. %         X= x(3)=v(3)-a1*x(2)-a2*x(1);
    
66.          
    
67. %           d=RT(loop2+delay1);
    
68. %         Y(loop2)=W'*U;
    
69.         e=d-W'*U;
    
70. %         e=d-W'*U;
    
71.         W=W+step_len*U*e;
    
72.          
    
73.         Wm1(loop2,loop1)=W(1);
    
74.         Wm2(loop2,loop1)=W(2);
    
75. %         W1(loop2)=W(1);
    
76. %         W2(loop2)=W(2);
    
77.         Square_err(loop1,loop2)=e.^2;
    
78.     end
    
79.      
    
80.      
    
81. end
    
82. 
    
83. 
    
84. 
    
85. %%%%%%%%%%%%%%作图，显示结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
86. Square_err_mean=sum(Square_err(:,:))/repeat_times;
    
87. % Wm2
    
88. W1=Wm1(:,repeat_times);
    
89. W2=Wm2(:,repeat_times);
    
90. meanW1=mean(Wm1');
    
91. meanW2=mean(Wm2');
    
92. Square_err;
    
93. %semilogy(Square_err_mean);
    
94. figure(1);
    
95. plot(Square_err_mean);
    
96. hold on;
    
97. plot(Square_err(1,:),'r');
    
98. legend('20次重复试验平均误差','单次误差');
    
99. xlabel('迭代次数');
    
100. ylabel('均方误差');
     
101. title('步长为0.075，500次迭代均方误差变化曲线');
     
102. grid on;
     
103. hold on;
     
104. 
     
105. figure(2);
     
106. plot(W1);
     
107. hold on;
     
108. plot(W2,'r');
     
109. hold on;
     
110. plot(meanW1,'k');
     
111. hold on;
     
112. plot(meanW2,'--k');
     
113. legend('W1','W2','20次重复试验W1','20次重复试验W2');
     
114. xlabel('迭代次数');
     
115. ylabel('抽头权值');
     
116. title('步长为0.075，500次迭代权值学习曲线');
     
117. grid on;

复制代码

RLS 算法

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %RLS自适应滤波器性能分析
   
3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
4. clear all
   
5. %clc
   
6. %close all
   
7. %%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%%%%%%%%%%%
   
8. w=2.9;
   
9. % forget_factor=0.995;
   
10. forget_factor=0.5;
    
11. variance_v=0.001;
    
12. repeat_times=20;
    
13. iteration_times=500;
    
14. filt_len=2;
    
15. delay=fix(filt_len/2)-1;
    
16. %%%%%%%%%%%%%%%%%%%%变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
17. Chan_factor=(1+cos(2*pi*([1:1:3]-2)/w))/2;
    
18. Chan_factor1=[-0.195,0.95]
    
19. Square_err=zeros(repeat_times,iteration_times);
    
20. W1=[];
    
21. W2=[];
    
22. Wm1=[repeat_times,iteration_times];
    
23. Wm2=[repeat_times,iteration_times];
    
24. 
    
25. 
    
26. %%%%%%%%%%%%%外循环，重复做repeat_times次实验，平均值%%%%%%%%%%%%%%%%
    
27. for loop1=1:repeat_times
    
28.      
    
29.      
    
30. %%%%%%%%%%%%%%内循环变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%
    
31. Source=zeros(iteration_times,1);
    
32. W=zeros(filt_len,1);
    
33. P=eye(filt_len);
    
34. Pi=zeros(filt_len,1);
    
35. K=zeros(filt_len,1);
    
36. 
    
37. 
    
38. 
    
39. %%%%%%%%%%%%%内循环，做iteration_times次迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
40.    for i=1:iteration_times
    
41. %%%%%%%%%%%%%%信源，产生等概率分布的正负1%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
42.       if rand<0.5
    
43.           Source(i)=-1;
    
44.       else
    
45.           Source(i)=1;
    
46.       end
    
47.    end
    
48. %%%%%%%%%%%%%%信道，滤波和加噪声%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
49.    Revdata=filter(Chan_factor1,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    
50. %%%%%%%%%%%%%%均衡器，基于RLS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    
51.     for loop2=1:iteration_times-filt_len+1
    
52. 
    
53. %         U=Revdata(loop2:loop2+filt_len-1);
    
54. %         d=Source(loop2+delay);
    
55.         U=Source(loop2:loop2+filt_len-1);
    
56.         d=Revdata(loop2+1);
    
57.         Pi=P*U;
    
58.         K=Pi/(forget_factor+U'*Pi);
    
59.         e=d-W'*U;
    
60.         W=W+K*e;
    
61.         P=P/forget_factor-K*U'*P/forget_factor;
    
62.          
    
63.         Wm1(loop2,loop1)=W(1);
    
64.         Wm2(loop2,loop1)=W(2);
    
65. 
    
66.         Square_err(loop1,loop2)=e.^2;
    
67.     end
    
68. end
    
69. %%%%%%%%%%%%%%作图，显示结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
70. Square_err_mean=sum(Square_err(:,:))/repeat_times;
    
71. 
    
72. 
    
73. W1=Wm1(:,repeat_times)
    
74. W2=Wm2(:,repeat_times)
    
75. meanW1=mean(Wm1');
    
76. meanW2=mean(Wm2');
    
77. %semilogy(Square_err_mean);
    
78. figure(1);
    
79. plot(Square_err_mean);
    
80. hold on;
    
81. plot(Square_err(1,:),'r');
    
82. legend('20次重复试验平均误差','单次误差');
    
83. xlabel('迭代次数');
    
84. ylabel('均方误差');
    
85. title('遗忘因子为0.8，500次迭代均方误差变化曲线');
    
86. grid on;
    
87. hold on;
    
88. 
    
89. plot(Square_err_mean);
    
90. hold off;
    
91. 
    
92. figure(2);
    
93. 
    
94. plot(W1);
    
95. hold on;
    
96. plot(W2,'r');
    
97. hold on;
    
98. plot(meanW1,'k');
    
99. hold on;
    
100. plot(meanW2,'--k');
     
101. legend('W1','W2','20次重复试验W1','20次重复试验W2');
     
102. xlabel('迭代次数');
     
103. ylabel('抽头权值');
     
104. title('遗忘因子为0.8，500次迭代权值学习曲线');
     
105. grid on;

复制代码

---

作者: 破晓    时间: 2015-8-20 18:08

1. #include "math.h"    
2. #include "stdio.h"   
3. int     N=5;                           /* FIR滤波器阶数N */   
4. float   u=0.001;                       /* 自适应步长u */   
5. float   pi=3.1415926;   
6. int     i,j;   
7. float   x[504];                       
8. float   y[500];   
9. float   d[500];   
10. float   e[500];   
11. float   X[5];   
12. float   W[5];   
13. void main()   
14. {      
15.     for (i=0;i<504;i++)                /* 输入信号为正弦 */   
16.     {    
17.         x[i]=sin((9*i-36)*pi/180);     
18.     }   
19.        
20.                                         
21.     for (i=0;i<500;i++)                /* 期望信号等于输入信号 */     
22.     {    
23.         d[i]=x[i+N-1];     
24.     }   
25.        
26.     for (i=0;i<N;i++)                  /* 滤波器权矢量W初始化为0 */   
27.     {      
28.         W[i]=0;   
29.     }   
30.        
31.     for (i=0;i<500;i++)                /* 迭代运算 */   
32.     {      
33.         y[i]=0;   
34.         for(j=0;j<N;j++)            /* 更新X */   
35.         {   
36.             X[j]=x[i+N-j-1];       
37.         }   
38.         for (j=0;j<N;j++)                   
39.         {    
40.             y[i]=y[i]+W[j]*X[j];       /* 计算y */               
41.         }   
42.         e[i]=d[i]-y[i];    
43.         for(j=0;j<N;j++)   
44.         {   
45.             W[j]=W[j]+2*u*e[i]*X[j];   /* 更新W */   
46.         }      
47.     }      
48. }   

---

欢迎光临 守望者--AIR技术交流 (http://www.airmyth.com/)