第一篇：数字信号处理第四章介绍
第四章 线性时不变离散时间系统的频域分析
一、传输函数和频率响应 例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确定下面这个多项式的所有根是否都在单位圆内：
由此看出，都在单位圆内。
第二篇：数字信号处理第七章介绍
第七章 数字滤波器设计
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.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].从以上结果中观察到椭圆滤波器的阶数最低，并且符合要求。
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([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。= 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。
第三篇：数字信号处理第六章介绍
第六章
数字滤波器结构
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. 课程设计报告 课程名称： 数字信号处理 课题名称： 语音信号的处理与滤波 姓 名： 学 号： 院 系： 专业班级： 指导教师： 完成日期： 2013年7月2日 目录 第1部分 课程设计报告………………………………………3 一．设计目的……………………………………………3 二．设计内容……………………………………………3 三．设计原理……………………………………………3 四．具体实现……………………………………………5 1.录制一段声音…………………………………5 2.巴特沃斯滤波器的设计………………………8 3.将声音信号送入滤波器滤波…………………13 4.语音信号的回放………………………………19 5.男女语音信号的频谱分析……………………19 6.噪声的叠加和滤除……………………………22 五． 结果分析……………………………………………27 第2部分 课程设计总结………………………………28 一． 参考文献……………………………………………28 第1部分 课程设计报告 一．设计目的 综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 二．设计内容 录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 三．设计原理 1.在Matlab软件平台下，利用函数wavrecord()，wavwrite(),wavread(),wavplay()对语音信号进行录制，存储，读取，回放。 2.用y=fft(x)对采集的信号做快速傅立叶变换，并用[h1,w]=freqz(h)进行DTFT变换。 3.掌握FIR DF线性相位的概念，即线性相位对h(n)、H()及零点的约束，了解四种FIR DF的频响特点。 4.在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。 5.抽样定理 连续信号经理想抽样后时域、频域发生的变化（理想抽样信号与连续信号频谱之间的关系） 理想抽样信号能否代表原始信号、如何不失真地还原信号即由离散信号恢复连续信号的条件（抽样定理） 理想采样过程描述： 时域描述： ˆa(t)xa(t)T(t)xa(t)(tnT)xa(nT)(tnT)xnnT(t)频域描述：利用傅氏变换的性质，时域相乘频域卷积，若 n(tnT)ˆa(t)Xa(j)xXa(j)xa(t)T(j)T(t) 则有 ˆ(j)1X(j)(j)XaaT2121ˆXa(j)Xa(jjk)Xa(jjks)TkTTkˆ(j)与X(j)的关系：理想抽样信号的频谱是连续信号频谱的Xaa 周期延拓,重复周期为s(采样角频率)。如果： X(j)Xa(j)a0s/2s/2即连续信号是带限的，且信号最高频率不超过抽样频率的二分之一，则可不失真恢复。 奈奎斯特采样定理：要使实信号采样后能够不失真还原，采样频率必须大于信号最高频率的两倍：s2h 或 fs2fh 四．具体实现 1.录制一段声音 1.1录制并分析 在MATLAB中用wavrecord、wavread、wavplay、wavwrite对声音进行录制、读取、回放、存储。 程序如下： Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号 wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound.wav');%存储语音信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(2);subplot(2,1,1);plot(x);%录制信号的时域波形 title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格 y0=fft(x);%快速傅立叶变换 figure(2);subplot(2,1,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题 grid;%加网格 图1.1 原始信号的时域与频谱图 1.2滤除无效点 针对实际发出声音落后录制动作半拍的现象，如何拔除对无效点的采样的问题: 出现这种现象的原因主要是录音开始时，人的反应慢了半拍，导致出现了一些无效点，而后而出现的无效的点，主要是已经没有声音的动作，先读取声音出来，将原始语音信号时域波形图画出来，根据己得到的信号，可以在第二次读取声音的后面设定采样点，取好有效点，画出滤除无效点后的语音信号时域波形图，对比可以看出。这样就可以解决这个问题。 x=wavread('C:UsersacerDesktop数字信号sound.wav', 7 [4000,24000]);%从4000点截取到24000结束 plot(x);%画出截取后的时域图形 title('截取后的声音时域图形');%标题 xlabel('频率');ylabel('振幅');grid;%画网格 图1.2 去除无效点 2.巴特沃斯滤波器的设计 2.1设计巴特沃思低通滤波器 MATLAB程序如下。滤波器图如图3.3所示。 %低通滤波 fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');gridon;9 图2.1 巴特沃思低通滤波器 2.2设计巴特沃思高通滤波器 MATLAB程序如下。滤波器图如图3.5所示。%高通滤波 fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);10 [N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');grid on; 图2.2巴特沃思高通滤波器 2.3设计巴特沃思带通滤波器 MATLAB程序如下。滤波器图如图3.7所示。%带通滤波 fp=[1200,3000];fs=[1000,3200];Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h));title('巴特沃斯带通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');grid on;12 图2.3巴特沃思带通滤波器 3．将声音信号送入滤波器滤波 x=wavread('C:UsersacerDesktop数字信号sound.wav');%播放原始信号 wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(4);subplot(4,2,1);plot(x);%录制信号的时域波形 title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格 y0=fft(x);%快速傅立叶变换 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题 grid;%加网格 3．1低通滤波器滤波 fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi /width);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h); ys=fftfilt(h,x);%信号送入滤波器滤波，ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(4);subplot(4,2,2);%在四行两列的第二个窗口显示图形 plot(ys);%信号的时域波形 title('低通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%网格 subplot(4,2,4);%在四行两列的第四个窗口显示图形 plot(f, abs(fftwave(n1)));%绘制模值 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('低通滤波器滤波后信号的频谱图');%标题 grid;%加网格 wavplay(ys,8000);%播放滤波后信号 3．2高通滤波器滤波 fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs; width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi/width);h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));[h1,w]=freqz(h);ys=fftfilt(h,x);%将信号送入高通滤波器滤波 subplot(4,2,5);%在四行两列的第五个窗口显示图形 plot(ys);%信号的时域波形 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('高通滤波后信号的时域波形');%标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%网格 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 subplot(4,2,7);%在四行两列的第七个窗口显示图形 plot(f,abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('高通滤波器滤波后信号的频谱图');%标题 grid;%加网格 wavplay(ys,8000);%播放滤波后信号 3．3带通滤波器 fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2;wn2=(ws2+wc2)/2;wn=[wn1 wn2];n=ceil(12.8/width*pi);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);ys1= fftfilt(h,x);%将信号送入高通滤波器滤波 figure(4);subplot(4,2,6);%在四行两列的第六个窗口显示图形 plot(ys1);%绘制后信号的时域的图形 title('带通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示纵坐标表示的意义 grid;%网格 fftwave=fft(ys1);%对滤波后的信号进行快速傅立叶变换 subplot(4,2,8);%在四行两列的第八个窗口显示图形 plot(f, abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('带通滤波器滤波后信号的频谱图');%加标题 grid;%网格 wavplay(ys1,8000);%播放滤波后信号 图形如下： 原始信号的时域波形幅值/A0-1012x 10原始信号的频谱图34幅值/A1低通滤波后信号的时域波形0.50-0.5012频率w/pi3400.51频率w/pi高通滤波后信号的时域波形幅值/A0幅值/A0幅值/Ax 10高通滤波器滤波后信号的频谱图5012频率w/pi34幅值/A0.20-0.2幅值/A2001000x 10低通滤波器滤波后信号的频谱图200100000.51频率w/pi带通滤波后信号的时域波形0.50-0.501234频率w/pix 10带通滤波器滤波后信号的频谱图50幅值 00.5频率w/pi1000.5频率w/pi1 分析：三个滤波器滤波后的声音与原来的声音都发生了变化。其中低 通的滤波后与原来声音没有很大的变化，其它两个都又明显的变化 4.语音信号的回放 sound(xlow,Fs,bits);%在Matlab中，函数sound可以对声音进行回放，其调用格式: sound(xhigh, Fs,bits);%sound(x, Fs, bits);sound(xdaitong, Fs,bits);5.男女语音信号的频谱分析 5.1 录制一段异性的声音进行频谱分析 Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号 wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound2.wav');%存储语音信号 5.2 分析男女声音的频谱 x=wavread(' C:UsersacerDesktop数字信号sound2.wav ');%播放原始信号,解决落后半拍 wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2; f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(1);subplot(2,2,1);plot(x);%录制信号的时域波形 title('原始女生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格 y0=fft(x);%快速傅立叶变换 subplot(2,2,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始女生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格 [y,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');% 对语音信号进行采样 wavplay(y,fs);%播放原始信号 N=length(y);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 subplot(2,2,3);plot(y);%录制信号的时域波形 title('原始男生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格 y0=fft(y);%快速傅立叶变换 subplot(2,2,4);%在四行两列的第三个窗口显示图形 plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始男生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格 5.3男女声音的频谱图 原始女生信号的时域波形0.50-0.5-1150100原始女生信号的频谱图幅值/A幅值 012345000x 10原始男生信号的时域波形0.50.5频率w/pi原始男生信号的频谱图1300200幅值/A0幅值 012x 1034100-0.5000.5频率w/pi1 图5.3男女声音信号波形与频谱对比 分析：就时域图看，男生的时域图中振幅比女生的高，对于频谱图女生的高频成分比较多 6.噪声的叠加和滤除 6.1录制一段背景噪声 Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号 wavwrite(x,Fs,'C:UsersacerDesktop数字信号噪音.wav');%存储语音信号 6.2 对噪声进行频谱的分析 [x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对语音信号进行采样 wavplay(x1,fs);%播放噪声信号 N=length(x1);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,1);plot(x1);%信号的时域波形 title('噪声信号的时域波形');grid;ylabel('幅值/A');y0=fft(x1);%快速傅立叶变换 subplot(3,2,2);plot(f,abs(y0(n1)));%噪声信号的频谱图 ylabel('幅值');title('噪声信号的频谱图'); 6.3原始信号与噪音的叠加 fs=8000;[x,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');%对录入信号进行采样 [x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对噪声信号进行采样 yy=x+x1;%将两个声音叠加 6.4叠加信号的频谱分析： wavplay(yy,fs);%播放叠加后信号 N=length(yy);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,3);plot(yy,'LineWidth',2);%信号的时域波形 title('叠加信号的时域波形');xlabel('时间/t');ylabel('幅值/A');grid;y0=fft(yy);%快速傅立叶变换 subplot(3,2,4);plot(f,abs(y0(n1)));%叠加信号的频谱图 title('叠加信号的频谱图');xlabel('频率w/pi');ylabel('幅值/db');grid; 6.5 设计一个合适的滤波器将噪声滤除 fs=18000;%采样频率 Wp=2*1000/fs;%通带截至频率 Ws=2*2000/fs;%阻带截至频率 Rp=1;%最大衰减 Rs=100;%最小衰减 [N,Wn]=buttord(Wp,Ws,Rp,Rs);%buttord函数（n为阶数，Wn为截至频率） [num,den]=butter(N,Wn);%butter函数（num为分子系数den为分母系数） [h,w]=freqz(num,den);%DTFT变换 ys=filter(num,den,yy);%信号送入滤波器滤波，ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(5);subplot(3,2,5);plot(ys);%信号的时域波形 title('低通滤波后信号的时域波形');%加标题 ylabel('幅值/A');%显示标表示的意义 grid;%网格 subplot(3,2,6);plot(f, abs(fftwave(n1)));%绘制模值 title('低通滤波器滤波后信号的频谱图');%标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%加网格 wavplay(ys,8000);%播放滤波后信号 grid;图形如下： 噪声信号的时域波形1100噪声信号的频谱图幅值/A0-1幅值0123450000.5叠加信号的频谱图1x 10叠加信号的时域波形10-101时间/t2200幅值/db34幅值/A100000.5频率w/pi1x 10低通滤波后信号的时域波形0.5低通滤波器滤波后信号的频谱图200幅值/A0-0.5幅值/A012x 1034100000.5频率w/pi1 图6.1噪音的叠加与滤除前后频谱对比 7.结果分析 1.录制刚开始时，常会出现实际发出声音落后录制动作半拍，可在[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav'，[100 10000])，窗的长度可根据需要定义。 2.语音信号通过低通滤波器后，把高频滤除，声音变得比较低沉。当通过高通滤波器后，把低频滤除，声音变得比较就尖锐。通过带通滤波器后，声音比较适中。 3.通过观察男生和女生图像知：时域图的振幅大小与性别无关，只与说话人音量大小有关，音量越大，振幅越大。频率图中，女生高 27 频成分较多。 4.叠加噪声后，噪声与原信号明显区分,但通过低通滤波器后，噪声没有滤除，信号产生失真。原因可能为噪声与信号频率相近无法滤除。 第2部分 课程设计总结 通过本次课程设计，使我们对数字信号处理相关知识有了更深刻的理解，尤其是对各种滤波器的设计。在设计的过程中遇到了很多问题，刚刚开始时曾天真的认为只要把以前的程序改了参数就可以用了，可是问题没有我想象中的那么简单，单纯的搬程序是不能解决问题的。通过查阅资料和请教同学收获了很多以前不懂的理论知识。再利用所学的操作，发现所写的程序还是没有能够运行，通过不断地调试，运行，最终得出了需要的结果。整个过程中学到了很多新的知识，特别是对Matlab的使用终于有些了解。在以后的学习中还需要深入了解这方面的内容。在这次的课程设计中让我体会最深的是：知识来不得半点的马虎。也认识到自己的不足，以后要进一步学习。 八．参考文献 [1]数字信号处理教程(第三版)程佩青 清华大学出版社 [2]MATLAB信号处理 刘波 文忠 电子工业出版社 [3]MATLAB7.1及其在信号处理中的应用 王宏 清华大学出版社 [4]MATLAB基础与编程入门 张威 西安电子科技大学出版社 [5] 数字信号处理及其MATLAB实验 赵红怡 张常 化学工业出版社 [6]MATLAB信号处理详解 陈亚勇等 人民邮电出版社 [7] 数字信号处理 钱同惠 机械工业出版社 29 JIANGSU UNIVERSITY OF TECHNOLOGY 数字信号处理实验报告 学院名称： 电气信息工程学院 专 业： 班 级： 姓 名： 学 号： 指导老师： 张维玺（教授） 2013年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数字滤波器的设计出来的是两种窗的图形，通过两种窗的比较，我了解了他们各自的特点，幅频和相频特性。 最后，感谢张老师对我的谆谆教导！第四篇：数字信号处理课程设计..
第五篇：数字信号处理实验报告
点击咨询 