数字信号处理实验
实验报告
(第五章:连续时间信号的数字处理)
学 院: 信息工程学院 专 业: 通信工程 班 级:
学 号:
学生姓名:
Q7.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.
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。
Q7.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.
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 Q7.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].
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 Q7.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].
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 程序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(' 巴特沃斯带阻滤波器的增益响应' );
Q7.5通过运行程序P7. 1来设计巴特沃兹带阻滤波器。写出所产生
的传输函数的准确表达式。滤波器的指标是什么,你的设计符合指标吗,使用MA TLAB ,计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式是:
滤波器参数是:
Wp1=0.2π,Ws1=0.4π,Ws2=0.6π,Wp2=0.8π,Rp=0.4dB,Rs=50dB.
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.6修改程序P7.1来设计符合习题Q7.1所给指标的切比雪夫1型
低通滤波器。写出所产生的传输函数的准确表达式。你的设计符合指标吗? 使用MATLAB, 计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.7修改程序P7.1来设计符合习题Q7.2所给指标的切比雪夫2型
高通滤波器。写出所产生的传输函数的准确表达式。你的设计符合指标吗? 使用MATLAB, 计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.8修改程序P7.1来设计符合习题Q7.3所给指标的椭圆带通滤波
器。写出所产生的传输函数的准确表达式。你的设计符合指标吗,使用MATLAB ,计算井绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
通过截短由式(7,16)、式(7.18)、式((7.20)、式(7.22)、式((7.24)和式((7.26)给出的理想滤波器的冲激响应,来设计得到有限冲激响应滤波器,然后计算它们的频率响应,可以说明吉布斯现象的发生。低通滤波器的截短的冲激响应系数可在MA TLAB 中使用函数的sinc 二产生。该函数通过简单的修改,也可用于产生一种高通、带通或带阻滤波器的截短的冲激响应系数
习题:
Q7.9使用函数sinc 编写一个MATLAB 程序,以产生截止频率在
Wc= 0.4π处、长度分别为81,61,41和21的四个零相位低通滤波器的冲激响应系数,然后计算并画出它们的幅度响应。使用冒号“:”运算符从长度为81的滤波器的冲激响应系数中抽出较短长度滤波器的冲激响应系数。在每一个滤波器的截止频率两边研究频率响应的摆动行为。波纹的数量与滤波器的长度之间有什么关系? 最大波纹的高度与滤波器的长度之间有什么关系? 你将怎样修改上述程序以产生一个偶数长度的零相位低通滤波器的冲激响应系数?
长度为81时幅度响应如下:
长度分别为61,41和21的幅度响应如下:
从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。 波纹的数量与滤波器的长度之间的关系——波纹的数量减少与长度成正比。
最大波纹的高度与滤波器的长度之间的关系——最大波纹的高度与长度无关。
Q7.10使用函数sinc 编写一个MATLAB 程序,以产生一个截止频率
在Wc= 0.4π处、长度为45的零相位高通滤波器的冲激响应
系数,计算并画出其幅度响应。在每一个滤波器的截止频率
两边研究频率响应的摆动行为。你将怎样修改上述程序以产
生一个偶数长度的零相位高通滤波器的冲激响应系数?
长度为45时幅度响应如下:
从中可以观察到由于吉布斯现象产生的幅度响应摆动行为。 在这种情况下你不能改变长度。原因:这是一个零相位滤波器,这意味着它也是一个线性相位滤波器,因为零相是一种特殊的线性相位的子集。现在,理想的有限脉冲响应长度甚至有对称的中点h[n]。使其成了一个线性相位FIR 滤波器。二型滤波器不可能是高通滤波器,因为必须在z=-1处有零点,意味着w=+-π。 Q7.11编写一个MATLAB 程序,以产生长度分别为81,61,41和21
的四个零相位微分器的冲激响应系数,计算并画出它们的幅
度响应。下面的代码段显示了怎样产生一个长度为2M+1的
微分器。
n=1:M;
b=cos(pi*n)./n;
num=[-fliplr(b) 0 b];
对于每种情况, 研究微分器的频率响应的摆动行为。波纹的数量
与微分器的长度之间有什么关系,最大波纹的高度与滤波器的长度之间有什么关系?
幅度响应分别如下:
从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。 波纹的数量与微分器的长度之间的关系——两者成正比。
最大波纹的高度与滤波器的长度之间的关系——两者间没有关系。 Q7.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];
对于每种情况,研究希尔伯特变换器的频率响应的摆动行为。波纹的数量与希尔伯特变换器的长度之间有什么关系? 最大波纹的高度与滤波器的长度之间有什么关系?
幅度响应如下:
从中可以观察到由于吉布斯现像产生的幅度响应的摆动行为。 波纹的数量与希尔伯特变换器的长度之间的关系——两者成正比。
最大波纹的高度与滤波器的长度之间的关系——两者无关系。 Q7.13 线性相位低通FIR 滤波器的阶数估算,参数如下:
ωp =2 kHz, ωs =2.5 kHz, δp = 0.005, δs = 0.005, FT = 10kHz 使用 kaiord 的结果为N = 46
使用 ceil 命令的目的是朝正方向最接近整数方向取整。
使用nargin 命令的目的是表明函数M 文件体内变量的数目。 Q7.14 (a)线性相位FIR 滤波器的阶数估算,其中采样频率改为F T = 20
kHz ,则结果为 N=91。
(b) 线性相位FIR 滤波器阶数的估计,其中通带波纹改成δp = 0.002和
δ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,FT = 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.020.01 0.02],12000),此时N=35.可以从中任意选择。
和7.17和7.18的结果比较我们可以观察到通过firpmord 来估算的阶数在另外两个的中间,在设计Parks-McClellan 时更准确。 Q7.20 使用matlab 程序设计并画出线性相位FIR 滤波器增益和相位
反应,使用fir1如下。通过使用函数kaiserord. 来估计滤波器阶数,输出结果为滤波器的系数。低通滤波器满足7.20所要求的规格的系数如下:
增益和相位响应如下:
从增益图像我们可以知道这个设计不能满足规格. 这个滤波器满足规格的阶数为N=66.为了满足规格,图如下:
Q7.21 汉宁窗:
布莱克曼窗:
切比雪夫窗:
Q7.22 程序如下:
clear;
Fp = 2*10^3;
Fs = 2.5*10^3;
FT = 10*10^3;
Rp = 0.005;
Rs = 0.005;
N = kaiord(Fp,Fs,Rp,Rs,FT)
Wp = 2*Fp/FT;
Ws = 2*Fs/FT;
F = [0 Wp Ws 1];
A = [1 1 0 0];
h = firpm(N,F,A);
disp('Numerator Coefficients are ' );
disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
MagH = abs(Hz);
T1 = 1.005*ones(1,length(w2));
T2 = 0.995*ones(1,length(w2));
T3 = 0.005*ones(1,length(w2));
figure(4);
plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid;
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' );
title('Unwrapped Phase Response' );
低通滤波器系数:
0.0028 -0.0022 -0.0046 -0.0006 0.0053 0.0019 -0.0073 -0.0058 0.0079 0.0106 -0.0069 -0.0170 0.0032 0.0243 0.0045 -0.0319 -0.0182 0.0390 0.0422 -0.0448 -0.0924 0.0486 0.3136 0.4501 0.3136 0.0486 -0.0924 -0.0448 0.0422 0.0390 -0.0182 -0.0319 0.0045 0.0243 0.0032 -0.0170 -0.0069 0.0106 0.0079 -0.0058 -0.0073 0.0019 0.0053 -0.0006 -0.0046 -0.0022 0.0028
增益和相位响应:
从图中可以看出此时的滤波器不满足指标。欲满足指标,应调节N=47. Q7.23 用凯泽窗设计一个有限冲激响应低通滤波器。
程序:
clear;
Wp = 0.31;
Ws = 0.41;
Wn = Wp + (Ws-Wp)/2;
As = 50;
Ds = 10^(-As/20);
Dp = Ds;
if As > 21
N = ceil((As-7.95)*2/(14.36*(abs(Wp-Ws)))+1)
else
N = ceil(0.9222*2/abs(Wp-Ws)+1)
end
if As > 50
BTA = 0.1102*(As-8.7);
elseif As >= 21
BTA = 0.5842*(As-21)^0.4+0.07886*(As-21);
else
BTA = 0;
end
Win = kaiser(N+1,BTA);
h = fir1(N,Wn,Win);
disp('Numerator Coefficients are ' );
disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' );
title('Unwrapped Phase Response' );
低通滤波器系数:
0.0003 0.0008 0.0003 -0.0011 -0.0017 0.0000 0.0026 0.0027 -0.0010 -0.0049 -0.0035 0.0033 0.0080 0.0034 -0.0074 -0.0119 -0.0018 0.0140 0.0161 -0.0027 -0.0241 -0.0201 0.0127 0.0406 0.0236 -0.0354 -0.0754 -0.0258 0.1214 0.2871 0.3597 0.2871 0.1214 -0.0258 -0.0754 -0.0354 0.0236 0.0406 0.0127 -0.0201 -0.0241 -0.0027 0.0161 0.0140 -0.0018 -0.0119 -0.0074 0.0034 0.0080 0.0033 -0.0035 -0.0049 -0.0010 0.0027 0.0026 0.0000 -0.0017 -0.0011 0.0003 0.0008 0.0003 增益和相位响应如下:
从图中可以看出设计的滤波器满足要求。N=60. Q7.24 用函数kaiserord 和firl 重做习题Q7.23
程序:
clear;
Wp = 0.31;
Ws = 0.41;
As = 50;
Ds = 10^(-As/20);
F = [Wp Ws];
A = [1 0];
DEV = [Ds Ds];
[N,Wn,BTA,Ftype] = kaiserord(F,A,DEV);
Win = kaiser(N+1,BTA);
h = fir1(N,Wn,Ftype,Win);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
参数如下:
Wp = 0.31; Ws = 0.41; As = 50 dB
增益和相位响应如下:
从图中可以看出设计的滤波器满足要求。N=59. Q7.25 用fir2设计一个95阶有限冲激响应滤波器。
程序:
clear;
N = 95;
A = [0.4 0.4 1.0 1.0 0.8 0.8];
F = [0 0.25 0.3 0.45 0.5 1.0];
h = fir2(N,F,A);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
plot(w2/pi,abs(Hz));grid;
xlabel('\omega /\pi' ); ylabel('|H(e^{j\omega})|' ); title('|H(e^{j\omega})|' );
figure(3);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(4);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
幅度响应:
从幅度响应中可以看出,此滤波器满足指标。 Q7.26 使用remez 设计有限冲激响应带通滤波器。 程序:
clear;
F = [1200 1800 3600 4200];
A = [0 1 0];
DEV = [0.02 0.1 0.02];
Fs = 12000;
Dp = 0.1;
Ds = 0.02;
[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs); F2 = 2*[0 1200 1800 3600 4200 6000]/Fs; A2 = [0 0 1 1 0 0];
wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds]; h = firpm(N,F2,A2,wgts);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' ); title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
增益响应: 相位响应:
从增益响应的图像中可以看出,此滤波器满足指标。N=37. Q7.27 用remez 设计具有如下指标的有限冲激响应带通滤波器。
程序:
clear;
Fs1 = 1500;
Fp1 = 1800;
Fp2 = 3000;
Fs2 = 4200;
Fs = 12000;
Dp = 0.1;
Ds = 0.02;
F = [Fs1 Fp1 Fp2 Fs2];
A = [0 1 0];
DEV = [Ds Dp Ds];
[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs);
ws1 = 2*Fs1/Fs;
wp1 = 2*Fp1/Fs;
wp2 = 2*Fp2/Fs;
ws2 = 2*Fs2/Fs;
F2 = [0 ws1 wp1 wp2 ws2 1];
A2 = [0 0 1 1 0 0];
wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds];
h = firpm(N,F2,A2,wgts);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
figure(4);
hold on ;
tmpY = 0:1.4/4:1.4;
tmpX = ones(1,length(tmpY))*wp1;
plot(tmpX,tmpY, 'r-' );
tmpX = ones(1,length(tmpY))*wp2;
plot(tmpX,tmpY, 'r-' );
tmpX = ones(1,length(tmpY))*ws1;
plot(tmpX,tmpY, 'g-' );
tmpX = ones(1,length(tmpY))*ws2;
plot(tmpX,tmpY, 'g-' );
tmpY = ones(1,length(w))*(Dp);
plot(w/pi,tmpY, 'r-' );
tmpY = ones(1,length(w))*(Ds);
plot(w/pi,tmpY, 'g-' );
plot(w2/pi,abs(Hz));grid;
hold off ;
增益响应:
该滤波器不满足设计指标。用
估计滤波器的阶数N=73。 从增益响应的图像可以看出,kaiserord
数字信号处理实验
实验报告
(第五章:连续时间信号的数字处理)
学 院: 信息工程学院 专 业: 通信工程 班 级:
学 号:
学生姓名:
Q7.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.
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。
Q7.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.
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 Q7.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].
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 Q7.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].
从以上结果中观察到椭圆滤波器的阶数最低,并且符合要求。 程序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(' 巴特沃斯带阻滤波器的增益响应' );
Q7.5通过运行程序P7. 1来设计巴特沃兹带阻滤波器。写出所产生
的传输函数的准确表达式。滤波器的指标是什么,你的设计符合指标吗,使用MA TLAB ,计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式是:
滤波器参数是:
Wp1=0.2π,Ws1=0.4π,Ws2=0.6π,Wp2=0.8π,Rp=0.4dB,Rs=50dB.
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.6修改程序P7.1来设计符合习题Q7.1所给指标的切比雪夫1型
低通滤波器。写出所产生的传输函数的准确表达式。你的设计符合指标吗? 使用MATLAB, 计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.7修改程序P7.1来设计符合习题Q7.2所给指标的切比雪夫2型
高通滤波器。写出所产生的传输函数的准确表达式。你的设计符合指标吗? 使用MATLAB, 计算并绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
Q7.8修改程序P7.1来设计符合习题Q7.3所给指标的椭圆带通滤波
器。写出所产生的传输函数的准确表达式。你的设计符合指标吗,使用MATLAB ,计算井绘制滤波器的未畸变的相位响应及群延迟响应。
表达式如下:
设计的滤波器增益响应如下:
从图中可以总结出设计符合指标。
滤波器的未畸变的相位响应及群延迟响应如下:
通过截短由式(7,16)、式(7.18)、式((7.20)、式(7.22)、式((7.24)和式((7.26)给出的理想滤波器的冲激响应,来设计得到有限冲激响应滤波器,然后计算它们的频率响应,可以说明吉布斯现象的发生。低通滤波器的截短的冲激响应系数可在MA TLAB 中使用函数的sinc 二产生。该函数通过简单的修改,也可用于产生一种高通、带通或带阻滤波器的截短的冲激响应系数
习题:
Q7.9使用函数sinc 编写一个MATLAB 程序,以产生截止频率在
Wc= 0.4π处、长度分别为81,61,41和21的四个零相位低通滤波器的冲激响应系数,然后计算并画出它们的幅度响应。使用冒号“:”运算符从长度为81的滤波器的冲激响应系数中抽出较短长度滤波器的冲激响应系数。在每一个滤波器的截止频率两边研究频率响应的摆动行为。波纹的数量与滤波器的长度之间有什么关系? 最大波纹的高度与滤波器的长度之间有什么关系? 你将怎样修改上述程序以产生一个偶数长度的零相位低通滤波器的冲激响应系数?
长度为81时幅度响应如下:
长度分别为61,41和21的幅度响应如下:
从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。 波纹的数量与滤波器的长度之间的关系——波纹的数量减少与长度成正比。
最大波纹的高度与滤波器的长度之间的关系——最大波纹的高度与长度无关。
Q7.10使用函数sinc 编写一个MATLAB 程序,以产生一个截止频率
在Wc= 0.4π处、长度为45的零相位高通滤波器的冲激响应
系数,计算并画出其幅度响应。在每一个滤波器的截止频率
两边研究频率响应的摆动行为。你将怎样修改上述程序以产
生一个偶数长度的零相位高通滤波器的冲激响应系数?
长度为45时幅度响应如下:
从中可以观察到由于吉布斯现象产生的幅度响应摆动行为。 在这种情况下你不能改变长度。原因:这是一个零相位滤波器,这意味着它也是一个线性相位滤波器,因为零相是一种特殊的线性相位的子集。现在,理想的有限脉冲响应长度甚至有对称的中点h[n]。使其成了一个线性相位FIR 滤波器。二型滤波器不可能是高通滤波器,因为必须在z=-1处有零点,意味着w=+-π。 Q7.11编写一个MATLAB 程序,以产生长度分别为81,61,41和21
的四个零相位微分器的冲激响应系数,计算并画出它们的幅
度响应。下面的代码段显示了怎样产生一个长度为2M+1的
微分器。
n=1:M;
b=cos(pi*n)./n;
num=[-fliplr(b) 0 b];
对于每种情况, 研究微分器的频率响应的摆动行为。波纹的数量
与微分器的长度之间有什么关系,最大波纹的高度与滤波器的长度之间有什么关系?
幅度响应分别如下:
从中可以观察到由于吉布斯现象产生的幅度响应的摆动行为。 波纹的数量与微分器的长度之间的关系——两者成正比。
最大波纹的高度与滤波器的长度之间的关系——两者间没有关系。 Q7.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];
对于每种情况,研究希尔伯特变换器的频率响应的摆动行为。波纹的数量与希尔伯特变换器的长度之间有什么关系? 最大波纹的高度与滤波器的长度之间有什么关系?
幅度响应如下:
从中可以观察到由于吉布斯现像产生的幅度响应的摆动行为。 波纹的数量与希尔伯特变换器的长度之间的关系——两者成正比。
最大波纹的高度与滤波器的长度之间的关系——两者无关系。 Q7.13 线性相位低通FIR 滤波器的阶数估算,参数如下:
ωp =2 kHz, ωs =2.5 kHz, δp = 0.005, δs = 0.005, FT = 10kHz 使用 kaiord 的结果为N = 46
使用 ceil 命令的目的是朝正方向最接近整数方向取整。
使用nargin 命令的目的是表明函数M 文件体内变量的数目。 Q7.14 (a)线性相位FIR 滤波器的阶数估算,其中采样频率改为F T = 20
kHz ,则结果为 N=91。
(b) 线性相位FIR 滤波器阶数的估计,其中通带波纹改成δp = 0.002和
δ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,FT = 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.020.01 0.02],12000),此时N=35.可以从中任意选择。
和7.17和7.18的结果比较我们可以观察到通过firpmord 来估算的阶数在另外两个的中间,在设计Parks-McClellan 时更准确。 Q7.20 使用matlab 程序设计并画出线性相位FIR 滤波器增益和相位
反应,使用fir1如下。通过使用函数kaiserord. 来估计滤波器阶数,输出结果为滤波器的系数。低通滤波器满足7.20所要求的规格的系数如下:
增益和相位响应如下:
从增益图像我们可以知道这个设计不能满足规格. 这个滤波器满足规格的阶数为N=66.为了满足规格,图如下:
Q7.21 汉宁窗:
布莱克曼窗:
切比雪夫窗:
Q7.22 程序如下:
clear;
Fp = 2*10^3;
Fs = 2.5*10^3;
FT = 10*10^3;
Rp = 0.005;
Rs = 0.005;
N = kaiord(Fp,Fs,Rp,Rs,FT)
Wp = 2*Fp/FT;
Ws = 2*Fs/FT;
F = [0 Wp Ws 1];
A = [1 1 0 0];
h = firpm(N,F,A);
disp('Numerator Coefficients are ' );
disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
MagH = abs(Hz);
T1 = 1.005*ones(1,length(w2));
T2 = 0.995*ones(1,length(w2));
T3 = 0.005*ones(1,length(w2));
figure(4);
plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid;
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' );
title('Unwrapped Phase Response' );
低通滤波器系数:
0.0028 -0.0022 -0.0046 -0.0006 0.0053 0.0019 -0.0073 -0.0058 0.0079 0.0106 -0.0069 -0.0170 0.0032 0.0243 0.0045 -0.0319 -0.0182 0.0390 0.0422 -0.0448 -0.0924 0.0486 0.3136 0.4501 0.3136 0.0486 -0.0924 -0.0448 0.0422 0.0390 -0.0182 -0.0319 0.0045 0.0243 0.0032 -0.0170 -0.0069 0.0106 0.0079 -0.0058 -0.0073 0.0019 0.0053 -0.0006 -0.0046 -0.0022 0.0028
增益和相位响应:
从图中可以看出此时的滤波器不满足指标。欲满足指标,应调节N=47. Q7.23 用凯泽窗设计一个有限冲激响应低通滤波器。
程序:
clear;
Wp = 0.31;
Ws = 0.41;
Wn = Wp + (Ws-Wp)/2;
As = 50;
Ds = 10^(-As/20);
Dp = Ds;
if As > 21
N = ceil((As-7.95)*2/(14.36*(abs(Wp-Ws)))+1)
else
N = ceil(0.9222*2/abs(Wp-Ws)+1)
end
if As > 50
BTA = 0.1102*(As-8.7);
elseif As >= 21
BTA = 0.5842*(As-21)^0.4+0.07886*(As-21);
else
BTA = 0;
end
Win = kaiser(N+1,BTA);
h = fir1(N,Wn,Win);
disp('Numerator Coefficients are ' );
disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' );
title('Unwrapped Phase Response' );
低通滤波器系数:
0.0003 0.0008 0.0003 -0.0011 -0.0017 0.0000 0.0026 0.0027 -0.0010 -0.0049 -0.0035 0.0033 0.0080 0.0034 -0.0074 -0.0119 -0.0018 0.0140 0.0161 -0.0027 -0.0241 -0.0201 0.0127 0.0406 0.0236 -0.0354 -0.0754 -0.0258 0.1214 0.2871 0.3597 0.2871 0.1214 -0.0258 -0.0754 -0.0354 0.0236 0.0406 0.0127 -0.0201 -0.0241 -0.0027 0.0161 0.0140 -0.0018 -0.0119 -0.0074 0.0034 0.0080 0.0033 -0.0035 -0.0049 -0.0010 0.0027 0.0026 0.0000 -0.0017 -0.0011 0.0003 0.0008 0.0003 增益和相位响应如下:
从图中可以看出设计的滤波器满足要求。N=60. Q7.24 用函数kaiserord 和firl 重做习题Q7.23
程序:
clear;
Wp = 0.31;
Ws = 0.41;
As = 50;
Ds = 10^(-As/20);
F = [Wp Ws];
A = [1 0];
DEV = [Ds Ds];
[N,Wn,BTA,Ftype] = kaiserord(F,A,DEV);
Win = kaiser(N+1,BTA);
h = fir1(N,Wn,Ftype,Win);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
参数如下:
Wp = 0.31; Ws = 0.41; As = 50 dB
增益和相位响应如下:
从图中可以看出设计的滤波器满足要求。N=59. Q7.25 用fir2设计一个95阶有限冲激响应滤波器。
程序:
clear;
N = 95;
A = [0.4 0.4 1.0 1.0 0.8 0.8];
F = [0 0.25 0.3 0.45 0.5 1.0];
h = fir2(N,F,A);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
plot(w2/pi,abs(Hz));grid;
xlabel('\omega /\pi' ); ylabel('|H(e^{j\omega})|' ); title('|H(e^{j\omega})|' );
figure(3);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(4);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
幅度响应:
从幅度响应中可以看出,此滤波器满足指标。 Q7.26 使用remez 设计有限冲激响应带通滤波器。 程序:
clear;
F = [1200 1800 3600 4200];
A = [0 1 0];
DEV = [0.02 0.1 0.02];
Fs = 12000;
Dp = 0.1;
Ds = 0.02;
[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs); F2 = 2*[0 1200 1800 3600 4200 6000]/Fs; A2 = [0 0 1 1 0 0];
wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds]; h = firpm(N,F2,A2,wgts);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' ); title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
增益响应: 相位响应:
从增益响应的图像中可以看出,此滤波器满足指标。N=37. Q7.27 用remez 设计具有如下指标的有限冲激响应带通滤波器。
程序:
clear;
Fs1 = 1500;
Fp1 = 1800;
Fp2 = 3000;
Fs2 = 4200;
Fs = 12000;
Dp = 0.1;
Ds = 0.02;
F = [Fs1 Fp1 Fp2 Fs2];
A = [0 1 0];
DEV = [Ds Dp Ds];
[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs);
ws1 = 2*Fs1/Fs;
wp1 = 2*Fp1/Fs;
wp2 = 2*Fp2/Fs;
ws2 = 2*Fs2/Fs;
F2 = [0 ws1 wp1 wp2 ws2 1];
A2 = [0 0 1 1 0 0];
wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds];
h = firpm(N,F2,A2,wgts);
disp('Numerator Coefficients are ' );disp(h);
[g, w] = gain(h,[1]);
figure(1);
plot(w/pi,g);grid;
axis([0 1 -80 5]);
xlabel('\omega /\pi' ); ylabel('Gain in dB' );
title('Gain Response' );
w2 = 0:pi/511:pi;
Hz = freqz(h,[1],w2);
figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\omega /\pi' ); ylabel('Phase (rad)' );
title('Phase Response' );
figure(3);
UPhase = unwrap(Phase);
plot(w2/pi,UPhase);grid;
xlabel('\omega /\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' );
figure(4);
hold on ;
tmpY = 0:1.4/4:1.4;
tmpX = ones(1,length(tmpY))*wp1;
plot(tmpX,tmpY, 'r-' );
tmpX = ones(1,length(tmpY))*wp2;
plot(tmpX,tmpY, 'r-' );
tmpX = ones(1,length(tmpY))*ws1;
plot(tmpX,tmpY, 'g-' );
tmpX = ones(1,length(tmpY))*ws2;
plot(tmpX,tmpY, 'g-' );
tmpY = ones(1,length(w))*(Dp);
plot(w/pi,tmpY, 'r-' );
tmpY = ones(1,length(w))*(Ds);
plot(w/pi,tmpY, 'g-' );
plot(w2/pi,abs(Hz));grid;
hold off ;
增益响应:
该滤波器不满足设计指标。用
估计滤波器的阶数N=73。 从增益响应的图像可以看出,kaiserord