下载文档第一篇：数字信号处理第七章介绍
第七章 数字滤波器设计
7.1：无限冲激响应滤波器的阶数的估计
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.202_.3/使用这些值得到切比雪夫2型低通滤波器最低阶数N=5，相应的标准通带边缘频率Wn是0.4000.4.使用这些值得到椭圆低通滤波器最低阶数N=8，相应的标准通带边缘频率Wn是0.202_.从以上结果中观察到椭圆滤波器的阶数最低，并且符合要求。
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].从以上结果中观察到椭圆滤波器的阶数最低，并且符合要求。
7.2：无限冲激响应滤波器设计
程序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来设计巴特沃兹带阻滤波器。写出所产生的传输函数的准确表达式。滤波器的指标是什么，你的设计符合指标吗，使用MATLAB，计算并绘制滤波器的未畸变的相位响应及群延迟响应。答：表达式是：
滤波器参数是：
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.3：吉布斯现象
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];对于每种情况，研究希尔伯特变换器的频率响应的摆动行为。波纹的数量与希尔伯特变换器的长度之间有什么关系?最大波纹的高度与滤波器的长度之间有什么关系? 答：幅度响应如下：
从中可以观察到由于吉布斯现像产生的幅度响应的摆动行为。波纹的数量与希尔伯特变换器的长度之间的关系——两者成正比。最大波纹的高度与滤波器的长度之间的关系——两者无关系。7.4：有限冲激响应滤波器的阶数估计
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滤波器的阶数估算，其中采样频率改为FT = 20 kHz，则结果为 N=91。(b)线性相位FIR(c)线性相位FIR从上述结果和7.13的对比我们可以观察到：
滤波器阶数和采样频率的关系为–对于一个给定的模拟过渡带宽，采样频率的增加导致估算阶数也相应增加，朝下一个整数取整。
其中模拟过渡带宽|Fp-Fs|和Δω的关系：Δω=2pi*|Fp-Fs|/FT。因此增加FT会减小Δω。
滤波器阶数和通带波纹宽度的关系为估计的阶数大致和log（底数为10）成比例的扩散。滤波器阶数和过渡带宽度的关系为在舍入的时候，阶数随着过渡带宽成比例的改变。有两个因素增加过渡带宽来分割顺序。
Q7.15 线性相位FIR低通滤波器阶数的估算，其中滤波器满足7.13给的规格，使用kaiserord的结果为N=54 正确结果：kaiserord([202_ 2500],[1 0],[0.005 0.005],10000)将上述结果和7.13比较我们观察到：用凯瑟来估算阶数是较小的。因为凯瑟使用了一个不同的近似估计。这种估计经常和FIR设计的凯瑟窗一起用。
Q7.16 线性相位FIR低通滤波器的阶数估算满足的规格和7.13中的一样，使用remezord函数的结果为N=47.正确结果：firpmord([202_ 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。= 0.002p
= 0.002 结果为 N=57。s
s = 2.3 kHz，结果为N=76.使用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时更准确。
7.5有限冲激响应滤波器设计
Q7.20 使用matlab程序设计并画出线性相位FIR滤波器增益和相位反应，使用fir1如下。通过使用函数kaiserord.来估计滤波器阶数，输出结果为滤波器的系数。低通滤波器满足7.20所要求的规格的系数如下：
增益和相位响应如下:
从增益图像我们可以知道这个设计不能满足规格.这个滤波器满足规格的阶数为N=66.为了满足规格，图如下：
Q7.21 汉宁窗：
布莱克曼窗：
切比雪夫窗：
Q7.22 程序如下：
% Program Q7_22 % Use Parks-McClellan to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.13.%Compute and plot the gain function.%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 using Parks-McClellan Wp = 2*Fp/FT;% These freqs are normalized: they go Ws = 2*Fs/FT;% zero to one, not zero to pi.F = [0 Wp Ws 1];A = [1 1 0 0];h = firpm(N,F,A);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;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));T3 = 0.005*ones(1,length(w2));figure(4);plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid;% Find and plot the phase 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 用凯泽窗设计一个有限冲激响应低通滤波器。程序：
% Program Q7_23 % Use Kaiser window to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.% % It is clear from the statement of the question that Mitra % wants us to use(7.36)and(7.37)for this problem.That % isn't the greatest thing to try because kaiserord already does % exactly what we need....but that's Q7_24!So here goes!%Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.23.Wp = 0.31;Ws = 0.41;Wn = Wp +(Ws-Wp)/2;As = 50;Ds = 10^(-As/20);Dp = Ds;%Kaiser window design has equal ripple in % passband and stopband.% estimate order using(7.37)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 % Use(7.36)to get Beta 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);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [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');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase 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 程序：
% Use Kaiser window to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.Use kaiserord and fir1.%Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.23.Wp = 0.31;Ws = 0.41;As = 50;Ds = 10^(-As/20);% Design the Filter 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);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [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');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase 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阶有限冲激响应滤波器。
程序：
% Program Q7_25 % Use fir2 to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.%Print out the numerator coefficients % for the transfer function.%-Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.17.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);N % firpm setup 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);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [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');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase 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设计具有如下指标的有限冲激响应带通滤波器。
程序：
% Program Q7_27 % Use kaiserord and firpm to design the linear phase bandpass % FIR Digital Filter specified in Q7.27.% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.27.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);% firpm setup 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);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [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');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase 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);% Add lines to the plot to help determine if the spec was met.hold on;tmpY = 0:1.4/4:1.4;tmpX = ones(1,length(tmpY))*wp1;plot(tmpX,tmpY, 'r-');% vertical line at passband edge freq tmpX = ones(1,length(tmpY))*wp2;plot(tmpX,tmpY, 'r-');% vertical line at passband edge freq tmpX = ones(1,length(tmpY))*ws1;plot(tmpX,tmpY, 'g-');% vertical line at stopband edge freq tmpX = ones(1,length(tmpY))*ws2;plot(tmpX,tmpY, 'g-');% vertical line at stopband edge freq tmpY = ones(1,length(w))*(Dp);plot(w/pi,tmpY, 'r-');% horizontal line at Dp tmpY = ones(1,length(w))*(Ds);plot(w/pi,tmpY, 'g-');% horizontal line at Ds % now plot the Frequency Response plot(w2/pi,abs(Hz));grid;hold off;增益响应：
从增益响应的图像可以看出，该滤波器不满足设计指标。用kaiserord估计滤波器的阶数N=73。
第二篇：数字信号处理第四章介绍
第四章 线性时不变离散时间系统的频域分析
一、传输函数和频率响应 例4.1传输函数分析 Q4.1 clear;M = input('Enter the filter length M: ');w = 0:2*pi/1023:2*pi;num =(1/M)*ones(1,M);den = [1];h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,abs(h));grid title('Magnitude Spectrum |H(e^{jomega})|')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));grid title('Phase Spectrum arg[H(e^{jomega})]')xlabel('omega /pi');ylabel('Phase in radians');
M=2
M=10
M=15
幅度谱为偶对称，相位谱为奇对称，这是一个低通滤波器。M越大，通带越窄且过渡带越陡峭。
Q4.2使用修改后的程序P3.1，计算并画出当w=[0,pi]时传输函数因果线性时不变离散时间系统的频率响应。它表示哪种类型的滤波器？ w = 0:pi/511:pi;
num = [0.15 0-0.15];den = [1-0.5 0.7];如下图1这是一个带通滤波器。的图1
图2 Q4.3对下面的传输函数重做习题Q4.2：，式（4.36）和式（4.37）给出的两个滤波器之间的区别是什么？你将选择哪一个滤波器来滤波，为什么？ w = 0:pi/511:pi;num = [0.15 0-0.15];den = [0.7-0.5 1];如上图2也是一个带通滤波器，这两个滤波器的幅度谱是一样的，相位谱不太一样，我会选择第一个带通滤波器，因为它的相位谱更加平滑，相位失真小。
Q4.4 使用MATLAB计算并画出当w=[0,pi]时因果线性时不变离散时间系统的群延迟。系统的传输函数为clf;w = 0:pi/511:pi;num = [1-1.2 1];den = [1-1.3 1.04-0.222];h= grpdelay(num,den,w);plot(w/pi,h);xlabel('w/pi');ylabel('群延迟')。
Q4.5 使用Q3.50中编写的程序，分别计算并画出式（4.36）和式（4.37）确定的两个滤波器的冲激响应中的前一百个样本。讨论你的结果。clf;num = [0.15 0-0.15];den = [0.7-0.5 1];L = input('输入样本数 L: ');[g t] = impz(num,den,L);stem(t,g);title(['前 ',num2str(L),' 脉冲响应的样本']);xlabel('时间序号 n');ylabel('h[n]');
（4.36）式（4.37）式
由图可知：这些情节由impz给生成的因果的脉冲响应实现的H(z)。我们观察到Q4.3因果滤波器与H(z)在(4.36)稳定,这意味着H[n]是绝对可和,我们看到交替和指数衰减的脉冲响应。在另一方面,因果编档人员与H(z)在(4.37)极点以外的单位圆,是不稳定的。不足为奇的是,相应的h[n]上图显示与n指数增长。
Q4.6 传输函数的极零点图同样能分析线性时不变离散时间系统的性质。使用命令zplane可以很容易地得到系统的极零点图。使用zplane分别生成式（4.36）和式（4.37）确定的两个滤波器的极零点图。讨论你的结果。clf;num = [0.15 0-0.15];den = [1-0.5 0.7];[z p k] = tf2zpk(num,den);disp('Zeros:');disp(z);disp('Poles:');disp(p);input('Hit
式（4.37）
由图可知：过滤器在(4.36)在单位圆和两极因此它的因果实现稳定；较低的图显示过滤器(4.37)极点在单位圆外,其因果关系的实现是不稳定的。
二、传输函数的类型 例4.2滤波器 Q4.7 clf;fc = 0.25;n = [-6.5:1:6.5];y = 2*fc*sinc(2*fc*n);k = n+6.5;stem(k,y);title('N = 14');axis([0 13-0.2 0.6]);xlabel('Time index n');ylabel('Amplitude');grid;
图1 图2 如图1低通有限冲激滤波器的长度为14，决定滤波器长度的语句为n = [-6.5:1:6.5]，而控制截止频率的参数是fc = 0.25。Q4.8 fc = 0.45;n = [-9.5:1:9.5];y = 2*fc*sinc(2*fc*n);k = n+9.5;stem(k,y);title('N = 20');axis([0 19-0.2 0.7]);xlabel('Time index n');ylabel('Amplitude');grid;修改参数fc和n，得到如上图2，可知低通有限冲激滤波器的长度变为20.Q4.9 clf;fc = 0.65;n = [-7.0:1:7.0];y = 2*fc*sinc(2*fc*n);k = n+7.0;stem(k,y);title('N = 14');axis([0 14-0.4 1.4]);xlabel('Time index n');ylabel('Amplitude');grid;
Q4.10 clear;N = input('Enter the filter time shift N: ');No2 = N/2;fc = 0.25;n = [-No2:1:No2];y = 2*fc*sinc(2*fc*n);w = 0:pi/511:pi;h = freqz(y, [1], w);plot(w/pi,abs(h));grid;title(strcat('|H(e^{jomega})|, N=',num2str(N)));xlabel('omega /pi');ylabel('Amplitude');
上图依次分别为N=5,10,30,100的四幅图，从这四幅图可以看出随着阶数N的增大，低通滤波器的过渡带越来越窄，阻带衰减越来越快，滤波器越来越接近理想低通滤波器。Q4.11 clf;M = 2;num = ones(1,M)/M;[g,w] = gain(num,1);plot(w/pi,g);grid axis([0 1-50 0.5])xlabel('omega /pi');ylabel('Gain in dB');title(['M = ',num2str(M)])
可以验证3dB截止频率在π/2处。Q4.12 clear;K = input('Enter the number of sections K: ');Hz = [1];for i=1:K;Hz = conv(Hz,[1 1]);end;Hz =(0.5)^K * Hz;[g,w] = gain(Hz,1);ThreedB =-3*ones(1,length(g));t1 = 2*acos((0.5)^(1/(2*K)))*ones(1,512)/pi;t2 =-50:50.5/511:0.5;plot(w/pi,g,w/pi,ThreedB,t1,t2);grid;axis([0 1-50 0.5])xlabel('omega /pi');ylabel('Gain in dB');title(['K = ',num2str(K),';Theoretical omega_{c} = ',num2str(t1(1))]);
Q4.13 clear;M = input('Enter the filter length M: ');n = 0:M-1;num =(-1).^n.* ones(1,M)/M;[g,w] = gain(num,1);plot(w/pi,g);grid;axis([0 1-50 0.5]);xlabel('omega /pi');ylabel('Gain in dB');title(['M = ', num2str(M)]);
其3dB截止频率约为0.82pi Q4.14 设计一个在0.45pi处具有3dB截止频率wc的一阶无限冲激响应低通滤波器和一阶无限冲激响应高通滤波器。用MATLAB计算并画出它们的增益响应，验证设计的滤波器是否满足指标。用MATLAB证明两个滤波器是全通互补和功率互补的。
Q4.15 级联10个式（4.15）所示一阶无限冲激响应低通滤波器，设计一个在0.3pi处具有3dB截止频率wc的无限冲激响应低通滤波器。把它与一个具有相同截止频率的一阶无限冲激响应低通滤波器的增益响应作比较。
Q4.16 设计一个中心频率wo在0.61pi处、3dB带宽为0.51pi的二阶无限冲激响应带通滤波器。由于式（4.20）是α的二次方程，为了产生相同的3dB带宽，参数α将有两个数值，得到的传输函数HBP（z）也会有两个不同的表达式。使用函数zplane可产生两个传输函数的极零点图，从中可以选择一个稳定的传输函数。用MATLAB计算并画出你所设计的滤波器的增益响应，并验证它确实满足给定的条件。用设计的稳定无限冲激响应带通滤波器的传输函数的参数α和β，生成一个二阶无限冲激响应带阻滤波器的传输函数HBS（z）。用MATLAB证明HBP(z)和HBS(z)都是全通互补和功率互补的。
Q4.17 用MATLAB计算并画出一个梳状滤波器的幅度响应，该梳状滤波器是在L取不同值的情况下，由式（4.40）给出的原型有限冲激响应低通滤波器得到的。证明新滤波器的幅度响应在k=0,1,2,3......，L-1.处有L个极小值，在处有L个极大值，Q4.18 用MATLAB计算并画出一个梳状滤波器的幅度响应，该梳状滤波器是在L取不同值的情况下，由式（4.42）在M=2时给出的原型有限冲激响应低通滤波器得到的。确定这种梳状滤波器冲激响应的极大值和极小值的位置。
从这些情节我们观察,梳状滤波器极距为1kπ/L,山峰为(2k+1)π/L.Q4.19 clf;b = [1-8.5 30.5-63];num1 = [b 81 fliplr(b)];num2 = [b 81 81 fliplr(b)];num3 = [b 0-fliplr(b)];num4 = [b 81-81-fliplr(b)];n1 = 0:length(num1)-1;n2 = 0:length(num2)-1;subplot(2,2,1);stem(n1,num1);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 1 FIR Filter');subplot(2,2,2);stem(n2,num2);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 2 FIR Filter');subplot(2,2,3);stem(n1,num3);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 3 FIR Filter');subplot(2,2,4);stem(n2,num4);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 4 FIR Filter');pause subplot(2,2,1);zplane(num1,1);title('Type 1 FIR Filter');subplot(2,2,2);zplane(num2,1);title('Type 2 FIR Filter');subplot(2,2,3);zplane(num3,1);title('Type 3 FIR Filter');subplot(2,2,4);zplane(num4,1);title('Type 4 FIR Filter');disp('Zeros of Type 1 FIR Filter are');disp(roots(num1));disp('Zeros of Type 2 FIR Filter are');disp(roots(num2));disp('Zeros of Type 3 FIR Filter are');disp(roots(num3));disp('Zeros of Type 4 FIR Filter are');disp(roots(num4));
1型有限冲激响应滤波器的零点是 Zeros of Type 1 FIR Filter are 2.9744 2.0888 0.9790 + 1.4110i 0.97900.4784i 0.4787 0.3362 2型有限冲激响应滤波器的零点是 Zeros of Type 2 FIR Filter are 3.7585 + 1.5147i 3.75852.6623i-1.0000 0.0893 + 0.3530i 0.08930.0922i 3型有限冲激响应滤波器的零点是 Zeros of Type 3 FIR Filter are 4.7627 1.6279 + 3.0565i 1.62790.2549i 0.2100 4型有限冲激响应滤波器的零点是 Zeros of Type 4 FIR Filter are 3.4139 1.6541 + 1.5813i 1.65410.9973i 1.0000 0.3159 + 0.3020i 0.31592.0140i-1.2659 + 2.0135i-1.26590.3559i 0.2457 + 0.2126i 0.24572.0392i-1.0101 + 2.1930i-1.01010.3762i 0.2397 + 0.1934i 0.23971.2263i 1.0000 0.6576 + 0.7534i 0.65760.7803i 4型有限冲激响应滤波器的零点是 Zeros of Type 4 FIR Filter are 2.0841 + 2.0565i 2.08411.9960i 1.0000-0.2408 + 0.3197i-0.24080.2399i Q4.21 用MATLAB 确定如下传输函数是否是有界实函数：一个有界实函数，求一个与
有着相同幅度的有界实函数。
它若不是由下图可知：H1(z)不是有界实函数。
故H2(z)为
Q4.22 用MATLAB 确定如下传输函数是否是有界实函数：有界实函数，求一个与
有着相同幅度的有界实函数。
它若不是一个使用zplane我们观察到的G1(z)在单位圆,因此传递函数是稳定的。
Q4.23 用MATLAB产生如下两个因果系统传输函数的极零点图：，研究生成的极零点图，你可以推断它们的稳定性么？
用Q4.6的程序做H1(z)
，Q4.24 用程序P4.4检测Q4,23中两个传输函数的稳定性。这两个传输函数哪一个是稳定的？ % Program P4_4 clf;den = input('分母系数 = ');ki = poly2rc(den);disp('稳定性测试参数是');disp(ki);
由此我们可以总结出H1(z)稳定，H2(z)不稳定
Q4.25 用程序P4.4确定下面这个多项式的所有根是否都在单位圆内：
由此看出，都在单位圆内。
Q4.26 用程序P4.4确定下面这个多项式的所有根是否都在单位圆内：
由此看出，都在单位圆内。
第三篇：数字信号处理第六章介绍
第六章
数字滤波器结构
6.1：级联的实现
num = input('分子系数向量 = ');den = input('分母系数向量 = ');[z,p,k] = tf2zp(num,den);sos = zp2sos(z,p,k)Q6.1使用程序P6.1，生成如下有限冲激响应传输函数的一个级联实现: H1(z)=2+10z^(-1)+23z^(-2)+34z^(-3)+31z^(-4)+16 z^(-5)+4z^(-6)画出级联实现的框图。H1(z)是一个线性相位传输函数吗? 答：运行结果：
sos = zp2sos(z,p,k)Numerator coefficient vector = [2,10,23,34,31,16,4] Denominator coefficient vector = [1] sos = 2.0000 6.0000 4.0000 1.0000 0 0 1.0000 1.0000 2.0000 1.0000 0 0 1.0000 1.0000 0.5000 1.0000 0 0
级联框图：
H1（z）不是一个线性相位传输函数，因为系数不对称。
Q6.2使用程序P6.1，生成如下有限冲激响应传输函数的一个级联实现: H2(z)=6+31z^(-1)+74z^(-2)+102z^(-3)+74z^(-4)+31 z^(-5)+6z^(-6)画出级联实现的框图。H2(z)是一个线性相位传输函数吗?只用4个乘法器生成H2(z)的一级联实现。显示新的级联结构的框图。
Numerator coefficient vector = [6,31,74,102,74,31,6] Denominator coefficient vector = [1] sos = 6.0000 15.0000 6.0000 1.0000 0 0 1.0000 2.0000 3.0000 1.0000 0 0 1.0000 0.6667 0.3333 1.0000 0 0 级联框图：
H2（z）是一个线性相位传输函数。只用四个乘法器生成级联框图：
6.2：级联和并联实现
Q6.3使用程序P6.1生成如下因果无限冲激响应传输函数的级联实现： 画出级联实现的框图。答：
Numerator coefficient vector = [3,8,12,7,2,-2] Denominator coefficient vector = [16,24,24,14,5,1] sos = 0.1875-0.0625 0 1.0000 0.5000 0 1.0000 2.0000 2.0000 1.0000 0.5000 0.2500 1.0000 1.0000 1.0000 1.0000 0.5000 0.5000
级联实现框图：
Q6.4使用程序P6.1生成如下因果无限冲激响应传输函数的级联实现：
画出级联实现的框图。答：级联实现框图：
程序P6.2生成两种类型的并联实现 num = input('分子系数向量 = ');den = input('分母系数分量 = ');[r1,p1,k1] = residuez(num,den);[r2,p2,k2] = residue(num,den);disp('并联I型')disp('留数是');disp(r1);disp('极点在');disp(p1);disp('常数');disp(k1);disp('并联II型')disp('留数是');disp(r2);disp('极点在');disp(p2);disp('常数');disp(k2);Q6.5使用程序P6.2生成式(6.27)所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图。答：并联I型框图：
并联II型框图：
Q6.6使用程序P6.2生成式(6.28)所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图。答：并联I型框图：
并联II型框图：
6.3：全通传输函数的实现
Q6.7使用程序P4.4生成如下全通传输函数的级联格型实现:
As(z)是一个稳定的传输函数吗? 答：
运行结果：
k(5)= 0.0625 k(4)= 0.2196 k(3)= 0.4811 k(2)= 0.6837 k(1)= 0.6246
2从{ki}的值我们可以得到传输函数A5(z)是稳定的，因为对所有的1 Q6.8使用程序P4.4生成如下全通传输函数的级联格型实现: A6(z)足一个稳定的传输函数吗? 答：得到A6(z)的{ki}值如下： k(6)= 0.0278 k(5)= 0.1344 k(4)= 0.3717 k(3)= 0.5922 k(2)= 0.7711 k(1)= 0.8109 从{ki}的值可以得到传输函数A6(z)是稳定的，因为反馈系数的平均幅值小于整体。 Q6.9 使用l型和2型全通项生成式(6.29)所示全通传输函数的典范级联实现。显示实现的框图。在最终的结构中，乘法器的总数是多少? 答：全通因子如下所示： 使用1型和2型全通项生成所示全通函数的典范级联实现，实现的结构框图如下： 整体结构中乘法器的总数是5.Q6.10 用zp2sos 我们可以得到 A6(z)的因子如下： sos = 0.0278 0.0556 0.1111 1.0000 0.5000 0.2500 1.0000 2.0000 3.0000 1.0000 0.6667 0.3333 1.0000 3.0000 3.0000 1.0000 1.0000 0.3333 从上面因子可以分解 A6(z)为低阶的全通因子: 使用2型的全通项生成A6(z)的典范级联实现框图如下： 整体结构中乘法器的总数是6。6.4：无限冲激响应传输函数的Gary-Markel实现 num = input('分子系数向量 = ');den = input('分母系数向量 = ');N = length(den)-1;% 分母多项式的阶数 k = ones(1,N);a1 = den/den(1);alpha = num(N+1:-1:1)/den(1);for ii = N:-1:1, alpha(N+2-ii:N+1)= alpha(N+2-ii:N+1)-alpha(N-ii+1)*a1(2:ii+1);k(ii)= a1(ii+1);a1(1:ii+1)=(a1(1:ii+1)-k(ii)*a1(ii+1:-1:1))/(1-k(ii)*k(ii));end disp('格型参数是');disp(k)disp('前馈乘法器是');disp(alpha)Q6.11 使用程序 P6_3我们通过IIR将Q6.3给的正向传输函数H1(z)的 Gray-Markel级联格型实现参数如下： 晶格参数和前馈乘数分别如下： 对应Gray-Markel的结构框图如下： 使用程序P6_3，从这些格型参数可以得到传输函数H1(z)是稳定的，因为所有格型参数的平方值比整体的小。 Q6.12 使用程序 P6_3我们通过IIR将Q6.4给的正向传输函数H2(z)的 Gray-Markel级联格型实现参数如下： 对应Gray-Markel的结构框图如下： 使用程序P6_3，从这些格型参数可以得到传输函数H2(z)是稳定的，因为所有格型参数的平方值比整体的小。 Q6.13使用函数tf2latc编写出一个MATLAB程序，以生成一个因果无限冲激响应传输函数的GrayMarkel实现。用该程序实现式(6.27)所示的传输函数。你的结果与习题6.11中得到的结果相符吗？使用函数1atc2tf由向量k和alpha确定传输函数。所得到的传输函数和式(6.27)给出的传输函数相同吗？ 答：程序如下： format long num = input('Numerator coefficient vector = ');den = input('Denominator coefficient vector = ');num = num/den(1);% normalize upstairs and down by d0.den = den/den(1);% here is the lattice/ladder realization from the transfer fcn: [k,alpha] = tf2latc(num,den)% now check inversion disp('Check of Lattice/Ladder Inversion:');[num2,den2] = latc2tf(k,alpha)运行结果如下： k = 0.62459686089013 0.68373782742919 0.48111942348398 0.21960784313725 0.06250000000000 alpha =-0.01982100623522-0.09085169508677 0.***849 0.16053921568627 0.31250000000000-0.12500000000000 结果与习题6.11中得到的结果相符。Q6.14使用在习题6.13中生成的程序，实现式(6.28)给出的传输函数。你的结果与习题6.12中得到的结果相符吗？使用函数latc2tf由向量k和alpha确定传输函数。所得到的传输函数和式(6.28)给出的传输函数相同吗? 答：运行结果： k = 0.81093584641352 0.77112772506402 0.592*** 0.37169052478550 0.***293 0.02777777777778 alpha =-0.01112037033486 0.02345313662512-0.0***79-0.04739265773254 0.***485 0.20370370370370 0.11111111111111 与题6.12中得到的结果相符。 6.5：无限冲激响应传输函数的并联全通实现 Q6.15 生成下式给出的只阶因果有界实低通1型切比雪夫传输函数G(z)的全通和的分解。使用zplane获得G(z)的零极点分布图： G(z)全通和的分解： G(z)的功率补充传输函数H(z)的表达式如下： 两个全通传输函数的阶数是1和2.Q6.15 生成一个五阶因果有界实低通椭圆传输函数G(z)的全通和的分解。使用zplane获得G(z)的零极点分布图： G(z)全通和的分解： G(z)的功率补充传输函数H(z)的表达式如下： 两个全通传输函数的阶数是3和2. 目 录 摘要...........................................................................................................................................1 1 绪论..............................................................................................................................................2 1.1 DSP系统特点和设计基本原则......................................................................................2 1.2 国内外研究动态.............................................................................................................2 2系统设计........................................................................................................................................3 3硬件设计........................................................................................................................................5 3.1 硬件结构...........................................................................................................................5 3.2 硬件电路设计...................................................................................................................7 3.2.1 总输入电路...........................................................................................................7 3.2.2 总输出电路...........................................................................................................7 3.2.3 语音输入电路.......................................................................................................9 3.2.4 语音输出电路.......................................................................................................9 实验结果及分析.........................................................................................................................10 4.1 实验结果.........................................................................................................................10 4.2 实验分析.........................................................................................................................12 5 总结与心得体会.........................................................................................................................13 参考文献.........................................................................................................................................14 致谢................................................................................................................................................15 摘要 基于DSP的语音信号处理系统，该系统采用TMS320VC5509作为主处理器，TLV320AIC23B作为音频芯片，在此基础上完成系统硬件平台的搭建和软件设计，从而实现对语音信号的采集、滤波和回放功能，它可作为语音信号处理的通用平台。 语音是人类相互之间进行交流时使用最多、最自然、最基本也是最重要的信息载体。在高度信息化的今天，语音信号处理是信息高速公路、多媒体技术、办公自动化、现代通信及智能系统等新兴领域应用的核心技术之一。通常这些信号处理的过程要满足实时且快速高效的要求，随着DSP技术的发展，以DSP为内核的设备越来越多，为语音信号的处理提供了良好的平台。本文设计了一个基于TMS320VC5509定点的语音信号处理系统，实现对语音信号的采集、处理与回放等功能，为今后复杂的语音信号处理算法的研究和实时实现提供一个通用平台。 关键词：语音处理；DSP；TMS320VC5509;TLV320AIC23B 1 绪论 语音是人类相互间所进行的通信的最自然和最简洁方便的形式，语音通信是一种理想的人机通信方式。语音通信的研究涉及到人工智能、数字信号处理、微型计算机技术、语言声学、语言学等许多领域，所以说语音的通信是一个多学科的综合研究领域，其研究成果具有重要的学术价值。另外通过语音来传递信息是人类最重要的、最有效、最常用的交换信息的形式。语言是人类特有的功能，声音是人类常用的工具，是相互传递信息的主要手段。同时也是众构成思想交流和感情沟通的最主要的途径。 1.1 DSP系统特点和设计基本原则 DSP（digital signal processor）是一种独特的微处理器，是以数字信号来处理大量信息的器件。其工作原理是接收模拟信号，转换为0或1的数字信号。再对数字信号进行修改、删除、强化，并在其他系统芯片中把数字数据解译回模拟数据或实际环境格式。它不仅具有可编程性，而且其实时运行速度可达每秒数以千万条复杂指令程序，远远超过通用微处理器，是数字化电子世界中日益重要的电脑芯片。它的强大数据处理能力和高运行速度，是最值得称道的两大特色。 1.2 国内外研究动态 语音信号处理作为一个重要的研究领域，已经有很长的研究历史。但是它的快速发展可以说是从1940年前后Dudley的声码器和Potter等人的可见语音开始的；20世纪60年代中期形成的一系列数字信号处理的理念和技术基础；到了80年代，由于矢量量化、隐马尔可夫模型和人工神经网络等相继被应用于语音信号处理，并经过不断改进与完善，使得语音信号处理技术产生了突破性的进展。一方面，对声学语音学统计模型的研究逐渐深入，鲁棒的语音识别、基于语音段的建模方法及隐马尔可夫模型与人工神经网络的结合成为研究的热点。另一方面，为了语音识别实用化的需要，讲者自适应、听觉模型、快速搜索识别算法以及进一步的语言模型的研究等课题倍受关注。 在通信越来越发达的当今世界，尤其最近几十年，语音压缩编码技术在移动 通信、IP电话通信、保密通信、卫星通信以及语音存储等很多方面得到了广泛的应用。因此，语音编码一直是通信和信号处理的研究热点，并其取得了惊人的进展，目前在PC机上的语音编码已经趋于成熟，而如何在嵌入式系统中实时实现语音压缩编码则是近些年来语音信号处理领域的研究热点之一。 2系统设计 在实际生活中，当声源遇到物体时会发生反射，反射的声波和声源声波一起传输，听者会发现反射声波部分比声源声波慢一些，类似人们面对山体高声呼喊后可以在过一会儿听到回声的现象。声音遇到较远物体产生的反射会比遇到较近的反射波晚些到达声源位置，所以回声和原声的延迟随反射物体的距离大小改变。同时，反射声音的物体对声波的反射能力，决定了听到的回声的强弱和质量。另外，生活中的回声的成分比较复杂，有反射、漫反射、折射，还有回声的多次反射、折射效果。 当已知一个数字音源后，可以利用计算机的处理能力，用数字的方式通过计算模拟回声效应。简单的讲，可以在原声音流中叠加延迟一段时间后的声流，实现回声效果。当然通过复杂运算，可以计算各种效应的混响效果。如此产生的回声，我们称之为数字回声。 本次实验的程序流程图如下： 图2.1 程序流程图 本次实验的系统框图如下： 图2.2 系统框图 3硬件设计 3.1 硬件结构 图3.1是系统的硬件结构框图, 系统主要包括VC5509和A IC23 两个模块。 图3.1系统硬件结构框图 利用VC5509 的片上外设I2C(Inter-Integrated Circuit, 内部集成电路)模块配置AIC23 的内部寄存器;通过VC5509 的McBSP(Multi channel Buffered Serial Ports, 多通道缓存串口)接收和发送采样的音频数据。控制通道只在配置AIC23 的内部寄存器时工作, 而当传输音频数据时则处于闲置状态。 AIC23通过麦克风输入或者立体声音频输入采集模拟信号, 并把模拟信号转化为数字信号, 存储到DSP的内部RAM中,以便DSP处理。 当DSP完成对音频数据的处理以后, AIC23再把数字信号转化为模拟信号, 这样就能够在立体声输出端或者耳机输出端听到声音。 AIC23能够实现与VC5509 DSP的McBSP端口的无缝连接, 使系统设计更加简单。接口的原理框图, 如下图所示。 图3.2 AIC23与VC5509接口原理图 系统中A IC23的主时钟12 MHz直接由外部的晶振提供。MODE接数字地, 表示利用I2 C控制接口对AIC23传输控制数据。CS接数字地, 定义了I2 C总线上AIC23的外设地址, 通过将CS接到高电平或低电平, 可以选择A IC23作为从设备在I2 C总线上的地址。SCLK和SDIN是AIC23控制端口的移位时钟和数据输入端,分别与VC5509的I2C模块端口SCL和SDA相连。 收发时钟信号CLKX1和CLKR1由A IC23的串行数据输入时钟BCLK提供, 并由A IC23的帧同步信号LRCIN、LRCOUT启动串口数据传输。DX1和DR1分别与A IC23 的D IN 和DOUT 相连, 从而完成VC5509与AIC23间的数字信号通信。 3.2 硬件电路设计 3.2.1 总输入电路 图3.3 总输入电路 从左到右各部分电路为： 话筒，开关，语音输入电路，UA741高增益放大电路，有源二阶带 通滤波器。 3.2.2 总输出电路 图3.4 总输出电路 从左到右各部分电路为： LM386高频功率放大器及其外围器件连接电路，语音输出电路，开关，扬声器。 3.2.3 语音输入电路 图3.5语音输入电路 3.2.4 语音输出电路 图3.6 语音输出电路 语音信号通道包括模拟输入和模拟输出两个部分。模拟信号的输入输出电路如图所示。上图中MICBIAS 为提供的麦克风偏压，通常是3/4 AVDD，MICIN为麦克风输入，可以根据需要调整输入增益。下图中LLINEOUT 为左声道输出，RLINEOUT为右声道输出。用户可以根据电阻阻值调节增益的大小，使语音输入输出达到最佳效果。从而实现良好的模拟语音信号输入与模拟信号的输出。4 实验结果及分析 4.1 实验结果 按“F5”键运行，注意观察窗口中的bEcho=0，表示数字回声功能没有激活。这时从耳机中能听到麦克风中的输入语音放送。将观察窗口中bEcho的取值改成非0值。这时可从耳机中听到带数字回声道语音放送。 分别调整uDelay和uEffect的取值，使他们保持在0-1023范围内，同时听听耳机中的输出有何变化。 当uDelay和uEffect的数值增大时，数字回声的效果就会越加的明显。 图4.1 修改前程序图 图4.2 修改前程序图 图4.3 频谱分析 图4.4 左声道及右声道波形 4.2 实验分析 所以，从本实验可知当已知一个数字音源后，可以利用计算机的处理能力，用数字的方式通过计算模拟回声效应。简单的讲，可以在原声音流中叠加延迟一段时间后的声流，实现回声效果。当然通过复杂运算，可以计算各种效应的混响效果。 声音放送可以加入数字回声，数字回声的强弱和与原声的延迟均可在程序中设定和调整。5 总结与心得体会 通过本次课程设计，我明白了细节决定成败这句话的道理，在实验中，有很多注意的地方，都被忽视了，导致再花费更多的时间去修改，这严重影响了试验的进度。同时，在本次实验中我了解了ICETEK – VC5509 – A板上语音codec芯片TLV320AIC23的设计和程序控制原理，并进一步掌握了数字回声产生原理、编程及其参数选择、控制，以及了解了VC5509DSP扩展存储器的编程使用方法。 这一学期的理论知识学习加上这次课程设计，使我对DSP有了更加深刻的了解，对数字信号的处理功能，软硬件相结合，语音信号的采集与放送等等方面都有了很深的了解，相信本次课程设计，无论是对我以后的学习，还是工作等方面都有一个很大的帮助。因此，本次课程设计让我受益匪浅。 参考文献 [1]李利.DSP原理及应用[M].北京:中国水利水电出版社,202_.[2]王安民,陈明欣,朱明.TMS320C54xxDSP实用技术[M].北京:清华大学出版社,202_ [3]彭启琮,李玉柏.DSP技术[M].成都:电子科技大学出版社,1997 [4]李宏伟,等.基于帧间重叠谱减法的语音增强方法[J].解放军理工大学学报,202_(1):41~44 [5]TexasInstrumentsIncorporated.TMS320C54x系列DSP的CPU与外设[M].梁晓雯,裴小平,李玉虎,译.北京:清华大学出版社,202_ [6]赵力.语音信号处理[M].北京:机械工业出版社,202_比较图4和图5,可以看到1200Hz以上的频谱明显得到了抑制。 [7]江涛,朱光喜.基于TMS320VC5402的音频信号采集与系统处理[J].电子技术用,202_,28(7):70~72[8]TexasInstrumentsIncorporated:TMS320VC5402Datasheet,202_ 致谢 在本次课程设计的即将完成之际，笔者的心情无法平静，本文的完成既是笔者孜孜不倦努力的结果，更是指导老师樊洪斌老师亲切关怀和悉心指导的结果。在整个课程设计的选题、研究和撰写过程中，老师都给了我精心的指导、热忱的鼓励和支持，他的精心点拨为我开拓了研究视野，修正了写作思路，对课程设计的完善和质量的提高起到了关键性的作用。另外，导师严谨求实的治学态度、一丝不苟的工作作风和高尚的人格魅力，都给了学生很大感触，使学生终生受益。在此，学生谨向老师致以最真挚的感激和最崇高的敬佩之情。 另外，还要感谢这段时间来陪我一起努力同学，感谢我们这个小团队，感谢每一个在学习和生活中所有给予我关心、支持和帮助的老师和同学们，几年来我们一起学习、一起玩耍，共同度过了太多的美好时光。我们始终是一个团结、友爱、积极向上的集体。 JIANGSU UNIVERSITY OF TECHNOLOGY 数字信号处理实验报告 学院名称： 电气信息工程学院 专 业： 班 级： 姓 名： 学 号： 指导老师： 张维玺（教授） 202_年12月20日 实验一 离散时间信号的产生 一、实验目的 数字信号处理系统中的信号都是以离散时间形态存在的，所以对离散时间信号的研究是数字信号的基本所在。而要研究离散时间信号，首先需要产生出各种离散时间信号。使用MATLAB软件可以很方便地产生各种常见的离散时间信号，而且它还具有强大绘图功能，便于用户直观地处理输出结果。 通过本实验，学生将学习如何用MATLAB产生一些常见的离散时间信号，实现信号的卷积运算，并通过MATLAB中的绘图工具对产生的信号进行观察，加深对常用离散信号和信号卷积和运算的理解。 二、实验原理 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用x（n）来表示，自变量必须是整数。常见的离散信号如下：（1）单位冲激序列δ（n） 如果δ（n）在时间轴上延迟了k个单位，得到δ（n-k），即长度为N的单位冲激序列δ（n）可以通过下面的MATLAB命令获得。 n=-（N-1）：N-1 x=[zeros(1,N-1)1 zeros(1,N-1)]； stem(n,x)延迟K个采样点的长度为N的单位冲激序列δ（n-k）（k n=0：N-1 y=[zeros(1,M)1 zeros(1,N-M-1)]； stem（n，y） （2）单位阶跃序列u（n） 如果u（n）在时间轴上延迟了k个单位，得到u（n-k），即长度为N的单位阶跃序列u（n）可以通过下面的MATLAB命令获得。 n=-（N-1）：N-1 x=[zeros(1,N-1)ones(1,N)]； stem(n,x)延迟的单位阶跃序列可以使用类似于单位冲激序列的方法获得。（3）矩形序列 矩形序列有一个重要的参数，就是序列的宽度N。矩形序列与u（n）之间的关系为矩形序列等= u（n）— u（n-N）。 因此，用MATLAB表示矩形序列可利用上面的单位阶跃序列组合而成。（4）正弦序列x（n） 这里，正弦序列的参数都是实数。与连续的正弦信号不同，正弦序列的自变量n必须为整数。可以证明，只有当2π/w为有理数时，正弦序列具有周期性。 长度为N的正弦序列x（n）可以通过下面的MATLAB命令获得。n=0：N-1 x=A*cos（2*pi*f*n/Fs+phase）（5）单边实指数序列x（n） 长度为N的实指数序列x（n）可以通过下面的MATLAB命令实现。n=0：N-1 x=a.^n stem(n,x)单边指数序列n的取值范围为n>=0。当|a|>1时，单边指数序列发散；当|a|<1时，单边指数序列收敛。当a>0时，该序列均取正值；当a<0时，序列在正负摆动。 （6）负指数序列x（n） 当a=0时，得到虚指数序列x（n）。 与连续负指数信号一样，我们将负指数序列实部和虚部的波形分开讨论，得到如下结论： 1）当a>0时，负指数序列x（n）的实部和虚部分别是按指数规律增长的正弦振荡序列； 2)当a<0时，负指数序列x（n）的实部和虚部分别是按指数规律衰减的正弦振荡序列； 3)当a=0时，负指数序列x（n）即为虚指数序列，其实部和虚部分别是等幅的正弦振荡序列； 长度为N的实指数序列x（n）可以通过下面的MATLAB命令实现。n=0：N-1 x=exp（（a.+j*w）*n）stem(n,real(x))或 stem(n,imag(x)) 三、实验内容及分析 1n01、编制程序产生单位冲激序列n“并绘出其图及n”学号后两位0n0形。程序:(1)N=4; n=-(N-1):N-1; x=[zeros(1,N-1)1 zeros(1,N-1)];stem(n,x); title('单位冲激序列'); grid on; (2)N=6; M=1;%学号01 n=-(N-1):N-1; y=[zeros(1,N-M+1)1 zeros(1,N-M-1)];stem(n,y); title('单位冲激序列');grid on; 分析:在上图的基础上向右平移了1个单位。 1n02、编制程序产生单位阶跃序列un、un“学号后两位”及 0n0unun“学号后两位”，并绘出其图形。程序: 4 (1)N=5; n=-(N-1):N-1; x=[zeros(1,N-1)ones(1,N)];stem(n,x); title('单位阶跃序列');grid on; (2)N=6; M=1;%学号01 n=-(N-1):N-1; x=[zeros(1,N-M+1)ones(1,N-M)];stem(n,x); title('单位阶跃序列');grid on; 分析:在上图的基础上平移了1个单位.(3)N=6; M=1;%学号01 n=-(N-1):N-1; x=[zeros(1,N-1)ones(1,N)];y=[zeros(1,N-M+1)ones(1,N-M)];z=x-y;stem(n,z); title('单位阶跃序列');grid on; 2 3、编制程序产生正弦序列xncos2n、xncosn及 学号后两位xnsin2n并绘出其图形。 程序:(1)N=5; A=1; w=2*pi;phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信号');grid on; 分析:该序列具有周期性,且输出为余弦信号.(2)N=5; A=1; w=2*pi/1;%学号01 phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信号');grid on; ; 分析:该序列具有周期性,且输出为余弦信号.(3)N=5; A=1; w=2*pi;phi=0; n=0:0.05:N-1;x=A*sin(w*n+phi);stem(n,x);title('正弦信号');grid on; 分析:该序列具有周期性,且输出为正弦信号.4、编制程序产生复正弦序列xne(2j学号后两位）n，并绘出其图形。N=3; n=0:0.2:N-1; w=1;%学号01 x=exp((2+j*w)*n);subplot(2,1,1) stem(n,real(x)),title('实部');grid on;subplot(2,1,2) stem(n,imag(x)),title('虚部');grid on; 5、编制程序产生指数序列xnan，并绘出其图形。其中a=学号后两位、a=1/“学号后两位”。 (1)N=10; n=0:N-1; a=1;%学号01 x=a.^n;stem(n,x);title('指数序列');grid on; (2)N=10; n=0:N-1; a=1;%学号01 x=a.^(-n);stem(n,x);title('指数序列');grid on; 实验三 离散时间信号的频域分析 一、实验目的 信号的频域分析是信号处理中一种有效的工具。在离散信号的频域分析中，通常将信号表示成单位采样序列的线性组合，而在频域中，将信号表示成复变量或的线性组合。通过这样的表示，可以将时域的离散序列映射到频域以便于进一步的处理。 在本实验中，将学习利用MATLAB计算离散时间信号的DTFT和DFT，并加深对其相互关系的理解。 二、实验原理 （１）DTFT和DFT的定义及其相互关系。 （２）使用到的MATLAB命令有基于DTFT离散时间信号分析函数以及求解序列的DFT函数。 三、实验内容及分析 （１）编程计算并画出下面DTFT的实部、虚部、幅度和相位谱。 X(e)jw0.05180.1553e11.2828ex(n)cosjwjw0.1553ej2w1.0388ej2w0.0518ej3w0.3418ej3w （２）计算32点序列 5n16，0≦n≦31的32点和64点DFT，分别绘出幅度谱图形，并绘出该序列的DTFT图形。 ３－１ clear; x=[0.0518,-0.1553,0.1553,0.0518];y=[1,1.2828,1.0388,0.3418];w=[0:500]*pi/500 H=freqz(x,y,w); magX=abs(H);angX=angle(H);realX=real(H);imagX=imag(H);subplot(221);plot(w/pi,magX);grid; xlabel('frequency in pi unit');ylabel('magnitude');title('幅度 part');axis([0 0.9 0 1.1]); subplot(223);plot(w/pi,angX);grid; xlabel('frequency in pi unit');ylabel('radians');title('相位 part');axis([0 1-3.2 3.2]); subplot(222);plot(w/pi,realX);grid; xlabel('frequency in pi unit');ylabel('real part');title('实部 part');axis([0 1-1 1]); subplot(224);plot(w/pi,imagX);grid; xlabel('frequency in pi unit');ylabel('imaginary');title('虚部 part');axis([0 1-1 1.1]); ３－２ N=32;n=0:N-1; xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('32点');figure N=64;n=0:N-1; xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('64点'); （１） （２） 实验四 离散时间LTI系统的Z域分析 一、实验目的 本实验通过使用MATLAB函数对离散时间系统的一些特性进行仿真分析，以加深对离散时间系统的零极点、稳定性，频率响应等概念的理解。学会运用MATLAB分析离散时间系统的系统函数的零极点；学会运用MATLAB分析系统函数的零极点分布与其时域特性的关系；学会运用MATLAB进行离散时间系统的频率特性分析。 二、实验原理 离散时间系统的系统函数定义为系统零状态响应的Z变化与激励的Z变化之比。 在MATLAB中系统函数的零极点可通过函数roots得到，也可借助函数tf2zp得到，tf2zp的语句格式为 [Z,P,K]=tf2zp(B,A)其中，B与A分别表示H（z）的分子与分母多项式的系数向量。它的作用是将H（z）的有理分式表示式转换为零极点增益形式。 若要获得系统函数H（z）的零极点分布图，可直接应用zplane函数，其语句格式为 Zplane（B，A） 其中，B与A分别表示H（z）的分子和分母多项式的系数向量。它的作用是在z平面上画出单位圆、零点与极点。 离散系统中z变化建立了时域函数h（n）与z域函数H（z）之间的对应关系。因此，z变化的函数H（z）从形式可以反映h（n）的部分内在性质。可根据系统的传递函数H（z）求单位冲激响应h（n）的函数impz、filter等。 利用系统的频率响应，可以分析系统对各种频率成分的响应特性，并推出系统的特性（高通、低通、带通、带阻等）。 MATLAB提供了求离散时间系统频响特性的函数freqz，调用freqz的格式主要有两种。一种形式为 [H,w]= reqz(B,A,N)其中，B与A分别表示H（z）分子和分母多项式的系数向量；N为正整数，默认值为512；返回值w包含[0，π]范围内的N个频率等分点；返回值H则是离散时间系统频率响应在0~π范围内N个频率处的值。另一种形式为 [H，w]= freqz（B，A，N，‘whole’） 与第一种方式不同之处在于角频率的范围由[0，π]扩展到[0，2π]。 三、实验内容与结果分析 已知LTI离散时间系统，要求由键盘实现系统参数输入，并绘出幅频和相频响应曲线和零极点分布图，进而分析系统的滤波特性和稳定性。 （一）程序 b=[0.0528,0.797,0.1295,0.1295,0.797,0.0528]; a=[1,-1.8107,2.4947,-1.8801,0.9537,-0.2336];w=[0:20:500]*pi/500; x1=0.0528+0.797*exp(-1*j*w)+0.1295*exp(-2*j*w)+0.1295*exp(-3*j*w)+0.797*exp(-4*j*w)+0.0528*exp(-5*j*w); x2=1-1.8107*exp(-1*j*w)+2.4947*exp(-2*j*w)+1.8801*exp(-3*j*w)+0.9537*exp(-4*j*w)+0.2336*exp(-5*j*w);x22=x2+(x2==0)*eps;x=x1./x22;magx=abs(x); angx=angle(x).*180/pi; subplot(2,2,3);zplane(b,a);title('零极点图');subplot(2,2,2);stem(w/pi,magx);title('幅度部分');ylabel('振幅');subplot(2,2,4);stem(w/pi,angx); xlabel('以pi为单位的频率');title('相位部分');ylabel('相位'); （二）波形图 图4-1 幅频、相频响应曲线、零极点分布图 实验六 IIR数字滤波器的设计 一、实验目的 从理论上讲，任何的线性是不变（LTI）离散时间系统都可以看做一个数字滤波器，因此设计数字滤波器实际就是设计离散时间系统。数字滤波器你包括IIR（无限冲激响应）和FIR（有限冲激响应）型，在设计时通常采用不同的方法。 本实验通过使用MATLAB函数对数字滤波器进行设计和和实现，要求掌握IIR数字巴特沃斯滤波器、数字切比雪夫滤波器的设计原理、设计方法和设计步骤；能根据给定的滤波器指标进行滤波器设计；同时也加深学生对数字滤波器的常用指标和设计过程的理解。 二、实验原理 在IIR滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器，然后根据一定的准则将此模拟原型滤波器转换为数字滤波器。 IIR滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定主要是在模拟原型滤波器中进行的。 IIR数字滤波器的设计方法如下：（1）冲激响应不变法。（2）双线性变化法。 一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。冲激响应不变法一个重要特点是频率坐标的变化是线性的，因此如果模拟滤波器的频率响应带限于折叠频率的话，则通过变换后滤波器的频率响应可不失真地反映原响应与频率的关系。 与冲激响应不变法比较，双线性变化的主要优点是靠频率的非线性关系得到s平面与z平面的单值一一对应关系，整个值对应于单位圆一周。所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。 MATLAB提供了一组标准的数字滤波器设计函数，大大简化了滤波器的设计工程。 （1）butter。 （2）cheby1、cheby2。 三、实验内容及分析 利用MATLAB编程方法或利用MATLAB中fdatool工具设计不同功能的IIR数字滤波器。 1、基于chebyshev I型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求参数为通带截止频率p0.4；通带最大衰减Ap1dB；阻带截止频率s0.4；阻带最小衰减As35dB。 程序: wp=0.2*pi; %通带边界频率 ws=0.4*pi; %阻带截止频率 rp=1; %通带最大衰减 rs=35; %阻带最小衰减 Fs=1000; %¼ÙÉè³éÑùÂö³å1000hz [N,Wn]=cheb1ord(wp,ws,rp,rs,'s'); [Z,P,K]=cheby1(N,rp,Wn,'s');[H,W]=zp2tf(Z,P,K); figure(1);freqs(H,W);[P,Q]=freqs(H,W);figure(2);plot(Q*Fs/(2*pi),abs(P));grid on; xlabel('频率/Hz');ylabel('幅度'); 2、基于Butterworth型模拟滤波器原型使用双线性变换方法设计数字滤波器的，要求参数为截止频率p0.4；通带最大衰减Ap1dB；阻带截止频率s0.25；阻带最小衰减AS40dB。程序: wp=0.4*pi;ws=0.25*pi;rp=1;rs=40;fs=500;ts=1/fs;wp1=wp*ts;ws1=ws*ts; wp2=2*fs*tan(wp1/2);ws2=2*fs*tan(ws1/2); [N,Wn]=buttord(wp2,ws2,rp,rs,'s');[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[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 on;xlabel('频率')ylabel('幅度')subplot(2,1,2); plot(W/pi,20*log10(abs(H)));grid on;xlabel('频率');ylabel('幅度(dB)'); 实验七 FIR数字滤波器的设计 一、实验目的 掌握用窗函数设计FIR数字滤波的原理及其设计步骤；熟悉线性相位数字滤波器的特性。学习编写数字滤波器的设计程序的方法，并能进行正确编程；根据给定的滤波器指标，给出设计步骤。 二、实验原理 如果系统的冲激响应h（n）为已知，则系统的输入输出关系为 y（n）=x（n）*h（n） 对于低通滤波器，只要设计出低通滤波器的冲激响应函数，就可以由式得到系统的输出了。 但是将h（n）作为滤波器的脉冲响应有两个问题：一是它是无限长的；二是它是非因果的。对此，采取两项措施：一是将h（n）截短；二是将其右移。 设计时，要根据阻带的最小衰减和过渡带宽度来选择恰当的窗函数类型和窗口长度N。常用的窗函数有矩形窗、海明窗和布莱克曼窗等。 窗函数设计FIR滤波器步骤如下： （1）给定理想频率响应的幅频特性和相频特性； （2）求理想单位脉冲响应，在实际计算中，可对理想频率响应采样。（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗口长度N；（4）求FIR滤波器单位脉冲响应； （5）分析幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。 三、实验内容及分析 1、分别用海明窗和布莱克曼窗设计一个48阶的FIR带通滤波器，通带为Wn0.450.55。程序1:海明窗设计 N=48; Window=hamming(N+1);w1=0.45;w2=0.55;ws=[w1,w2]; b=fir1(N,ws/pi,Window);freqz(b,1,512);title('海明窗');grid on; 程序2:莱克曼窗设计 N=48; Window=blackman(N+1);w1=0.45;w2=0.55;ws=[w1,w2]; b=fir1(N,ws/pi,Window);freqz(b,1,512);title('布莱克曼窗');grid on; 2、用矩形窗设计一个线性相位高通滤波器。其中Hejwej00.3 00.3程序: N=9; alpha=(N-1)/2;Wc=0.7*pi;n=(0:8);i=n-alpha;i=i+(i==0)*eps; h=(-1).^n.*sin((i).*Wc)./((i).*pi);%矩形窗函数设计的系统脉冲响应 w=(0:1:500)*2*pi/500; H=h*exp(-j*n'*w);%矩形窗函数设计的频响 magH=abs(H);% 矩形窗函数设计的振幅 subplot(211);stem(n,h); axis([0,8,-0.4,0.4]);title('矩形窗设计h(n)');line([0,10],[0,0]);xlabel('n');ylabel('h');subplot(212);plot(w/pi,magH); xlabel('以pi为单位的频率');ylabel('H振幅');axis([0,2,0,1.7]);title('矩形窗设计振幅谱'); 实验心得体会： 这次实验使我进一步加深了对MATLAB软件的使用。从上次的信号系统实验的初步使用到这一次的深入了解，有了更深刻的认识。对这种语言环境也有了新的了解。 在实验的过程中，我对数字滤波器的整个过程有了很好的理解和掌握。IIR数字滤波器的设计让我知道了巴特沃思滤波器和切比雪夫滤波器的频率特性，还有双线性变换及脉冲响应不变法设计的滤波器的频率特性。做这两个实验的时候程序有点困难，但经过细心的改写图形最终出来了。FIR数字滤波器的设计出来的是两种窗的图形，通过两种窗的比较，我了解了他们各自的特点，幅频和相频特性。 最后，感谢张老师对我的谆谆教导！第四篇：数字信号处理课程设计
第五篇：数字信号处理实验报告
