《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 156|回复: 0
收起左侧

[源码] 基于LMS算法的自适应Notch滤波器

[复制链接]
发表于 2020-2-6 19:23:05 | 显示全部楼层 |阅读模式
产生中心频率为15khz的CW信号,然后加入带限白噪声(频带范围5kHz~15kHz),构建基于LMS算法的单通道自适应notch滤波器,对上述信号进行滤波,最终给出滤波输出信号的包络、相位及瞬时频率。

  1. %产生cw信号
  2. f0=15e3;
  3. fs=60e3;
  4. ts=1/fs;
  5. SNR=1;
  6. t=0:ts:1-ts;
  7. s=0.5*cos(2*pi*f0*t)+0.5*sin(2*pi*f0*t);
  8. cw=[zeros(1,length(s)) s zeros(1,length(s))];

  9. %产生带限白噪声
  10. N=length(cw);
  11. t1=(1:length(cw))*ts;
  12. noise=randn(1,length(cw));
  13. b=fir1(256,[5e3,15e3]/(fs/2));          %带通滤波器
  14. noise1=filter(b,1,noise);                  %白噪声通过滤波器
  15. noise2=noise1/std(noise1);              %归一化
  16. a=0.5/(10^(SNR/10));                         %根据信噪比计算噪声参数
  17. noise3=noise2*a;

  18. %叠加信号
  19. d=cw+noise3;

  20. %notch滤波器
  21. w1(1)=0.1;                       %权矢量初值
  22. w2(1)=0.1;                     
  23. e=zeros(1,N);
  24. y=zeros(1,N);
  25. u=0.001;                          %步长
  26. x1=cos(2*pi*(f0)*t1);
  27. x2=sin(2*pi*(f0)*t1);
  28. for i=1:N                           %LMS算法
  29.     y(i)=w1(i)*x1(i)+w2(i)*x2(i);
  30.     e(i)=d(i)-y(i);
  31.     w1(i+1)=w1(i)+2*u*e(i)*x1(i);
  32.     w2(i+1)=w2(i)+2*u*e(i)*x2(i);
  33.    
  34. end

  35. %画图
  36. f=-fs/2:fs/length(cw):fs/2-fs/length(cw);
  37. figure(1)
  38. plot(t1,cw);
  39. title('原信号时域波形 ');
  40. xlabel('时间')
  41. ylabel('幅度')

  42. figure(2)
  43. subplot(2,1,1)
  44. plot(t1,d);
  45. title('叠加信号时域波形 ');
  46. xlabel('时间')
  47. ylabel('幅度')
  48. subplot(2,1,2)
  49. plot(t1,y);
  50. title('滤波输出信号时域波形 ');
  51. xlabel('时间')
  52. ylabel('幅度')

  53. figure(3)
  54. subplot(2,1,1)
  55. plot(f,fftshift(abs(fft(cw))),'r');
  56. title('原信号频谱图');
  57. xlabel('频率')
  58. ylabel('幅度')
  59. subplot(2,1,2)
  60. plot(f,fftshift(abs(fft(d))),'b');
  61. title('叠加信号频谱图');
  62. xlabel('频率')
  63. ylabel('幅度')


  64. figure(4)
  65. plot(t1,e);
  66. title('残差时域波形 ');
  67. xlabel('时间')
  68. ylabel('幅度')

  69. figure(5)
  70. subplot(2,1,1)
  71. plot(t1,w1(1:N));
  72. title('权w1时域波形 ');
  73. xlabel('时间')
  74. ylabel('幅度')
  75. subplot(2,1,2)
  76. plot(t1,w2(1:N));
  77. title('权w2时域波形 ');
  78. xlabel('时间')
  79. ylabel('幅度')

  80. figure(6)
  81. subplot(2,1,1)
  82. plot(t1,abs(hilbert(y)));
  83. title('滤波输出信号包络 ');
  84. xlabel('时间')
  85. ylabel('幅度')
  86. subplot(2,1,2)
  87. plot(f,fftshift(abs(fft(y))));
  88. title('滤波输出信号频谱图 ');
  89. xlabel('频率')
  90. ylabel('幅度')

  91. xiangwei=atan(w2(1:N)./w1(1:N));
  92. w11=2*pi*f0-(atan(w2(2:end)./w1(2:end))-atan(w2(1:end-1)./w1(1:end-1)));
  93. figure(7)
  94. subplot(2,1,1)
  95. plot(t1,xiangwei);
  96. title('滤波输出信号相位 ');
  97. xlabel('时间')
  98. ylabel('相位')
  99. subplot(2,1,2)
  100. plot(t1,w11);
  101. title('滤波输出信号瞬时频率 ');
  102. xlabel('时间')
  103. ylabel('频率')
复制代码
您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 蜀ICP备19014457号 

GMT+8, 2020-7-12 00:33 , Processed in 0.052255 second(s), 11 queries , Gzip On, MemCached On.

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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