守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

搜索
热搜: ANE FlasCC 炼金术
查看: 1185|回复: 1

[音频分析] Matlab LMS 算法和 RLS 算法实现

[复制链接]
  • TA的每日心情
    擦汗
    2018-4-10 15:18
  • 签到天数: 447 天

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    562
    贡献
    29
    金币
    51753
    钢镚
    1422

    开源英雄守望者

    发表于 2015-8-20 14:24:13 | 显示全部楼层 |阅读模式
    本帖最后由 破晓 于 2015-8-20 14:46 编辑

    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. variance_v=0.01;%方差
    12. repeat_times=20;%重复次数
    13. iteration_times=500;%迭代次数
    14. filt_len=2;%滤波器长度
    15. delay=fix(filt_len/2)-1;
    16. delay1=1;



    17. %%%%%%%%%%%%%%%%%%%%变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    18. Chan_factor=(1+cos(2*pi*([1:1:3]-2)/w))/2
    19. Square_err=zeros(repeat_times,iteration_times);
    20. % Chan_factor1=[-0.195 0.95]
    21. Chan_factor1=[-0.195 0.95]
    22. X=zeros(filt_len,1);
    23. W1=[];
    24. W2=[];
    25. Wm1=[repeat_times,iteration_times];
    26. Wm2=[repeat_times,iteration_times];


    27. %%%%%%%%%%%%%外循环,重复做repeat_timws次实验,平均值%%%%%%%%%%%%%%%%

    28. for loop1=1:repeat_times%循环次数
    29.      
    30. %%%%%%%%%%%%%%内循环变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%
    31. Source=zeros(iteration_times,1);
    32. W=zeros(filt_len,1);

    33. %%%%%%%%%%%%%内循环,做iteration_times次迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    34.    for i=1:iteration_times
    35. %%%%%%%%%%%%%%信源,产生等概率分布的正负1%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    36.       if rand<0.5
    37.           Source(i)=-1;
    38.       else
    39.           Source(i)=1;
    40.       end
    41.    end
    42. %%%%%%%%%%%%%%信道,滤波和加噪声%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    43. Revdata=filter(Chan_factor1,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    44. %  Revdata=filter(Chan_factor,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    45. %    Revdata=filter(Chan_factor1,1,Source);
    46.    
    47. % RT=filter(Chan_factor1,1,Source);
    48. %%%%%%%%%%%%%%均衡器,基于LMS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    49.     for loop2=1:iteration_times-filt_len+1
    50.          
    51. %         U=Revdata(loop2:loop2+filt_len-1);
    52. %          y(n)=W'*U=W'*x(n),U=x(n)=输入
    53. %          d为输出
    54.         U=Source(loop2:loop2+filt_len-1);
    55.         d=Revdata(loop2+1);
    56. %         d=Source(loop2+delay);
    57. %         X= x(3)=v(3)-a1*x(2)-a2*x(1);
    58.          
    59. %           d=RT(loop2+delay1);
    60. %         Y(loop2)=W'*U;
    61.         e=d-W'*U;
    62. %         e=d-W'*U;
    63.         W=W+step_len*U*e;
    64.          
    65.         Wm1(loop2,loop1)=W(1);
    66.         Wm2(loop2,loop1)=W(2);
    67. %         W1(loop2)=W(1);
    68. %         W2(loop2)=W(2);
    69.         Square_err(loop1,loop2)=e.^2;
    70.     end
    71.      
    72.      
    73. end



    74. %%%%%%%%%%%%%%作图,显示结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    75. Square_err_mean=sum(Square_err(:,:))/repeat_times;
    76. % Wm2
    77. W1=Wm1(:,repeat_times);
    78. W2=Wm2(:,repeat_times);
    79. meanW1=mean(Wm1');
    80. meanW2=mean(Wm2');
    81. Square_err;
    82. %semilogy(Square_err_mean);
    83. figure(1);
    84. plot(Square_err_mean);
    85. hold on;
    86. plot(Square_err(1,:),'r');
    87. legend('20次重复试验平均误差','单次误差');
    88. xlabel('迭代次数');
    89. ylabel('均方误差');
    90. title('步长为0.075,500次迭代均方误差变化曲线');
    91. grid on;
    92. hold on;

    93. figure(2);
    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.075,500次迭代权值学习曲线');
    105. 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. %%%%%%%%%%%%%外循环,重复做repeat_times次实验,平均值%%%%%%%%%%%%%%%%
    25. for loop1=1:repeat_times
    26.      
    27.      
    28. %%%%%%%%%%%%%%内循环变量及数组初始化%%%%%%%%%%%%%%%%%%%%%%%%
    29. Source=zeros(iteration_times,1);
    30. W=zeros(filt_len,1);
    31. P=eye(filt_len);
    32. Pi=zeros(filt_len,1);
    33. K=zeros(filt_len,1);



    34. %%%%%%%%%%%%%内循环,做iteration_times次迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    35.    for i=1:iteration_times
    36. %%%%%%%%%%%%%%信源,产生等概率分布的正负1%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    37.       if rand<0.5
    38.           Source(i)=-1;
    39.       else
    40.           Source(i)=1;
    41.       end
    42.    end
    43. %%%%%%%%%%%%%%信道,滤波和加噪声%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    44.    Revdata=filter(Chan_factor1,1,Source)+randn(iteration_times,1)*sqrt(variance_v);
    45. %%%%%%%%%%%%%%均衡器,基于RLS算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    46.     for loop2=1:iteration_times-filt_len+1

    47. %         U=Revdata(loop2:loop2+filt_len-1);
    48. %         d=Source(loop2+delay);
    49.         U=Source(loop2:loop2+filt_len-1);
    50.         d=Revdata(loop2+1);
    51.         Pi=P*U;
    52.         K=Pi/(forget_factor+U'*Pi);
    53.         e=d-W'*U;
    54.         W=W+K*e;
    55.         P=P/forget_factor-K*U'*P/forget_factor;
    56.          
    57.         Wm1(loop2,loop1)=W(1);
    58.         Wm2(loop2,loop1)=W(2);

    59.         Square_err(loop1,loop2)=e.^2;
    60.     end
    61. end
    62. %%%%%%%%%%%%%%作图,显示结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    63. Square_err_mean=sum(Square_err(:,:))/repeat_times;


    64. W1=Wm1(:,repeat_times)
    65. W2=Wm2(:,repeat_times)
    66. meanW1=mean(Wm1');
    67. meanW2=mean(Wm2');
    68. %semilogy(Square_err_mean);
    69. figure(1);
    70. plot(Square_err_mean);
    71. hold on;
    72. plot(Square_err(1,:),'r');
    73. legend('20次重复试验平均误差','单次误差');
    74. xlabel('迭代次数');
    75. ylabel('均方误差');
    76. title('遗忘因子为0.8,500次迭代均方误差变化曲线');
    77. grid on;
    78. hold on;

    79. plot(Square_err_mean);
    80. hold off;

    81. figure(2);

    82. plot(W1);
    83. hold on;
    84. plot(W2,'r');
    85. hold on;
    86. plot(meanW1,'k');
    87. hold on;
    88. plot(meanW2,'--k');
    89. legend('W1','W2','20次重复试验W1','20次重复试验W2');
    90. xlabel('迭代次数');
    91. ylabel('抽头权值');
    92. title('遗忘因子为0.8,500次迭代权值学习曲线');
    93. grid on;
    复制代码


    守望者AIR技术交流社区(www.airmyth.com)
    回复

    使用道具 举报

  • TA的每日心情
    擦汗
    2018-4-10 15:18
  • 签到天数: 447 天

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    562
    贡献
    29
    金币
    51753
    钢镚
    1422

    开源英雄守望者

     楼主| 发表于 2015-8-20 18:08:30 | 显示全部楼层
    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技术交流社区(www.airmyth.com)
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    
    关闭

    站长推荐上一条 /4 下一条

    QQ|手机版|Archiver|网站地图|小黑屋|守望者 ( 京ICP备14061876号

    GMT+8, 2019-10-19 16:01 , Processed in 0.046446 second(s), 32 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

    快速回复 返回顶部 返回列表