南昌大学数字信号处理实验报告7讲解

南昌大学数字信号处理实验报告7讲解
南昌大学数字信号处理实验报告7讲解

实验七数字滤波器设计

一:实验目的

1.掌握数字巴特沃斯滤波器的设计原理和步骤

2.进一步学习用窗函数法设计FIR数字滤波器的原理及其设计步骤

7.1 无限冲激响应滤波器的阶数的估计

滤波器设计过程中的第一步是,选择接近所使用的滤波器的类型,然后由滤波器指标来估计传输函数的阶数。用来估计巴特沃兹滤波器的阶数的MATLAB命令是

[N,Wn] = buttord(Wp,WS,Rp,RS)

其中输入参数是归一化通带边界频率Wp、归一化阻带边界频率Ws、单位为dB 的通带波纹Rp和单位为dB的最小阻带衰减Rs,,由于抽样频率被假定为2Hz,Wp 和Ws均必须是0和1之间的一个数。输出数据是满足指标的最低阶数N和归一化截止频率Wn。若Rp =3dB,则Wn =Wp。buttord也可用于估计高通、带通和带阻巴特沃兹滤波器的阶数。对于高通滤波器设计,Wp>Ws。对于带通和带阻滤波器设计,Wp和Ws是指定边界频率的双元素向量,其中较低的边界频率是向量的第一个元素。在后面的情况中,Wn也是一个双元素向量。

习题:

1.用MATTAB确定一个数字无限冲激响应低通滤波器所有四种类型的最低阶

数。指标如下:40 kHz的抽样率,,4 kHz的通带边界频率,8 kHz的阻带边界频率,0.5 dB的通带波纹,40 dB的最小阻带衰减。评论你的结果。

答:标准通带边缘角频率Wp是:

标准阻带边缘角频率Ws是:

理想通带波纹Rp是0.5dB

理想阻带波纹Rs是40dB

(1)使用这些值得到巴特沃斯低通滤波器最低阶数N=8,相应的标准通带边缘频率Wn是0.2469.

(2)使用这些值得到切比雪夫1型低通滤波器最低阶数N=5,相应的标准通带边缘频率Wn是0.2000.

(3)使用这些值得到切比雪夫2型低通滤波器最低阶数N=5,相应的标准通带边缘频率Wn是0.4000.

(4)使用这些值得到椭圆低通滤波器最低阶数N=8,相应的标准通带边缘频率Wn是0.2000.

从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。

2.用MATLAB确定一个数字无限冲激响应高通滤波器所有四种类型的最低阶

数。指标如下:3500Hz的抽样率,1050 Hz的通带边界频率,600 Hz的阻带边界频率,1 dB的通带波纹,50 dB的最小阻带衰减。评论你的结果。

答:标准通带边缘角频率Wp是:

标准阻带边缘角频率Ws是:

理想通带波纹Rp是1dB

理想阻带波纹Rs是50dB

(1)使用这些值得到巴特沃斯高通滤波器最低阶数N=8,相应的标准通带边缘频率Wn是0.5646.

(2)使用这些值得到切比雪夫1型高通滤波器最低阶数N=5,相应的标准通带边缘频率Wn是0.6000.

(3)使用这些值得到切比雪夫2型高通滤波器最低阶数N=5,相应的标准通带边缘频率Wn是0.3429.

(4)使用这些值得到椭圆低通滤波器最低阶数N=4,相应的标准通带边缘频率Wn是0.6000.

从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。

3.用MATLAB确定一个数字无限冲激响应带通滤波器所有四种类型的最低阶

数。指标如下:7 kHz的抽样率,1.4 kHz和2.1 kHz的通带边界频率,1.05 kHz和

2.45 kHz的阻带边界频率,,0 .4 dB的通带波纹,50 dB的最小阻带衰减。评论

你的结果。

答:标准通带边缘角频率Wp是:

标准阻带边缘角频率Ws是:

理想通带波纹Rp是0.4dB

理想阻带波纹Rs是50dB

(1)使用这些值得到巴特沃斯带通滤波器最低阶数2N=18,相应的标准通带边缘频率Wn是[0.3835 0.6165].

(2)使用这些值得到切比雪夫1型带通滤波器最低阶数2N=12,相应的标准通带边缘频率Wn是[0.4000 0.6000].

(3)使用这些值得到切比雪夫2型带通滤波器最低阶数2N=12,相应的标准通带边缘频率Wn是[0.3000 0.7000].

(4)使用这些值得到椭圆带通滤波器最低阶数2N=8,相应的标准通带边缘频率Wn是[0.4000 0.6000].

从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。

4.用MATLAB确定一个数字无限冲激响应带阻滤波器所有四种类型的最低阶

数。指标如下:12 kHz的抽样率,2.1 kHz和4.5 kHz的通带边界频率,2.7 kHz和

3.9 kHz的阻带边界频率,0.6 dB的通带波纹,45 dB的最小阻带衰减。评论你

的结果。

答:标准通带边缘角频率Wp是:

标准阻带边缘角频率Ws是:

理想通带波纹Rp是0.6dB

理想阻带波纹Rs是45dB

(1)使用这些值得到巴特沃斯带阻滤波器最低阶数2N=18,相应的标准通带边缘频率Wn是[0.3873 0.7123].

(2)使用这些值得到切比雪夫1型带阻滤波器最低阶数2N=10,相应的标准通带边缘频率Wn是[0.3500 0.7500].

(3)使用这些值得到切比雪夫2型带阻滤波器最低阶数2N=10,相应的标准通带边缘频率Wn是[0.4500 0.6500].

(4)使用这些值得到椭圆带阻滤波器最低阶数2N=8,相应的标准通带边缘频率Wn是[0.3500 0.7500].

从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。

7.2无限冲激响应滤波器设计

程序P7.1说明巴特沃斯带阻滤波器的设计。

% 程序 P7_1

% 巴特沃斯带阻滤波器的设计

Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50;

% 估计滤波器的阶数

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs);

% 设计滤波器

[num,den] = butter(N1,Wn1,'stop');

% 显示传输函数

disp('分子系数是 ');disp(num);

disp('分母系数是 ');disp(den);

% 计算增益响应

[g, w] = gain(num,den);

% 绘制增益响应

plot(w/pi,g);grid

axis([0 1 -60 5]);

xlabel('\omega /\pi'); ylabel('增益, dB');

title('巴特沃斯带阻滤波器的增益响应');

习题:

5.通过运行程序P7. 1来设计巴特沃兹带阻滤波器。写出所产生的传输函数的准

确表达式。滤波器的指标是什么,你的设计符合指标吗,使用MATLAB,计算并绘制滤波器的未畸变的相位响应及群延迟响应。

答:表达式是:

滤波器参数是:

Wp1=0.2π,Ws1=0.4π,Ws2=0.6π,Wp2=0.8π,Rp=0.4dB,Rs=50dB. 设计的滤波器增益响应如下:

从图中可以总结出设计符合指标。

滤波器的未畸变的相位响应及群延迟响应如下:

6.修改程序P

7.1来设计符合习题Q7.1所给指标的切比雪夫1型低通滤波器。

写出所产生的传输函数的准确表达式。你的设计符合指标吗?使用MATLAB,计算并绘制滤波器的未畸变的相位响应及群延迟响应。

答:表达式如下:

设计的滤波器增益响应如下:

从图中可以总结出设计符合指标。

滤波器的未畸变的相位响应及群延迟响应如下:

7.修改程序P7.1来设计符合习题Q7.2所给指标的切比雪夫2型高通滤波器。

写出所产生的传输函数的准确表达式。你的设计符合指标吗?使用MATLAB,计算并绘制滤波器的未畸变的相位响应及群延迟响应。

答:表达式如下:

设计的滤波器增益响应如下:

从图中可以总结出设计符合指标。

滤波器的未畸变的相位响应及群延迟响应如下:

8.修改程序P7.1来设计符合习题Q7.3所给指标的椭圆带通滤波器。写出所产

生的传输函数的准确表达式。你的设计符合指标吗,使用MATLAB,计算井绘制滤波器的未畸变的相位响应及群延迟响应。

答:表达式如下:

设计的滤波器增益响应如下:

从图中可以总结出设计符合指标。

滤波器的未畸变的相位响应及群延迟响应如下:

7.3吉布斯现象

通过截短由式(7,16)、式(7.18)、式((7.20)、式(7.22)、式((7.24)和式((7.26)给出的理想滤波器的冲激响应,来设计得到有限冲激响应滤波器,然后计算它们的频率响应,可以说明吉布斯现象的发生。低通滤波器的截短的冲激响应系数可在MATLAB中使用函数的sinc二产生。该函数通过简单的修改,也可用于产生一种高通、带通或带阻滤波器的截短的冲激响应系数

习题:

9.使用函数sinc编写一个MATLAB程序,以产生截止频率在Wc= 0.4π处、长度

分别为81,61,41和21的四个零相位低通滤波器的冲激响应系数,然后计算并画出它们的幅度响应。使用冒号“:”运算符从长度为81的滤波器的冲激响

应系数中抽出较短长度滤波器的冲激响应系数。在每一个滤波器的截止频率两边研究频率响应的摆动行为。波纹的数量与滤波器的长度之间有什么关系?最大波纹的高度与滤波器的长度之间有什么关系?你将怎样修改上述程序以产生一个偶数长度的零相位低通滤波器的冲激响应系数?

答:长度为81时幅度响应如下:

长度分别为61,41和21的幅度响应如下:

从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。

波纹的数量与滤波器的长度之间的关系——波纹的数量减少与长度成正比。

最大波纹的高度与滤波器的长度之间的关系——最大波纹的高度与长度无关。

10.使用函数sinc编写一个MATLAB程序,以产生一个截止频率在Wc= 0.4π处、

长度为45的零相位高通滤波器的冲激响应系数,计算并画出其幅度响应。在每一个滤波器的截止频率两边研究频率响应的摆动行为。你将怎样修改上述程序以产生一个偶数长度的零相位高通滤波器的冲激响应系数?

答:长度为45时幅度响应如下:

从中可以观察到由于吉布斯现象产生的幅度响应摆动行为。

在这种情况下你不能改变长度。原因:这是一个零相位滤波器,这意味着它也是一个线性相位滤波器,因为零相是一种特殊的线性相位的子集。现在,理想的有限脉冲响应长度甚至有对称的中点h[n]。使其成了一个线性相位FIR 滤波器。二型滤波器不可能是高通滤波器,因为必须在z=-1处有零点,意味着w=+-π。

11.编写一个MATLAB程序,以产生长度分别为81,61,41和21的四个零相位微

分器的冲激响应系数,计算并画出它们的幅度响应。下面的代码段显示了怎样产生一个长度为2M+1的微分器。

n=1:M;

b=cos(pi*n)./n;

num=[-fliplr(b) 0 b];

对于每种情况,研究微分器的频率响应的摆动行为。波纹的数量与微分器的长度之间有什么关系,最大波纹的高度与滤波器的长度之间有什么关系?

答:幅度响应分别如下:

从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。

波纹的数量与微分器的长度之间的关系——两者成正比。

最大波纹的高度与滤波器的长度之间的关系——两者间没有关系。

12.编写一个MA11AB程序,以产生长度分别为81,61.41和21的四个离散时间

希尔伯特变换器的冲激响应系数,计算并画出它们的幅度响应。下面的代码段显示了怎样产生一个长度为2M十1的希尔伯特变换器。

n=1:M;

c=sin(pi*n)./2;

b=2*(c.*c)./(pi*n);

num=[-fliplr(b) 0 b];

对于每种情况,研究希尔伯特变换器的频率响应的摆动行为。波纹的数量与希尔伯特变换器的长度之间有什么关系?最大波纹的高度与滤波器的长度之间有什么关系?

答:幅度响应如下:

从中可以观察到由于吉布斯现像产生的幅度响应的摆动行为。

波纹的数量与希尔伯特变换器的长度之间的关系——两者成正比。

最大波纹的高度与滤波器的长度之间的关系——两者无关系。

项目7.4 FIR滤波器阶数估算

回答:

Q7.13 线性相位低通FIR滤波器的阶数估算,参数如下: ωp = 2 kHz, ωs = 2.5 kHz, δp = 0.005, δs = 0.005, 和F T = 10 kHz 使用kaiord的结果为-N = 46

使用ceil命令的目的是–朝正方向最接近整数方向取整

使用nargin命令的目的是–表明函数M文件体内变量的数目。

Q7.14(a)线性相位FIR滤波器的阶数估算,其中采样频率改为F T = 20 kHz ,则结果为-N=91.

(b) 线性相位FIR滤波器阶数的估计,其中通带波纹改成δp = 0.002 and δs = 0.002 结果为

-N=57

(c)线性相位FIR滤波器的阶数估算,其中阻带宽度改成ωs = 2.3 kHz ,结果为-N=76.

从上述结果和7.13的对比我们可以观察到:

滤波器阶数和采样频率的关系为–对于一个给定的模拟过渡带宽,采样频率的增加导致估

算阶数也相应增加,朝下一个整数取整。其中模拟过渡带宽|Fp-Fs|和Δω的关系:

Δω=2pi*|Fp-Fs|/FT。因此增加FT会减小Δω。

滤波器阶数和通带波纹宽度的关系为–估计的阶数大致和log(底数为10)成比例的扩散。

滤波器阶数和过渡带宽度的关系为–在舍入的时候,阶数随着过渡带宽成比例的改变,有两

个因素增加过渡带宽来分割顺序。

Q7.15线性相位FIR低通滤波器阶数的估算,其中滤波器满足7.13给的规格,使用kaiserord的结果为- N=54

正确结果:kaiserord([2000 2500],[1 0],[0.005 0.005],10000)

将上诉结果和7.13比较我们观察到–用凯瑟来估算阶数是较小的。因为凯瑟使用了一个不

同的近似估计。这种估计经常和FIR设计的凯瑟窗一起用。

Q7.16线性相位FIR低通滤波器的阶数估算满足的规格和7.13中的一样,使用remezord函数的结果为–N=47.

正确结果:firpmord([2000 2500],[1 0],[0.005 0.005],10000)

通过和7.13和7.15比较我们可以观察到–在这里,firpmord给了一个比凯尔更大比凯瑟更

小一点的结果。使用凯尔则更接近与一般情况。而使用凯瑟和firpmord则有专门的用途。

Q7.17 线性相位带通FIR滤波器的阶数估算满足如下规格:通带边界为1.8和3.6khz,阻带边界为

1.2khz到4.2khz,通带波纹δp = 0.01,阻带波纹δs = 0.02, 并且F T = 12 kHz, 使用kaiord

函数求得的结果为–它有一个矛盾,通带波纹δp= 0.1,得到的结果为:kaiord([1800

3600],[1200 4200],0.1,0.02,12000),然后N=20。但是当δp= 0.01时结果为:kaiord([1800

3600],[1200 4200],0.01,0.02,12000),得到的N=33。所以答案不唯一,可以选择其中一个。Q7.18线性相位带通FIR滤波器的阶数估算,其中FIR滤波器的规格和7.17一样,则使用kaiserord的结果为–同样,它也有矛盾。当使用δp= 0.1时,得到的结果为:

kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.1 0.02],12000),则N=37.当用δp= 0.01时,结

果为:kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.01 0.02],12000),此时N=45.我们可以

从两个结果中任选。

和7.17的结果比较我们观察到–通过kaiserord函数估计的阶数要更高,但如果你要设计

Kaiser窗的话则结果更精确。

Q7.19线性相位带通FIR滤波器的阶数估算,其中FIR滤波器的规格和7.17一样,使用函数remezord的结果为–当取δp= 0.01时,结果为firpmord([1200 1800 3600 4200],[0 1 0],[0.02 0.1 0.02],12000),此时N=22.而如果δp= 0.01,则结果为:firpmord([1200 1800 3600 4200],[0 1 0],[0.02 0.01 0.02],12000),此时N=35.可以从中任意选择。

和7.17和7.18的结果比较我们可以观察到–通过firpmord来估算的阶数在另外两个的中

间,在设计Parks-McClellan时更准确。

项目7.5 FIR滤波器的设计

回答:

Q7.20使用matlab程序设计并画出线性相位FIR滤波器增益和相位反应,使用fir1如下。通过使用函数kaiserord.来估计滤波器阶数,输出结果为滤波器的系数

% Program Q7_20

% Design a linear phase Lowpass FIR Digital Filter

% meeting the design specification given in Q7.13.

% - Print out the numerator coefficients

% for the transfer function.

% - Compute and plot the gain function.

% - Compute and plot the phase response.

% - Compute and plot the unwrapped phase response. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;

% Design spec as given in Q7.13.

Fp = 2*10^3;

Fs = 2.5*10^3;

FT = 10*10^3;

Rp = 0.005;

Rs = 0.005;

% Estimate the filter order and print to console

N = kaiord(Fp,Fs,Rp,Rs,FT)

% Design the filter; Hamming window by default

Wp = 2*Fp/FT; % These freqs are normalized: they go Ws = 2*Fs/FT; % zero to one, not zero to pi.

Wn = Wp + (Ws - Wp)/2;

h = fir1(N,Wn);

% Show the Numerator Coefficients

disp('Numerator Coefficients are ');disp(h);

% Compute and plot the gain response

[g, w] = gain(h,[1]); % same “gain” fcn as in Lab 4

figure(1);

plot(w/pi,g);grid;

%axis([0 1 -60 5]);

xlabel('\omega /\pi'); ylabel('Gain in dB');

title('Gain Response');

% Compute the frequency response

w2 = 0:pi/511:pi;

Hz = freqz(h,[1],w2);

% TEST: did we meet the spec?

MagH = abs(Hz);

T1 = 1.005*ones(1,length(w2));

T2 = 0.995*ones(1,length(w2));

相关主题
相关文档
最新文档