• 使用MATLAB验证IIR滤波参数

    　　基础使用教程：https://ww2.mathworks.cn/help/matlab/getting-started-with-matlab.html?s_tid=CRUX_lftnav 

    　　界面说明在文章末尾，太久没使用的下去看下我们直接开干 转载请注明文章来源 ，请勿用于任何商业用途

          对于滤波器设计，以前虽然学过相关的理论（现代数字信号处理和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，fdatool调出工具计算系数；来个示例来说明采用MATLAB的IIR函数设计过程，花一会儿的功夫就可以快速入门，so easy！废话不多说，直接上MATLAB IIR.m文件，最后附上效果图。
    % IIR滤波器设计  
% 目的：设计一个采样频率为500Hz、通带截止频率为5Hz、阻带截止频率为10Hz的低通滤波器，并要求通带最大衰减为1dB，阻带最小衰减为60dB。  
% 1. 产生信号（频率为1Hz和10Hz的正弦波叠加）  
Fs=500; % 采样频率500Hz  
t=0:1/Fs:1;  
sa=sin(5*pi*t); % 产生5Hz正弦波  
gssa = awgn(sa, 1);
sb=sin(10*pi*t); % 产生10Hz正弦波
gssb = awgn(sb, 1); 
s=gssa+gssb; % 信号叠加  
figure(1); % 画图  
subplot(2,1,1);
plot(s);
grid;
title('原始信号');xlabel('x');ylabel('Hz'); 


% 2. FFT分析信号频谱  
len = 500;
y=fft(s,len);  % 对信号做len点FFT变换  
f=Fs*(0:len/2 - 1)/len;  
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;  
title('原始信号频谱');xlabel('Hz');ylabel('幅值');  


% 3. IIR滤波器设计  
N=0; % 阶数  
Fp=5; % 通带截止频率50Hz  
Fc=10; % 阻带截止频率100Hz  
Rp=1; % 通带波纹最大衰减为1dB  
Rs=60; % 阻带衰减为60dB  
  
% 3.0 计算最小滤波器阶数  
na=sqrt(10^(0.1*Rp)-1);  
ea=sqrt(10^(0.1*Rs)-1);  
N=ceil(log10(ea/na)/log10(Fc/Fp));  

% 3.1 巴特沃斯滤波器  
Wn=Fp*2/Fs;  
[Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器  
[BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线  
Bf=filter(Bb,Ba,s); % 进行低通滤波  
By=fft(Bf,len);  % 对信号f1做len点FFT变换  
  
figure(2); % 画图  
subplot(2,1,1);plot(t,s,'blue',t,Bf,'red');grid;  
legend('10Hz原始信号','巴特沃斯滤波器滤波后');  
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;  
title('巴特沃斯低通滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.2 切比雪夫I型滤波器  
[C1b C1a]=cheby1(N,Rs,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器  
[C1H,C1W]=freqz(C1b,C1a); % 绘制频率响应曲线  
C1f=filter(C1b,C1a,s); % 进行低通滤波  
C1y=fft(C1f,len);  % 对滤波后信号做len点FFT变换  
  
figure(3); % 画图  
subplot(2,1,1);plot(t,s,'blue',t,C1f,'red');grid;  
legend('10Hz原始信号','切比雪夫I型滤波器滤波后');  
subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;  
title('切比雪夫I型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.3 切比雪夫II型滤波器  
[C2b C2a]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器  
[C2H,C2W]=freqz(C2b,C2a); % 绘制频率响应曲线  
C2f=filter(C2b,C2a,s); % 进行低通滤波  
C2y=fft(C2f,len);  % 对滤波后信号做len点FFT变换  
  
figure(4); % 画图  
subplot(2,1,1);
plot(t,s,'blue',t,C2f,'red');grid;  
legend('10Hz原始信号','切比雪夫II型滤波器滤波后');  
subplot(2,1,2);
plot(f,abs(C2y(1:len/2)));grid;  
title('切比雪夫II型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.4 椭圆滤波器  
[Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器  
[EH,EW]=freqz(Eb,Ea); % 绘制频率响应曲线  
Ef=filter(Eb,Ea,s); % 进行低通滤波  
Ey=fft(Ef,len);  % 对滤波后信号做len点FFT变换  
  
figure(5); % 画图  
subplot(2,1,1);plot(t,s,'blue',t,Ef,'red');grid;  
legend('10Hz原始信号','椭圆滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;  
title('椭圆滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.5 yulewalk滤波器  
fyule=[0 Wn Fc*2/Fs 1]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
myule=[1 1 0 0]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
[Yb Ya]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器  
[YH,YW]=freqz(Yb,Ya); % 绘制频率响应曲线  
Yf=filter(Yb,Ya,s); % 进行低通滤波  
Yy=fft(Yf,len);  % 对滤波后信号做len点FFT变换  
  
figure(6); % 画图  
subplot(2,1,1);plot(t,s,'blue',t,Yf,'red');grid;  
legend('10Hz原始信号','yulewalk滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;  
title('yulewalk滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 4. 各个滤波器的幅频响应对比分析  
A1 = BW*Fs/(2*pi);  
A2 = C1W*Fs/(2*pi);  
A3 = C2W*Fs/(2*pi);  
A4 = EW*Fs/(2*pi);  
A5 = YW*Fs/(2*pi);  
  
figure(7); % 画图  
subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;  
xlabel('频率／Hz');  
ylabel('频率响应幅度');  
legend('butter','cheby1','cheby2','ellip','yulewalk');
    0-20HZ
    % IIR滤波器设计  
% 目的：设计一个采样频率为1000Hz、通带截止频率为50Hz、阻带截止频率为100Hz的低通滤波器，并要求通带最大衰减为1dB，阻带最小衰减为60dB。  
  
clc;clear;close all;  
  
% 1. 产生信号（频率为10Hz和100Hz的正弦波叠加）  
Fs=1000; % 采样频率1000Hz  
t=0:1/Fs:1;  
s10=sin(20*pi*t); % 产生10Hz正弦波  
s100=sin(200*pi*t); % 产生100Hz正弦波  
s=s10+s100; % 信号叠加  
  
figure(1); % 画图  
subplot(2,1,1);plot(s);grid;  
title('原始信号');  
  
% 2. FFT分析信号频谱  
len = 512;  
y=fft(s,len);  % 对信号做len点FFT变换  
f=Fs*(0:len/2 - 1)/len;  
  
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;  
title('原始信号频谱')  
xlabel('Hz');ylabel('幅值');  
  
% 3. IIR滤波器设计  
N=0; % 阶数  
Fp=50; % 通带截止频率50Hz  
Fc=100; % 阻带截止频率100Hz  
Rp=1; % 通带波纹最大衰减为1dB  
Rs=60; % 阻带衰减为60dB  
  
% 3.0 计算最小滤波器阶数  
na=sqrt(10^(0.1*Rp)-1);  
ea=sqrt(10^(0.1*Rs)-1);  
N=ceil(log10(ea/na)/log10(Fc/Fp));  
  
% 3.1 巴特沃斯滤波器  
Wn=Fp*2/Fs;  
[Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器  
[BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线  
Bf=filter(Bb,Ba,s); % 进行低通滤波  
By=fft(Bf,len);  % 对信号f1做len点FFT变换  
  
figure(2); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;  
legend('10Hz原始信号','巴特沃斯滤波器滤波后');  
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;  
title('巴特沃斯低通滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.2 切比雪夫I型滤波器  
[C1b C1a]=cheby1(N,Rp,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器  
[C1H,C1W]=freqz(C1b,C1a); % 绘制频率响应曲线  
C1f=filter(C1b,C1a,s); % 进行低通滤波  
C1y=fft(C1f,len);  % 对滤波后信号做len点FFT变换  
  
figure(3); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;  
legend('10Hz原始信号','切比雪夫I型滤波器滤波后');  
subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;  
title('切比雪夫I型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.3 切比雪夫II型滤波器  
[C2b C2a]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器  
[C2H,C2W]=freqz(C2b,C2a); % 绘制频率响应曲线  
C2f=filter(C2b,C2a,s); % 进行低通滤波  
C2y=fft(C2f,len);  % 对滤波后信号做len点FFT变换  
  
figure(4); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;  
legend('10Hz原始信号','切比雪夫II型滤波器滤波后');  
subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;  
title('切比雪夫II型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.4 椭圆滤波器  
[Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器  
[EH,EW]=freqz(Eb,Ea); % 绘制频率响应曲线  
Ef=filter(Eb,Ea,s); % 进行低通滤波  
Ey=fft(Ef,len);  % 对滤波后信号做len点FFT变换  
  
figure(5); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;  
legend('10Hz原始信号','椭圆滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;  
title('椭圆滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.5 yulewalk滤波器  
fyule=[0 Wn Fc*2/Fs 1]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
myule=[1 1 0 0]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
[Yb Ya]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器  
[YH,YW]=freqz(Yb,Ya); % 绘制频率响应曲线  
Yf=filter(Yb,Ya,s); % 进行低通滤波  
Yy=fft(Yf,len);  % 对滤波后信号做len点FFT变换  
  
figure(6); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;  
legend('10Hz原始信号','yulewalk滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;  
title('yulewalk滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 4. 各个滤波器的幅频响应对比分析  
A1 = BW*Fs/(2*pi);  
A2 = C1W*Fs/(2*pi);  
A3 = C2W*Fs/(2*pi);  
A4 = EW*Fs/(2*pi);  
A5 = YW*Fs/(2*pi);  
  
figure(7); % 画图  
subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;  
xlabel('频率／Hz');  
ylabel('频率响应幅度');  
legend('butter','cheby1','cheby2','ellip','yulewalk');
    View Code
    % IIR滤波器设计  
% 目的：设计一个采样频率为1000Hz、通带截止频率为50Hz、阻带截止频率为100Hz的低通滤波器，并要求通带最大衰减为1dB，阻带最小衰减为60dB。  
  
clc;clear;close all;  
  
% 1. 产生信号（频率为10Hz和100Hz的正弦波叠加）  
Fs=1000; % 采样频率1000Hz  
t=0:1/Fs:1;  
s10=sin(20*pi*t); % 产生10Hz正弦波  
s100=sin(200*pi*t); % 产生100Hz正弦波  
s=s10+s100; % 信号叠加  
  
figure(1); % 画图  
subplot(2,1,1);plot(s);grid;  
title('原始信号');  
  
% 2. FFT分析信号频谱  
len = 512;  
y=fft(s,len);  % 对信号做len点FFT变换  
f=Fs*(0:len/2 - 1)/len;  
  
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;  
title('原始信号频谱')  
xlabel('Hz');ylabel('幅值');  
  
% 3. IIR滤波器设计  
N=0; % 阶数  
Fp=50; % 通带截止频率50Hz  
Fc=100; % 阻带截止频率100Hz  
Rp=1; % 通带波纹最大衰减为1dB  
Rs=60; % 阻带衰减为60dB  
  
% 3.0 计算最小滤波器阶数  
na=sqrt(10^(0.1*Rp)-1);  
ea=sqrt(10^(0.1*Rs)-1);  
N=ceil(log10(ea/na)/log10(Fc/Fp));  
  
% 3.1 巴特沃斯滤波器  
Wn=Fp*2/Fs;  
[Bb Ba]=butter(N,Wn,'low'); % 调用MATLAB butter函数快速设计滤波器  
[BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线  
Bf=filter(Bb,Ba,s); % 进行低通滤波  
By=fft(Bf,len);  % 对信号f1做len点FFT变换  
  
figure(2); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;  
legend('10Hz原始信号','巴特沃斯滤波器滤波后');  
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;  
title('巴特沃斯低通滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.2 切比雪夫I型滤波器  
[C1b C1a]=cheby1(N,Rp,Wn,'low'); % 调用MATLAB cheby1函数快速设计低通滤波器  
[C1H,C1W]=freqz(C1b,C1a); % 绘制频率响应曲线  
C1f=filter(C1b,C1a,s); % 进行低通滤波  
C1y=fft(C1f,len);  % 对滤波后信号做len点FFT变换  
  
figure(3); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,C1f,'red');grid;  
legend('10Hz原始信号','切比雪夫I型滤波器滤波后');  
subplot(2,1,2);plot(f,abs(C1y(1:len/2)));grid;  
title('切比雪夫I型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.3 切比雪夫II型滤波器  
[C2b C2a]=cheby2(N,Rs,Wn,'low'); % 调用MATLAB cheby2函数快速设计低通滤波器  
[C2H,C2W]=freqz(C2b,C2a); % 绘制频率响应曲线  
C2f=filter(C2b,C2a,s); % 进行低通滤波  
C2y=fft(C2f,len);  % 对滤波后信号做len点FFT变换  
  
figure(4); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,C2f,'red');grid;  
legend('10Hz原始信号','切比雪夫II型滤波器滤波后');  
subplot(2,1,2);plot(f,abs(C2y(1:len/2)));grid;  
title('切比雪夫II型滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.4 椭圆滤波器  
[Eb Ea]=ellip(N,Rp,Rs,Wn,'low'); % 调用MATLAB ellip函数快速设计低通滤波器  
[EH,EW]=freqz(Eb,Ea); % 绘制频率响应曲线  
Ef=filter(Eb,Ea,s); % 进行低通滤波  
Ey=fft(Ef,len);  % 对滤波后信号做len点FFT变换  
  
figure(5); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Ef,'red');grid;  
legend('10Hz原始信号','椭圆滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Ey(1:len/2)));grid;  
title('椭圆滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 3.5 yulewalk滤波器  
fyule=[0 Wn Fc*2/Fs 1]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
myule=[1 1 0 0]; % 在此进行的是简单设计，实际需要多次仿真取最佳值  
[Yb Ya]=yulewalk(N,fyule,myule); % 调用MATLAB yulewalk函数快速设计低通滤波器  
[YH,YW]=freqz(Yb,Ya); % 绘制频率响应曲线  
Yf=filter(Yb,Ya,s); % 进行低通滤波  
Yy=fft(Yf,len);  % 对滤波后信号做len点FFT变换  
  
figure(6); % 画图  
subplot(2,1,1);plot(t,s10,'blue',t,Yf,'red');grid;  
legend('10Hz原始信号','yulewalk滤波器滤波后');  
subplot(2,1,2);plot(f,abs(Yy(1:len/2)));grid;  
title('yulewalk滤波后信号频谱');  
xlabel('Hz');ylabel('幅值');  
  
% 4. 各个滤波器的幅频响应对比分析  
A1 = BW*Fs/(2*pi);  
A2 = C1W*Fs/(2*pi);  
A3 = C2W*Fs/(2*pi);  
A4 = EW*Fs/(2*pi);  
A5 = YW*Fs/(2*pi);  
  
figure(7); % 画图  
subplot(1,1,1);plot(A1,abs(BH),A2,abs(C1H),A3,abs(C2H),A4,abs(EH),A5,abs(YH));grid;  
xlabel('频率／Hz');  
ylabel('频率响应幅度');  
legend('butter','cheby1','cheby2','ellip','yulewalk');
    View Code
    1. 效果图

    转载请注明文章来源 – 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代码)

          

    　　界面说明

        

       

    简单的 变量创建和函数运行：

    简单的运算

     查看源码

    出错解决

    在matlab命令行中键入“mex -setup”，系统报错，显示：

    “错误使用 mex
    未找到支持的编译器或 SDK。您可以安装免费提供的 MinGW-w64 C/C++ 编译器；请参阅安装 MinGW-w64 编译
    器。有关更多选项，请访问 http://www.mathworks.com/support/compilers/R2016b/win64.html。”

    1、下载：MinGW-w64 C/C++
    http://tdm-gcc.tdragon.net/download

    2、安装MinGW-w64 C/C++ 编译器
    注意：按照默认路径安装

     3、在Windows 上将MW_MINGW64设置为系统环境变量
    打开计算机的"控制面板",打开”控制面板“中的”系统与安全“中的”系统“。也就是系统属性。

    点击高级系统设置>高级选项卡。点击环境变量。

    在系统变量下，选择新建。

    在“新建系统变量”对话框的“变量名称”字段中输入MW_MINGW64_LOC。
    在变量值字段中，键入MinGW-w64编译器安装的位置，例如“C:TDM-GCC-64”（注意：“C:TDM-GCC-64”是MinGW的安装目录）

    点击“确定”关闭对话框，然后关闭控制面板对话框。

    %-------------------------------------------------------------------------------------------------------------------------------------------%
  • 相关阅读:
    HTML5服务器发送事件(Server-Send Events)
    无人问津的排序（一）----希尔排序
    NB二人组（二）----归并排序
    40、常用字符串格式化有哪几种？
    39、请用代码简答实现stack
    38、一行代码实现删除列表中重复的值 ?
    37、如何在函数中设置一个全局变量 ?
    NB二人组（一）----堆排序
    快排
    LOW逼三人组（三）----插入排序
  • 原文地址：https://www.cnblogs.com/luckytimor/p/14765134.html
Copyright © 2020-2023  润新知