• IIR滤波器设计(调用MATLAB IIR函数来实现)


    转载请注明文章来源 – http://blog.csdn.net/v_hyx ,请勿用于任何商业用途
     
          对于滤波器设计,以前虽然学过相关的理论(现代数字信号处理和DSP设计),但一直不求甚解,也没用过。趁着最近使用了一下,就来重学一回,温故而知新。
          先来说说IIR滤波器设计,理论与原理参考如下博客,写得简明易懂,不错。
    http://blog.csdn.net/thnh169/article/details/9034483  [数字信号处理]IIR滤波器基础
    http://blog.csdn.net/thnh169/article/details/9069475  [数字信号处理]IIR滤波器的间接设计(C代码)
    http://blog.csdn.net/thnh169/article/details/9076283  [数字信号处理]IIR滤波器的直接设计(C代码)
     
          一般IIR设计可分三种:间接设计(原型转换设计)、直接设计、直接使用工具软件如MATLAB的IIR函数设计。前2种方法,上面的博客已经写得很清楚,理论比较多,设计还是很复杂的。但在实际工程应用中,多采用MATLAB的IIR函数或者FDATOOL工具进行,非常方便快捷。
          OK,来个示例来说明采用MATLAB的IIR函数设计过程,花一会儿的功夫就可以快速入门,so easy!废话不多说,直接上MATLAB IIR.m文件,最后附上效果图。
     
     
    IIR.m
      1 % IIR滤波器设计
      2 % 目的:设计一个采样频率为1000Hz、通带截止频率为50Hz、阻带截止频率为100Hz的低通滤波器,并要求通带最大衰减为1dB,阻带最小衰减为60dB。
      3 
      4 clc;clear;close all;
      5 
      6 % 1. 产生信号(频率为10Hz和100Hz的正弦波叠加)
      7 Fs=1000; % 采样频率1000Hz
      8 t=0:1/Fs:1;
      9 s10=sin(20*pi*t); % 产生10Hz正弦波
     10 s100=sin(200*pi*t); % 产生100Hz正弦波
     11 s=s10+s100; % 信号叠加
     12 
     13 figure(1); % 画图
     14 subplot(2,1,1);plot(s);grid;
     15 title('原始信号');
     16 
     17 % 2. FFT分析信号频谱
     18 len = 512;
     19 y=fft(s,len);  % 对信号做len点FFT变换
     20 f=Fs*(0:len/2 - 1)/len;
     21 
     22 subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
     23 title('原始信号频谱')
     24 xlabel('Hz');ylabel('幅值');
     25 
     26 % 3. IIR滤波器设计
     27 N=0; % 阶数
     28 Fp=50; % 通带截止频率50Hz
     29 Fc=100; % 阻带截止频率100Hz
     30 Rp=1; % 通带波纹最大衰减为1dB
     31 Rs=60; % 阻带衰减为60dB
     32 
     33 % 3.0 计算最小滤波器阶数
     34 na=sqrt(10^(0.1*Rp)-1);
     35 ea=sqrt(10^(0.1*Rs)-1);
     36 N=ceil(log10(ea/na)/log10(Fc/Fp));
     37 
     38 % 3.1 巴特沃斯滤波器
     39 Wn=Fp*2/Fs;
     40 [Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器
     41 [BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线
     42 Bf=filter(Bb,Ba,s); % 进行低通滤波
     43 By=fft(Bf,len);  % 对信号f1做len点FFT变换
     44 
     45 figure(2); % 画图
     46 subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
     47 legend('10Hz原始信号','巴特沃斯滤波器滤波后');
     48 subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
     49 title('巴特沃斯低通滤波后信号频谱');
     50 xlabel('Hz');ylabel('幅值');
     51 
     52 % 3.2 切比雪夫I型滤波器
     53 [C1b C1a]=cheby1(N,Rp,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器
     54 [C1H,C1W]=freqz(C1b,C1a); % 绘制频率响应曲线
     55 C1f=filter(C1b,C1a,s); % 进行低通滤波
     56 C1y=fft(C1f,len);  % 对滤波后信号做len点FFT变换
     57 
     58 figure(3); % 画图
     59 subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;
     60 legend('10Hz原始信号','切比雪夫I型滤波器滤波后');
     61 subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;
     62 title('切比雪夫I型滤波后信号频谱');
     63 xlabel('Hz');ylabel('幅值');
     64 
     65 % 3.3 切比雪夫II型滤波器
     66 [C2b C2a]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器
     67 [C2H,C2W]=freqz(C2b,C2a); % 绘制频率响应曲线
     68 C2f=filter(C2b,C2a,s); % 进行低通滤波
     69 C2y=fft(C2f,len);  % 对滤波后信号做len点FFT变换
     70 
     71 figure(4); % 画图
     72 subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;
     73 legend('10Hz原始信号','切比雪夫II型滤波器滤波后');
     74 subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;
     75 title('切比雪夫II型滤波后信号频谱');
     76 xlabel('Hz');ylabel('幅值');
     77 
     78 % 3.4 椭圆滤波器
     79 [Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器
     80 [EH,EW]=freqz(Eb,Ea); % 绘制频率响应曲线
     81 Ef=filter(Eb,Ea,s); % 进行低通滤波
     82 Ey=fft(Ef,len);  % 对滤波后信号做len点FFT变换
     83 
     84 figure(5); % 画图
     85 subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;
     86 legend('10Hz原始信号','椭圆滤波器滤波后');
     87 subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;
     88 title('椭圆滤波后信号频谱');
     89 xlabel('Hz');ylabel('幅值');
     90 
     91 % 3.5 yulewalk滤波器
     92 fyule=[0 Wn Fc*2/Fs 1]; % 在此进行的是简单设计,实际需要多次仿真取最佳值
     93 myule=[1 1 0 0]; % 在此进行的是简单设计,实际需要多次仿真取最佳值
     94 [Yb Ya]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器
     95 [YH,YW]=freqz(Yb,Ya); % 绘制频率响应曲线
     96 Yf=filter(Yb,Ya,s); % 进行低通滤波
     97 Yy=fft(Yf,len);  % 对滤波后信号做len点FFT变换
     98 
     99 figure(6); % 画图
    100 subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;
    101 legend('10Hz原始信号','yulewalk滤波器滤波后');
    102 subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;
    103 title('yulewalk滤波后信号频谱');
    104 xlabel('Hz');ylabel('幅值');
    105 
    106 % 4. 各个滤波器的幅频响应对比分析
    107 A1 = BW*Fs/(2*pi);
    108 A2 = C1W*Fs/(2*pi);
    109 A3 = C2W*Fs/(2*pi);
    110 A4 = EW*Fs/(2*pi);
    111 A5 = YW*Fs/(2*pi);
    112 
    113 figure(7); % 画图
    114 subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;
    115 xlabel('频率/Hz');
    116 ylabel('频率响应幅度');
    117 legend('butter','cheby1','cheby2','ellip','yulewalk');
     
    效果图

    转载请注明文章来源 – http://blog.csdn.net/v_hyx ,请勿用于任何商业用途
     
    参考资料:
    1. http://www.cnblogs.com/sunev/archive/2011/11/23/2260579.html 基于Matlab的FIR滤波器设计与实现
    2. http://blog.csdn.net/thnh169/article/details/9034483  [数字信号处理]IIR滤波器基础
    3. http://blog.csdn.net/thnh169/article/details/9069475  [数字信号处理]IIR滤波器的间接设计(C代码)
    4. http://blog.csdn.net/thnh169/article/details/9076283  [数字信号处理]IIR滤波器的直接设计(C代码)
  • 相关阅读:
    多线程
    socket
    十一、内置模块
    十、模块与包
    Javaweb学习笔记11 JDBC
    Javaweb学习笔记10 4种范围对象
    Javaweb学习笔记9 application
    Javaweb学习笔记8 session与cookie区别
    Javaweb学习笔记7 session
    Javaweb学习笔记6 jsp页面元素
  • 原文地址:https://www.cnblogs.com/IDoIUnderstand/p/3280716.html
Copyright © 2020-2023  润新知