IIR 数字滤波器设计与应用的 MATLAB 实现

准备前的工作——一段收集来的代码

一、输入信号时域&频域图像（非声音信号，而是函数）

f1=10;%第一个点频信号分量频率
f2=30;%第二个点频信号分量频率
f3=45;%第三个点频信号分量频率
fs=100;%采样率
T=2;%间隔长度
n=round(T*fs);%采样点个数
t=linspace(0,T,n);
y=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t)+randn(size(t));
figure;
subplot(2,1,1);
plot(t,y);
title('输入信号时域图像');
xlabel('t/s');
ylabel('V');
fft_y=fftshift(fft(y));
f=linspace(-fs/2,fs/2,n);
subplot(2,1,2);
plot(f,abs(fft_y));
title('输入信号频域图像');
xlabel('f/Hz');
ylabel('V');

https://blog.****.net/Qinlong_Stm32/article/details/84570524
文章介绍的简单易懂↑↑↑
——————————————————————————————

思考流程：

接下来，我需要对该信号进行处理：

首先，如何加指定频段加入固定幅值的噪声？频域图像的横坐标又该如何确定？
其次，如何设计IIR滤除指定频段的噪声？该频谱图又该如何确定横坐标？
最后，滤波前和滤波后的误差曲线如何画？误差均方值又该如何计算？

一、输入声音信号的时域图像和频域图像

clear;
clc;
%输入原始信号
[y,fs]=audioread('E:\son\dudu.wav'); %将声音放于matlab中
info=audioinfo('E:\son\dudu.wav')
%sound(y,fs);
T=1/fs; %采样时间
t=(0:length(y)-1)*T;%时间
f=(0:length(y)-1)*fs/length(y);
figure(1);
yz=y(:,1);%左声道
subplot(2,1,1);
plot(t,yz);%输入信号时域曲线
title('原始信号时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
n=length(yz);%进行变换的点数
y1=fft(yz,n); %对n点进行傅里叶变换到频域
F=fs/length(yz); %谱分辨率，频谱间隔
plot(f,abs(y1));%左声道频谱图
axis([0,1000,0,10000]);
title('原始信号频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
grid on

Note: 在plot频谱图的时候，我取的范围是（-fs/2,fs/2-F）,为什么不取（0，fs/2-F）呢？
因为我接下来要在100,500,900处加噪声，前者的取法会使中间频率集中在0附近，这样加噪声会对声音有影响；而后者会使频率集中在2×10^4处，前面的部分几乎振幅为0，加噪声后也没有影响。

如何选取坐标范围是根据需求或者要求来的。但需要注意的是，如果y需要选取0~1000hz的，那么对应的f也应该改变。

二、加100、500、900Hz的噪声

加上噪声以后会出现，这样的错误显示。原因是因为相加的时候维度搞错了。yz是1×n的列向量，x是n×1的行向量，无法相加。我一开始并没有想到这上面，这是一种基础的语法错误，我需要的是经验。

#加噪声
noise1=sin(2*pi*100*t);
noise2=sin(2*pi*500*t);
noise3=sin(2*pi*900*t);
x1=noise1+noise2+noise3;
x1=yz+x1';
figure(2);
subplot(2,1,1);
plot(t,x1);
title('加噪信号时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
nx=length(x1);
x=fft(x1,nx);
plot(f,abs(x));%左声道频谱图
axis([0,1000,0,10000]);
title('加噪信号频谱');
xlabel('F(Hz)');
ylabel('H(jw)');

三、设计IIR滤波器

数字滤波器的设计方法，三个步骤：
1、给出所需要的滤波器的技术指标；
2、设计一个使其逼近所需要的技术指标；
实现所设计的。
3、IIR滤波器的设计主要借助模拟滤波器z做转换来进行；FIR的设计主要建立在对理想滤波器频率特性作某种近似的基础上，有窗函数法、频率抽样法等；
线性相位滤波器通常采用FIR型，其单位脉冲响满足一定条件时，其相位特性是严格线性的。

一、IIR滤波器的设计步骤如下：

1)、按一定规则把给定的数字滤波器技术指标转换为模拟低通滤波器的技术指标；
2)、根据转换后的技术指标设计模拟低通滤波器系统函数H(s)；
3)、再按一定规则将转换成；若所设计的数字滤波器是低通的，那么完成整个设计工作；
4)、将高通，带通或带阻数字滤波器的技术指标转化为低通模拟滤波器的技术指标，然后从第(2)步开始设计低通滤波器H(s)，再将转换为所需的H(z)。

我现在需要确定基本参数，根据该参数设计滤波器。如何确定参数？——fdatool；确定参数以后如何写代码？

1.参数可以通过fdatool确定，其实我觉得直接自定义一个即可。衰减一个是0.1，一个是60dB。滤波器貌似都是这个标准(不知道为什么)
2.如何将滤波器图形转化成代码语言——可以参考书本，很简单的定义
3.我现在的问题是，如何将滤波器加到声音信号里实现滤波。

---

######我有点累了。。。。。现在搞出来一个半成品（因为不太确定这个图像和声音对不对总感觉很怪+对于代码有一些知识盲点未完全吃透），这些知识点需要攻破，比如每个滤波器的公式这么多选哪个才好？图像的输出的横坐标怎么确定？
######算了，先写FIR的代码吧。IIR晚点再继续写。

---

四、待滤波信号加滤波器

%%%%%%%###带阻滤波器100Hz####%%%%%
x2=filter(ww100,x1);
figure(3);
subplot(2,1,1);
plot(t,x2);
title('100HZ被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x2)));
axis([0 1000 0 10000]);
title('100Hz被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
%%%%%%%#####带阻滤波器500Hz#####%%%%%%%%%
x3=filter(ww500,x2);
figure(4);
subplot(2,1,1);
plot(t,x3);
title('100HZ、500HZ被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x3)));
axis([0 1000 0 10000]);
title('100Hz、500Hz被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
%%%%%%%%%########带阻滤波器900Hz######%%%%%%%%%
x4=filter(ww900,x3);
figure(5);
subplot(2,1,1);
plot(t,x4);
title('全被滤掉后的时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x4)));
axis([0 1000 0 10000]);
title('全部被滤后的频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
sound(x4,fs);

五、滤波前和滤波后的误差曲线

yend=yz-x4;  %误差
figure(6);
subplot(2,1,1);
plot(t,yend);  %误差时域曲线
title('误差时域');
xlabel('时间');
ylabel('振幅');
subplot(2,1,2);
yp=abs(y1)-abs(fft(x4));%此处要取绝对值，不然会将相位也算进去
plot(f,abs(yp));  %误差频谱分布
axis([0 1000 0 10000]);
title('误差频谱');
xlabel('F(Hz)');
ylabel('H(jw)');
S=std(yend)  %均方值

六、不同滤波器的定义

100Hz带阻滤波器

function Hd = ww100
Fs = 48000;  % Sampling Frequency
Fpass1 = 95;          % First Passband Frequency
Fstop1 = 96;          % First Stopband Frequency
Fstop2 = 102;         % Second Stopband Frequency
Fpass2 = 103;         % Second Passband Frequency
Apass1 = 0.1;         % First Passband Ripple (dB)
Astop  = 70;          % Stopband Attenuation (dB)
Apass2 = 0.1;         % Second Passband Ripple (dB)
match  = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its ELLIP method.
h  = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
                      Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);

500Hz带阻滤波器

function Hd = ww500
Fs = 48000;  % Sampling Frequency
Fpass1 = 497;         % First Passband Frequency
Fstop1 = 498;         % First Stopband Frequency
Fstop2 = 503;         % Second Stopband Frequency
Fpass2 = 504;         % Second Passband Frequency
Apass1 = 0.1;         % First Passband Ripple (dB)
Astop  = 70;          % Stopband Attenuation (dB)
Apass2 = 0.1;         % Second Passband Ripple (dB)
match  = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its ELLIP method.
h  = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
                      Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);

900Hz带阻滤波器

function Hd = ww900
Fs = 48000;  % Sampling Frequency
Fpass1 = 895;         % First Passband Frequency
Fstop1 = 897;         % First Stopband Frequency
Fstop2 = 903;         % Second Stopband Frequency
Fpass2 = 905;         % Second Passband Frequency
Apass1 = 0.1;         % First Passband Ripple (dB)
Astop  = 60;          % Stopband Attenuation (dB)
Apass2 = 0.1;         % Second Passband Ripple (dB)
match  = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its ELLIP method.
h  = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
                      Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);

心得

1、如果题目中给出的是3*pi rad，可以直接使用公式，也可以换算成Hz.如果换算成HZ,公式如下：

1rad/sec=1rad Hz
1Hz=2pi rad/sec
50Hz=100pi rad/sec

2、其实实现带通和带阻，不过是将通带频率和阻带频率换成了坐标范围

思考题

1）本实验中Butterworth，Chebyshev，Elliptic三种滤波器有什么差别？对滤波效果有什么影响？
相同阶数时：
巴特沃斯滤波器通带最平坦，阻带下降慢。
切比雪夫滤波器通带等纹波，阻带下降较快。
椭圆滤波器，椭圆滤波器在通带等纹波（阻带平坦或等纹波），阻带下降最快。
2）滤波器的阻带宽度可以设到多小？通带最小衰减可以设到多小？阻带最大衰减可以设到多大？
1、通带纹波是指在滤波器的频响中通带的最大幅值和最小幅值之间的差值，正常的纹波一般小于1db。不过也视情况而言，通带纹波会导致通带内的幅值大小有变化，一般要求越高，纹波越小越好。通带纹波和滤波器的阶数有关系，阶数越大纹波越小。
2、 阻带衰减：在通带中，有部分信号通，部分信号阻，而阻的部分不能不能全部阻断，只有部分衰减，部分留了下来，最小衰减描述了阻碍受阻信号的能力，衰减越大，则能力越好。
3）在本实验中，双线性变换与脉冲不变响应对滤波效果有什么影响？

双线性变换实现低通滤波器

%双线性变换前预畸
Fs=500;
wp=(100/Fs)*2*pi;
ws=(200/Fs)*2*pi;
Rp=2;
Rs=15;
wp2=2*Fs*tan(wp/2);
ws2=2*Fs*tan(ws/2);
%选择滤波器的最小阶数
[N,wn]=buttord(wp2,ws2,Rp,Rs,'s');%注意此处输入的是畸变后的指标，输出N为符合要求的模拟滤波器的最小阶数，wn为3dB带宽
%创建butterworth模拟滤波器
[Z,P,K]=buttap(N);
%把滤波器零极点模型转化为传递函数模型
[Bap,Aap]=zp2tf(Z,P,K);
%把模拟滤波器原型转换为截止频率为wn的模拟低通滤波器
[b,a]=lp2lp(Bap,Aap,wn);
%用双线性法实现模拟滤波器到数字滤波器的转换
[bz,az]=bilinear(b,a,Fs);
%绘制频率响应曲线
[H,W]=freqz(bz,az);
subplot(2,1,1);
plot(W/pi,abs(H));
grid;
xlabel('频率w/pi');
ylabel('幅度绝对值');
subplot(2,1,2);
plot(W/pi,20*log10((abs(H)+eps)/max(abs(H))));
grid;
xlabel('频率w/pi');
ylabel('幅度dB');
————————————————
版权声明：本文为****博主「Zhang_P_Y」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.****.net/lg1259156776/article/details/48603575/

双线性变换实现带通滤波器

https://wenku.baidu.com/view/3f4689221b37f111f18583d049649b6648d70994.html

扩

一、IIR设计双线性变换高通滤波器

https://blog.****.net/weixin_30415801/article/details/95458429

f_p     =   15000     ;   %阻带上截止频率
f_st    =   30000     ;   %通带下截止频率
R_p     =   0.4       ;   %通带允许的最大衰减
R_st    =   30      ;   %阻带允许的最小衰减
f_s     =   100000    ;   %采样频率
T_s     =   1 / f_s ;   %采样间隔
%1.将数字高通滤波器的频率参数变换为归一化的数字角频率参数
omega_p  = 2 * pi * f_p  / f_s;
omega_st = 2 * pi * f_st / f_s;
%2.预畸变处理,将归一化数字角频率参数变换成模拟高通滤波器的角频率参数
C = 2 / T_s ;
Omega_p   = C * tan( omega_p  / 2 );
Omega_st  = C * tan( omega_st / 2 ); 
%3.对模拟高通滤波器的角频率参数做归一化处理
lamda_p  = Omega_p  / Omega_p;
lamda_st = Omega_st / Omega_p;
%4.设计归一化模拟滤波器,得到归一化模拟高通滤波器的转移函数
[ N , Wn ] = buttord( Omega_p , Omega_st , R_p , R_st , 's' ); %选择模拟巴特沃斯滤波器的最小阶数
[ z , p , k ] = buttap(N); %创建巴特沃斯模拟低通滤波器
[ Bp , Ap ] = zp2tf(z,p,k); %由零点、极点、增益确定传输函数的分子与分母系数
%5.利用归一化模拟高通滤波器的转移函数确定模拟高通滤波器的转移函数
[ b , a ] = lp2hp(Bp,Ap,Wn);
%6.利用模拟高通滤波器的转移函数确定IIR数字滤波器的转移函数 
[ bz , az ] = bilinear(b,a,f_s);

figure(1);
freqz(bz,az);
title('高通滤波器幅度谱和相位谱特性');

设计步骤

带通滤波器

https://wenku.baidu.com/view/3f4689221b37f111f18583d049649b6648d70994.html
https://blog.****.net/weixin_30415801/article/details/95458429

f_s1     =   150     ;   %阻带上截止频率
f_s2    =  600;
f_p1    =   200     ;   %通带下截止频率
f_p2    =   500;
R_p     =   0.5       ;   %通带允许的最大衰减
R_s    =   40     ;   %阻带允许的最小衰减
f_s     =   2000    ;   %采样频率
T_s     =   1 / f_s ;   %采样间隔
%1.将数字高通滤波器的频率参数变换为归一化的数字角频率参数
omega_p1  = 2 * pi * f_p1  /f_s;
omega_s1 = 2 * pi * f_s1 /f_s;
omega_p2  = 2 * pi * f_p2  /f_s;
omega_s2 = 2 * pi * f_s2 /f_s;
%2.预畸变处理,将归一化数字角频率参数变换成模拟高通滤波器的角频率参数
C = 2 / T_s ;
Omega_p1   =C*tan( omega_p1  / 2 );
Omega_s1  = C*tan( omega_s1 / 2 ); 
Omega_p2   = C*tan( omega_p2  / 2 );
Omega_s2  = C*tan( omega_s2 / 2 );
W0=sqrt(Omega_p1*Omega_p2);
Omega_BW=Omega_p1-Omega_p2;
%3.对模拟高通滤波器的角频率参数做归一化处理
eta_sl      =   Omega_sl / Omega_BW;
eta_p1       =   Omega_p1  / Omega_BW;   
eta_p2       =   Omega_p2  / Omega_BW;
eta_s2      =   Omega_s2 / Omega_BW;
%4.设计归一化模拟滤波器,得到归一化模拟高通滤波器的转移函数
wph=[Omega_p1,Omega_p2 ];
wsh=[Omega_s1,Omega_s2];
[ N , Wn ] = buttord( wph , wsh , R_p , R_st , 's' ); %选择模拟巴特沃斯滤波器的最小阶数
[ z , p , k ] = buttap(N); %创建巴特沃斯模拟低通滤波器
%[ Bp , Ap ] = zp2tf(z,p,k); %由零点、极点、增益确定传输函数的分子与分母系数
%5.利用归一化模拟高通滤波器的转移函数确定模拟高通滤波器的转移函数
BB=real(poly(z));
AA=real(poly(p));
[ b , a ] = lp2bp(BB,AA,W0,Omega_BW);
%6.利用模拟高通滤波器的转移函数确定IIR数字滤波器的转移函数 
[ bz , az ] = bilinear(b,a,f_s);

figure(1);
freqz(bz,az);
title('带通滤波器幅度谱和相位谱特性');