(精品)《数字信号处理》上机实验指导

《(精品)《数字信号处理》上机实验指导》由会员分享，可在线阅读，更多相关《(精品)《数字信号处理》上机实验指导（28页珍藏版）》请在装配图网上搜索。

1、数字信号处理上机实验指导实验一、Z变换及离散时间系统分析（一）、实验目的1、通过本实验熟悉Z变换在离散时间系统分析中的地位和作用。2、掌握并熟练使用有关离散系统分析的MATLAB调用函数及格式，以深入理解离散时间系统的频率特性。（二）、实验内容及步骤 对于一个给定的LSI系统，其转移函数H(z)习惯被定义为H(z)=B(z)/A(z),即：公式中和分别是H(Z)分子与分母多项式的阶次，在有关MATLAB 的系统分析的文件中，分子和分母的系数被定义为向量，即并要求=1,如果1，则程序将自动的将其归一化为1。1、系统的阶跃响应调用格式为：y=filter(b,a,x),其中x,y,a,b都是向量。

2、例1 令求该系统的阶跃响应（y（n）。实现该任务的程序如下：clear; x=ones(100);% x(n)=1,n=1100；t=1:100;% t用于后面的绘图；b=.001836，.007344，.011016，.007374，.001836; % 形成向量b；a=1，-3.0544，3.8291，-2.2925，.55075; % 形成向量a；y=filter(b,a,x);% 求所给系统的输出，本例实际上是求所给系统的阶跃响应；plot(t,x,r.,t,y,k-);grid on;% 将x(n)（绿色）y(n)（黑色）画在同一个%图上；ylabel(x(n) and y(n)xl

3、abel(n)2、单位抽样响应h(n)调用格式为：h=impz(b，a，N) 或 h，t=impz(b，a，N)其中N是所需的h(n)的长度，前者绘图时n从1开始，而后者从0开始。例2、求上例所给系统的单位抽样响应h(n)。实现该任务的程序如下：clear;b=.001836，.007344，.011016，.007374，.001836;a=1，-3.0544，3.8291，-2.2925，.55075;h,t=impz(b，a，40); % 求单位抽样响应stem(t，h，.);grid on;3、求频率响应基本调用格式为：H，w=freqz(b，a，N,whole,Fs)其中N是频率轴的

4、分点数，建议N为2的整次幂；w是返回频率轴坐标向量，供绘图用；Fs是抽样频率，若Fs=1,频率轴给出归一化频率；whole指定计算的频率范围是从0Fs,缺省时是从0Fs/2。例3、求例1所给系统的频率响应。实现该任务的程序如下：clear all;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.8291,-2.2925,.55075;H,w=freqz(b,a,256,1);Hr=abs(H);%绝对值（幅值）；Hphase=angle(H);%相位角；Hphase=unwrap(Hphase); % 解卷绕subplot(211

5、)plot(w,Hr);grid on;ylabel( Amplitude Freq. Res.)subplot(212)plot(w,Hphase);grid on;ylabel( Phase Freq. Res.)4、离散系统的极零图调用格式：zplane(z,p) 或zplane(b,a)前者是在已知系统零点的列向量z和极点的列向量p的情况下画出的极零图，后者是在已知B(z)，A(z)的情况下的极零图。例4、显示例1系统及FIR系统的极零图。实现该任务的程序为：clear;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.82

6、91,-2.2925,.55075;subplot(221);zplane(b,a); % 求并画出所给系统的极零图，该系统为IIR系统；b=1 -1.7 1.53 -0.68;a=1;subplot(222);zplane(b,a); % 求并画出第二个系统的极零图，该系统为FIR系统（三）、作业给定系统，编程并绘出系统的单位阶跃响应y(n),频率响应。给出实验报告。clear; x=ones(100);% x(n)=1,n=1100；t=1:100;% t用于后面的绘图；b=0,0,-0.2; % 形成向量b；a=1,0,0.8; % 形成向量a；y=filter(b,a,x);% 求所给

7、系统的输出，本例实际上是求所给系统的阶跃响应；subplot(311)plot(t,x,r.,t,y,k-);grid on;% 将x(n)（绿色）y(n)（黑色）画在同一个%图上；ylabel(x(n) and y(n)xlabel(n)H,w=freqz(b,a,256,1);Hr=abs(H);%绝对值（幅值）；Hphase=angle(H);%相位角；Hphase=unwrap(Hphase); % 解卷绕subplot(312)plot(w,Hr);grid on;ylabel( Amplitude Freq. Res.)subplot(313)plot(w,Hphase);grid

8、 on;ylabel( Phase Freq. Res.)实验二、快速傅里叶变换（一）、实验目的1、通过本实验进一步加深对快速傅里叶变换的理解。2、会熟练运用fft,ifft,czt实现线性调频z变换。（二）、实验内容1、快速傅里叶变换（fft）调用格式为 X=fft(x) 或 X=fft(x,N)对前者，若x的长度是2的整次幂，则按该长度实现x的快速变换，否则，实现的是非2 的整次幂的变换；对后者，N应为2的整次幂，若x得长度小于N，则补零，若超过N，则舍弃N以后的数据。ifft的调用格式与之相同。例1、令x(n)是两个正弦信号及白噪声的叠加，试用fft文件对其作频谱分析。实现该任务的程序为

9、：clear all;% 产生两个正弦加白噪声； N=256;% 2的8次幂，进行8级蝶形运算； f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N); % 应用FFT 求频谱； subplot(3,1,1); plot(x(1:N/4);%作图向量x的第一到第N/4个值； f=-0.5:1/N:0.5-1/N; X=fft(x); % 快速傅里叶变换； y=ifft(X);% 快速傅里叶逆变化 subplot(3,1,2); plot(f,fftshift(ab

10、s(X); subplot(3,1,3); plot(real(y(1:N/4); %作图向量y的第一到第N/4个值的实部；附:fftshift函数fftshift函数就是一个交换函数： 交换规则如下： 如:x=1 2 3 4 5 6 7 8; y=fftshift(x); then y=5 6 7 8 1 2 3 4; 其在fft运算里的物理意义： 把0频(低频)周围的频谱搬移到中频范围(采样频率的一半)，只是形象化的展示FT变换后的低频成分(正负频率)。 其实质是把Fs/2的右边频谱平移到Fs/2的左边，把低频平移到Fs/2的右边，各图象间距不变。2、线性调频Z变换（CZT）CZT可用来计

11、算单位圆上任一段曲线上的Z变换，做DFT时输入的点数N和输出的点数可以不相等，从而达到频域“细化”的目的。CZT在单位圆上的Z变换就是傅里叶变换。其调用格式为： X=czt(x,M,W,A)式中x是待变换的时域信号x(n)，其长度设为N，M 是变换的长度，W确定变换的步长，A确定变换的起点。若M=N,A=1，则CZT变成DFT。CZT应用举例：例2、设x(n)由三个实正弦组成，频率分别是8Hz,8.22Hz,9Hz,抽样频率为40 Hz,时域取128点，作CZT变换和FFT变换。实现该任务的程序为：clear all；% 构造三个不同频率的正弦信号的叠加作为试验信号N=128;f1=8;f2=

12、8.22;f3=9;fs=40;stepf=fs/N;n=0:N-1;t=2*pi*n/fs;n1=0:stepf:fs/2-stepf;x=sin(f1*t)+sin(f2*t)+sin(f3*t);M=N;W=exp(-j*2*pi/M);% A=1时的czt变换A=1;Y1=czt(x,M,W,A);subplot(311)plot(n1,abs(Y1(1:N/2);grid on;% DTFTY2=abs(fft(x);subplot(312)plot(n1,abs(Y2(1:N/2);grid on;% 详细构造A后的czt M=60;f0=7.2;DELf=0.05;A=exp(j

13、*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);Y3=czt(x,M,W,A);n2=f0:DELf:f0+(M-1)*DELf;subplot(313);plot(n2,abs(Y3);grid on;（三）、作业设x(n)由三个实正弦组成，频率分别是8Hz,9Hz,10Hz,抽样频率为60 Hz,时域取256点，作CZT变换、IFFT变换和FFT变换，观察波形，更改参数，得出不同参数下的CZT变换波形。给出实验报告。clear all;N=256;% 2的8次幂，进行8级蝶形运算f1=3;f2=5;f3=7;fs=40;a1=1;a2=3;a3=5;stepf=fs

14、/N;n=0:N-1;t=2*pi*n/fs;n1=0:stepf:fs/2-stepf;x=a1*sin(f1*t)+a2*sin(f2*t)+a3*sin(f3*t);M=N;W=exp(-j*2*pi/M);% 三个不同频率的正弦信号； subplot(511); plot(n,x);grid on;ylabel(one) % 应用FFT 求频谱； X=fft(x); % 快速傅里叶变换； subplot(512); plot(n1,abs(X(1:N/2);grid on; ylabel(two) y=ifft(X);% 快速傅里叶逆变化 subplot(513); plot(real

15、(y(1:N);grid on;ylabel(three)% A=1时的czt变换A=1;Y1=czt(x,M,W,A);subplot(514)plot(n1,abs(Y1(1:N/2);grid on;ylabel(four)% 详细构造A后的czt M=60;f0=7.2;DELf=0.05;A=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);Y3=czt(x,M,W,A);n2=f0:DELf:f0+(M-1)*DELf;subplot(515);plot(n2,abs(Y3);grid on;ylabel(five)实验三、无限冲击响应数字滤波器设计

16、（一） 实验目的1、 要求掌握IIR数字滤波器的设计原理、设计方法和设计步骤；2、 能根据给定的滤波器指标进行滤波器设计；3、 掌握数字巴特沃斯滤波器、数字切比雪夫滤波器的设计原理和步骤；（二）、实验内容IIR数字滤波器的设计有多种方法，如频率变换法、数字域直接设计以及计算辅助设计等。下面只介绍频率变换设计法。首先考虑由模拟低通滤波器到数字低通滤波器的转换，其基本的设计过程如下:1、 将数字滤波器的技术指标转换为模拟滤波器的技术指标；2、 设计模拟滤波器G(S)；3、 将G(S)转换成数字滤波器H(Z)在低通滤波器的设计基础上，可以得到数字高通、带通、带阻滤波器的设计流程，以高通数字滤波器的设

17、计为例：1、 将数字高通滤波器的技术指标，通过转变为模拟高通的技术指标，作归一化处理后得；2、 利用频率变换关系，将模拟高通的技术指标转换为归一化的低通滤波器G(p)的技术指标，并有p=;3、 设计模拟低通滤波器G(p);4、 将G(p)转换为模拟高通滤波器的转移函数,p=;5、 将转换成数字高通滤波器的转移函数，s=(z-1)/(z+1)。以上5个步骤同样适用于数字带通、数字带阻滤波器的设计。只是在步骤2，3，4中频率转换的方法不同。例1、试设计一个数字高通滤波器，要求带通下限频率=0.8pi,阻带上限频率为0.44pi,通带衰减不大于3dB,阻带衰减不小于20 dB。实现该任务的MATLA

18、B程序为：clear allwp=.8*pi;ws=.44*pi;rp=3;rs=20;% 通带和阻带的衰减Fs=2000;% 抽样衰减%数字高通滤波器的技术指标% Firstly to finish frequency prewarping;wap=2*Fs*tan(wp/2)was=2*Fs*tan(ws/2);% 模拟高通滤波器的技术指标n,wn=buttord(wap,was,rp,rs,s);% 确定模拟低通滤波器的阶数% Note: s!z,p,k=buttap(n);% 设计模拟低通原型滤波器，z,p,k分别是设计出低通原型滤波器的极点，零点及增益b,a=zp2tf(z,p,k)

19、% 零、极点系数转换为传递函数% 巴特沃斯模拟低通原型滤波器频率响应bt,at=lp2hp(b,a,wap)%由低通滤波器转移函数产生模拟高通滤波器转移函数% Note: z=(2/ts)(z-1)/(z+1);2/ts=1,that is 2fs=1,fs=0.5;bz,az=bilinear(bt,at,Fs)% 由模拟高通滤波器H(s)得到数字高通滤波器H(z)h,w=freqz(bz,az,256,1);%得到频率响应figure(1)plot(w,20*log10(abs(h);grid;%得出幅频图ylabel( High-pass DF)% Directly to design

20、H(z) by butter.mn,wn=buttord(wp/pi,ws/pi,rp,rs);% 确定模拟低通滤波器的阶数b,a=butter(n,wp/pi,high);%直接设计巴特沃思滤波器h1,w1=freqz(b,a,256,1);h1=20*log10(abs(h1);figure(2)plot(w,h1);gridylabel( High-pass DF)例2、一个数字系统的抽样频率为1000Hz,已知该系统受到频率为100Hz的噪声的干扰，现设计一陷波滤波器去掉该噪声。要求3dB的带边频率为95Hz和105Hz,阻带衰减不小于14Db.阻带的下边和上边频率分别为99Hz和10

21、1Hz。实现程序为：fp=95 105;fs=99 101;%wp=.19*pi 0.21*pi;ws=.198*pi 0.202*pi;Fs=1000;rp=3;rs=14;wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;% Firstly to finish frequency prewarping;wap=2*Fs*tan(wp./2)was=2*Fs*tan(ws./2);n,wn=buttord(wap,was,rp,rs,s);% Note: s!z,p,k=buttap(n);b,a=zp2tf(z,p,k)bw=wap(2)-wap(1)w0=sqrt(wap(1)*w

22、ap(2)bt,at=lp2bs(b,a,w0,bw)%将模拟低通滤波器G(p)转换为实际带阻滤波器% Note: z=(2/ts)(z-1)/(z+1);bz1,az1=bilinear(bt,at,Fs)h,w=freqz(bz1,az1,256,Fs);figure(1)plot(w,20*log10(abs(h)% Directly to design H(z) by butter.mn,wn=buttord(wp/pi,ws/pi,rp,rs);b,a=butter(n,wp/pi,stop)h1,w1=freqz(b,a,256,Fs);figure(2)plot(w1,20*lo

23、g10(abs(h1)附：buttord.m本文件用来确定数字低通或模拟低通滤波器的阶数，其调用格式为(1) N,Wn= buttord(Wp, Ws, Rp, Rs)(2) N,Wn= buttord(Wp, Ws, Rp, Rs,s)格式（1）对应数字滤波器，式中Wp, Ws分别是通带和阻带的截止频率，实际上他们是归一化频率，其值在01之间，1对应抽样频率的一半。Rp, Rs分别是通带和阻带的衰减，单位是Db。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率。格式（2）对应模拟滤波器，式中各个变量的含义与格式1相同，但Wp, Ws及Wn的单位为rad/s,因此它们实际上是频率。but

24、ter.m本文件可用来直接设计巴特沃思数字滤波器。其调用格式为（1）B,A=butter（N，Wn）;（2）B,A=butter（N，Wn，high）;（3）B,A=butter（N，Wn，stop）;（4）B,A=butter（N，Wn，s）;格式13用来设计数字滤波器,所以B,A分别是H(z)的分子，分母的系数向量，Wn是通带截止频率，其值在01之间，1对应抽样频率的一半。若Wn是标量，则格式1用来设计低通数字滤波器，若Wn是1*2的向量，则格式1用来设计带通数字滤波器，2用来设计高通数字滤波器，3用来设计带阻数字滤波器，4用来设计模拟滤波器。lp2bp.m本文件用来将模拟低通滤波器G(p

25、)转换为实际的带通滤波器，其调用格式为 B,A= lp2bp(b,a,WO, BW)式中b,a分别是模拟低通原型滤波器G(p)的分子、分母多项式的系数向量，B,A分别是转换后的H(s)的分子、分母多项式的系数向量，WO 是带通滤波器的中心频率，BW是其带宽。（三）、作业设计一个数字带通滤波器，参数自定。给出实验报告。clear all;fp=100 400;fs=50 450;Fs=1000;% 抽样衰减 rp=3;rs=18;% 通带和阻带的衰减 %数字高通滤波器的技术指标wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;wap=2*Fs*tan(wp./2)was=2*Fs*tan(

26、ws./2);n,wn=buttord(wap,was,rp,rs,s);z,p,k=buttap(n);% 设计模拟低通原型滤波器，z,p,k分别是设计出低通原型滤波器的极点，零点及增益b,a=zp2tf(z,p,k)% 零、极点系数转换为传递函数% 巴特沃斯模拟低通原型滤波器频率响应bw=wap(2)-wap(1)w0=sqrt(wap(1)*wap(2)bt,at=lp2bp(b,a,w0,bw)bz1,az1=bilinear(bt,at,Fs)h,w=freqz(bz1,az1,256,Fs);%得到频率响应figure(1)plot(w,20*log10(abs(h)grid on

27、实验四、有限冲击响应数字滤波器设计(一)、实验目的1、了解无限冲击响应数字滤波器设计和有限冲击响应数字滤波器设计各自的特点，比较两者的优缺点。2、掌握用窗函数法设计FIR 数字滤波器的原理和方法。3、熟悉线性相位FIR 数字滤波器特性。4、了解各种窗函数对滤波特性的影响。（二）、实验原理及内容窗口法的基本思想：根据给定的滤波器技术指标，选择滤波器长度N和窗函数，使其具有最窄宽度的主瓣和最小的旁瓣。1、窗口法的设计步骤：（1）、给定理想频响函数；（2）、根据指标选择窗函数。确定窗函数类型的主要依据是过度带宽和阻带最小衰耗的指标；（3）、由求，加窗得；（4）、检验。由求，求是否在误差容限之内。2、

28、常用MATLAB文件介绍fir1.m本文件采用窗函数法设计FIR数字滤波器，其调用格式为(1) b=fir1(N,Wn)(2) b=fir1(N,Wn,high)(3) b=fir1(N,Wn,stop)式中N为滤波器的阶数，因此滤波器的长度为N+1; Wn是通带截止频率，其值在01之间，1对应抽样频率的一半；b是设计好的滤波器系数h(n).对格式1，若Wn是一标量，则可用来设计低通滤波器；若Wn是1*2的向量，则用来设计带通滤波器；若Wn是1*L的向量，则可用来设计L带通滤波器，这时，格式1 要改为b=fir1(N,Wn,DC-1) 或 b=fir1(N,Wn,DC-0)前者保证第一个带为通

29、带，后者保证第一个带为阻带。格式2用来设计高通滤波器，3用来设计带阻滤波器。在上述所有格式中，若不指定窗函数的类型，fir1自动选择汉明窗。例1、 设计一FIR低通滤波器，所希望的频率响应在00.25pi之间，在0.25pipi之间为0，令N=10，分别用矩形窗和汉明窗设计，观察其幅频响应的特性。实现该任务的程序为：clear all;N=10;b1=fir1(N,0.25,boxcar(N+1); % 用矩形窗作为冲激响应的窗函数b2=fir1(N,0.25,hamming(N+1);% 用Hamming窗作为冲激响应的窗函数M=128;h1=freqz(b1,1,M);h2=freqz(b

30、2,1,M);% 分别求两个滤波器的频率响应；t=0:10;subplot(221)stem(t,b2,.);hold on;% 绘制火柴梗图；plot(t,zeros(1,11);grid;%绘制1*11的零数组网格图；f=0:0.5/M:0.5-0.5/M;subplot(222)plot(f,abs(h1),b-,f,abs(h2),g-);grid;fir2.m本文件采用窗函数法设计具有任意幅频响应的FIR数字滤波器，其调用格式为b=fir2(N,F,M)其中F是频率向量，其值在01之间，M是和F相对应的所希望的幅频响应。如同fir1缺省时自动选用汉明窗。例2、 设计一多带滤波器，要求

31、理想幅频响应在归一化频率 0.20.3，0.60.8之间为1，其余处为0。实现程序如下（给出了两个不同长度滤波器以作比较）：clear all;f=0 0.19 0.2 0.3 0.31 0.59 0.6 0.8 0.81 1; % 给定频率轴分点；m=0 0 1 1 0 0 1 1 0 0;% 给定在这些频率分点上理想的幅频响应N1=30;N2=90;% 取两种不同的滤波器长度；b1=fir2(N1,f,m);b2=fir2(N2,f,m);% 得到两个滤波器；subplot(311);stem(b1,.);grid;subplot(312);stem(b2,.);grid;M=128;h1

32、,w=freqz(b1,1,M,1);h2,w=freqz(b2,1,M,1);subplot(313);plot(w,abs(h1),b-,w,abs(h2),g-);grid;（三）、作业设计一FIR低通滤波器，所希望的频率响应在00.3pi之间，在0.3pipi之间为0，分别取N=10，20，40，自行选择窗函数，观察其幅频响应的特性。给出实验报告。clear all;N=10;b1=fir1(N,0.3,boxcar(N+1);b2=fir1(N,0.3,hamming(N+1);M=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);t=0:10;subplo

33、t(221)stem(t,b2,.);hold on;plot(t,zeros(1,11);grid;f=0:0.5/M:0.5-0.5/M;subplot(222)plot(f,abs(h1),b-,f,abs(h2),g-);grid; clear all;N=20;b1=fir1(N,0.3,boxcar(N+1);b2=fir1(N,0.3,hamming(N+1);M=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);t=0:20;subplot(221)stem(t,b2,.);hold on;plot(t,zeros(1,21);grid;f=0:0.5

34、/M:0.5-0.5/M;subplot(222)plot(f,abs(h1),b-,f,abs(h2),g-);grid; clear all;N=40;b1=fir1(N,0.3,boxcar(N+1);b2=fir1(N,0.3,hamming(N+1);M=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);t=0:40;subplot(221)stem(t,b2,.);hold on;plot(t,zeros(1,41);grid;f=0:0.5/M:0.5-0.5/M;subplot(222)plot(f,abs(h1),b-,f,abs(h2),g-);grid; 28