数字信号处理实验(武汉科技大学)

1、数字信号处理实验指导书尉 宇 编写二00七年九月目 录实验一 离散信号的MATLAB实现（1）实验二 离散傅立叶变换(DFT) （8）实验三 快速傅立叶变换(FFT) （12）实验四 设计IIR数字滤波器（16）实验五 用窗函数法设计FIR数字滤波器（20）实验一 离散信号的MATLAB实现一、 实验目的1) 学习典型的离散时间信号的MATLAB实现方法；2) 学习离散时间序列的基本运算：相加、相乘、移位等；3) 学习噪声的产生；4) 掌握两个序列的卷积和相关运算5) 掌握离散系统单位脉冲响应的求解。二、实验内容1) 典型的离散信号的表示方法：用matlab产生单位抽样信号(n) 、单位阶跃序

2、列u(n) 矩形序列RN(n) 、实指数序列 anu(n)、正弦序列sin(n)。编写程序、并画出图形。2) 离散信号的基本运算：对序列x(n)=2,3,4,1,2,5 ,n=0,1,2,3,4,5,的移位、乘法、加法、翻转及尺度变换。编写程序、并画出图形。3) 噪声的产生：产生方差为1，2，0.5的白噪声。三、实验仪器、设备及材料PC机，打印机，白纸四、 实验原理与方法五、实验步骤1)典型的离散信号的表示方法：用matlab产生单位抽样信号(n) 、单位阶跃序列u(n) 矩形序列RN(n) 、实指数序列 anu(n)、正弦序列sin(n)。编写程序、并画出图形。n=-20:20;n0=0;n

3、1=10;w0=pi/4;x=(n-n0)=0;x1=(n-n0)>=0;x2=(n-n0)>=0& (n-n1)<=0;x3=0.9.n.*x1;x4=sin(w0*n);subplot(511);stem(n,x1); title('单位抽样信号');xlabel('n');ylabel('(n)');subplot(512);stem(n,x2); title('单位阶跃序列');xlabel('n');ylabel('u(n)');subplot(513);stem

4、(n,x3); title('矩形序列');xlabel('n');ylabel('RN(n)');subplot(514);stem(n,x4); title('实指数序列');xlabel('n');ylabel('anu(n)');subplot(515);stem(n,x5); title('正弦序列');xlabel('n');ylabel('sin(wn)');xlabel('n');2)离散信号的基本运算：对序列x(n)=

5、2,3,4,1,2,5 ,n=0,1,2,3,4,5,的移位、乘法、加法、翻转及尺度变换。编写程序、并画出图形。n=-10:10;k=2;N=length(n);x=zeros(size(n);x1=zeros(size(n);y=zeros(size(n);y1=zeros(size(n);y2=zeros(size(n);y3=zeros(size(n);y4=zeros(size(n);x(11:16)=2 3 4 1 2 5;x1(11:16)=1 2 3 4 5 6;y(1+k:N)=x(1:N-k);y1(1:N-k)=x(1+k:N);y2=x.*x1;y3=x+x1;y4=fl

6、iplr(x);y5=x(1:2:N);n1=fix(n(1:2:N)/2);subplot(421);stem(n,x);xlabel('n');ylabel('x(n)');subplot(422);stem(n,x1);xlabel('n');ylabel('x1(n)');subplot(423);stem(n,y);xlabel('n');ylabel('x(n-2)');subplot(424);stem(n,y1);xlabel('n');ylabel('x(n

7、+2)');subplot(425);stem(n,y2);xlabel('n');ylabel('x(n)*x1(n)');subplot(426);stem(n,y3);xlabel('n');ylabel('x(n)+x1(n)');subplot(427);stem(n,y4);xlabel('n');ylabel('x(-n)');subplot(428);stem(n1,y5);xlabel('n');ylabel('x(2n)');3)噪声的产生

8、：产生方差为1，2，0.5的白噪声。a=0.5 1 2;N=1000;x=a'*randn(1,N);subplot(311);plot(1:N,x(1,:);subplot(312);plot(1:N,x(2,:);subplot(313);plot(1:N,x(3,:);4)用matlab求系统的单位脉冲响应 、频率响应和零极点分布。编写程序、并画出图形。 二阶离散系统的差分方程: y(n)-2.5y(n-1)+y(n-2)=f(n)-2f(n-1)Y(z)=(1-2z-1)/(1-2.5z-1+z-2)程序：b=1 -2;a=1 -2.5 1;h,t=impz(b,a,40);H

9、,W=freqz(b,a,256,'whole',1); Hr=abs(H); Hphase=angle(H);Hphase=unwrap(Hphase); subplot(221);stem(t,h);xlabel('n');ylabel('h(n)'); subplot(222);plot(W,Hr);xlabel('w');ylabel('H(ejw)'); subplot(223);plot(W,Hphase);xlabel('w');ylabel('phi(w)'); su

10、bplot(224);zplane(b,a);xlabel('Re(z)');ylabel('IM(z)');六、 实验报告要求(1) 简述实验目的及实验原理。(2) 按实验步骤附上实验过程中的信号序列、 系统单位脉冲响应及系统响应序列的时域和幅频特性曲线， 并对所得结果进行分析和解释。(3) 总结实验中的主要结论。(4) 简要回答思考题七、实验注意事项无八、思考题讨论连续正弦信号离散化为正弦序列sin(n+/6)时单位换算关系。设y(t)=sin(2ft+/6),f=20HZ,试用fs=40,80,160HZ离散化，画图，并讨论离散后的结果。f=20;fs=4

11、0 80 160;w=2*pi*f./fs; n=0:100;x=sin(w'*n+pi/6);subplot(311);stem(n,x(1,:);xlabel('n,w=1/2');ylabel('sin(w*n)');subplot(312);stem(n,x(2,:);xlabel('n,w=1/4');ylabel('sin(w*n)');subplot(313);stem(n,x(3,:);xlabel('n,w=1/8');ylabel('sin(w*n)');只有当fs>

12、;3*f时，正弦抽样信号才能恢复正弦模拟信号。实验二 离散傅立叶变换(DFT)一、 实验目的1）掌握离散傅立叶变换(DFT)及其反变换的实现。2）掌握有限长序列的循环移位、循环卷积。二、实验内容1）试用公式直接编程实现离散序列的傅立叶变换。x(n)= p(n),n=-50:50并讨论w在单位园上取10，20，100个点时，对离散序列的傅立叶变换什么影响。2）直接编程实现离散序列的傅立叶变换wc=pi/4,x(n)=sin(wc*n); x(n)=cos(wc*n);3）循环卷积：试用FFT实现x1(n)=1,2,3,4,5，x2(n)=5,4,3,2,1的循环卷积.并用图表法计算验证结果。三、

13、实验仪器、设备及材料PC机，打印机，白纸四、 实验原理与方法1）复习DFT的定义、 性质和用DFT作频谱分析的有关内容。用matlab函数dfs(dft)和idfs(idft)实现dfs(dft)正反变换。2）计算期序列的DFS.3）为了说明高密度频谱和高分辨率频谱之间的区别，考察序列x(n)有限个样本的频谱。 当时，求x(n)的DFT.当时，求x(n)的DFT.4）利用矩阵相乘计算循环卷积。设计算其5点、6点、7点、8点、9点循环卷积。并分析线性卷积与循环卷积的关系。5） 复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。6） 编写主程序。

14、 图 1 给出了主程序框图，供参考。 本实验提供FFT子程序和通用绘图子程序。2.1 主程序框图五、实验步骤 1）试用公式直接编程实现离散序列的傅立叶变换。x(n)= p(n),n=-50:50并讨论w在单位园上取10，20，100个点时，对离散序列的傅立叶变换什么影响。n=-50:50;x=ones(size(n);ws1=(-pi:pi/10:(pi-pi/10)' W1=exp(-j*ws1*n); y1=W1*x' ws2=(-pi:pi/20:(pi-pi/20)'W2=exp(-j*ws2*n); y2=W2*x'ws3=(-pi:pi/100:(p

15、i-pi/100)'W3=exp(-j*ws3*n); y3=W3*x'subplot(311);plot(ws1,abs(y1);subplot(312);plot(ws2,abs(y2);subplot(313);plot(ws3,abs(y3);2）直接编程实现离散序列的傅立叶变换wc=pi/4,x(n)=sin(wc*n); x(n)=cos(wc*n);n=-10:10;wc=pi/4;x=sin(wc*n);w2=(-4*pi:pi/5000:4*pi)'y1=exp(-j*w2*n);y2=y1*x'Hr1=imag(y2);x=cos(wc*n)

16、;y1=exp(-j*w2*n);y2=y1*x'Hr2=real(y2);subplot(121);plot(w2/(2*pi),Hr1);xlabel('w/2pi');subplot(122);plot(w2/(2*pi),Hr2);xlabel('w/2pi');3）循环卷积：试用FFT实现x1(n)=1,2,3,4,5，x2(n)=5,4,3,2,1的循环卷积.并用图表法计算验证结果。x1=1 2 3 4 5;x2=5 4 3 2 1;y1=fft(x1);y2=fft(x2);y=y1.*y2;x=ifft(y);x=abs(x)六、 实验报

17、告要求1) 简述实验原理及目的。2) 结合实验中所得给定典型序列幅频特性曲线，总结公式直接编程实现离散序列的傅立叶变换的方法。3) 总结实验所得主要结论。4) 简要回答思考题。七、实验注意事项无八、思考题讨论循环卷积和线性卷积的区别。是否可用循环卷积来计算线性卷积。若能试用循环卷积计算上述x1(n), x2(n)的线性卷积。1）循环卷积和线性卷积的区别：略2）能。将x1和x2均补零加长为N1+N2-1点长的序列。再用循环卷积来计算线性卷积：x1=1 2 3 4 5 0 0 0 0;x2=5 4 3 2 1 0 0 0 0;y1=fft(x1);y2=fft(x2);y=y1.*y2;x=iff

18、t(y);x=abs(x)实验三 快速傅立叶变换(FFT)一、目的和要求1）掌握快速傅立叶变换(FFT)的基本思想。2）了解快速傅立叶反变换(IFFT)的实现。二、实验内容1）一被噪声污染的信号，分析其频率分量。如一个由50Hz和120Hz的正弦信号构成的信号，受到高斯白噪声的干扰，数据采样率为1000Hz,数据点数N=1024。试绘制信号的频谱图和无噪声信号的频谱图。比较分析两者的区别。2) 已知一模拟信号，现在采样率f=20HZ进行采样。用DFT计算当序列长度分别为L=100,L=20时，N=200点的幅度频谱样值并通过作图进行比较。再作L=100时的反变换，与原采样信号比较，有什么区别。

19、三、实验仪器、设备及材料PC机，打印机，白纸四、 实验原理设序列x(n)的长度为N，且满足 自然数按n的奇偶把x(n)分解为两个N/2点的子序列。偶数项为一组，奇数项为一组。则x(n)的DFT为：于是，一个N点的DFT被分解为两个N/2点的DFT了，这两个N/2点的DFT再按照上式合并成一个N点DFT。五、实验步骤1）一被噪声污染的信号，分析其频率分量。如一个由Hz和Hz的正弦信号构成的信号，受到高斯白噪声的干扰，数据采样率为Hz,数据点数N=1024。试绘制信号的频谱图和无噪声信号的频谱图。比较分析两者的区别。f1=50;f2=120;fs=1000;w1=2*pi*f1/fs;w2=2*p

20、i*f2/fs;N=1000;n=0:N-1;s=sin(w1*n)+sin(w2*n);x=s+randn(size(s);subplot(221);plot(n,s);xlabel('n');ylabel('AM');axis(0 N -5 5);subplot(222);plot(n,x);xlabel('n');ylabel('AM');axis(0 N -5 5);f=linspace(0,N/2,N/2);y=fft(s);y=abs(y);y1=y(1:N/2);y=fft(x);y=abs(y);y2=y(1:N/

21、2);subplot(223);plot(f,y1);xlabel('HZ');ylabel('AM');axis( 0 N/2 min(y1) max(y1);subplot(224);plot(f,y2);xlabel('HZ');ylabel('AM');axis( 0 N/2 min(y1) max(y1);比较分析两者的区别：略2) 已知一模拟信号，现在采样率f=20HZ进行采样。用DFT计算当序列长度分别为L=100,L=20时，N=200点的幅度频谱样值并通过作图进行比较。再作L=100时的反变换，与原采样信号比较，

22、有什么区别。fs=20;Ts=1/fs;N=200;L=100;n=0:L-1;sa=exp(-n*Ts);y=fft(sa,N);y1=abs(y);y1=y1(N/2:N) y1(1:N/2-1) ;y1=y1/max(y1);L=20;n=0:L-1;sa=exp(-n*Ts);y2=fft(sa,N);y3=abs(y2);y3=y3(N/2:N) y3(1:N/2-1) ;y3=y3/max(y3);f=linspace(0,N,N);f=(f-N/2)/N;subplot(221);plot(f,y1);xlabel('L=100，N=200FFT');ylabel

23、('AM');axis(-.5 .5 0 1);subplot(222);plot(f,y3);xlabel('L=20，N=200FFT');ylabel('AM');axis(-.5 .5 0 1);L=100;n=0:L-1;x1=ifft(y,L);x1=real(x1);x2=ifft(y2,L);x2=real(x2);subplot(223);plot(n,x1);xlabel('L=100 FFT，N=100 IFFT');ylabel('AM');axis(0 L-1 0 1);subplot(2

24、24);plot(n,x2);xlabel('L=20 FFT，N=100 IFFT');ylabel('AM');axis(0 L-1 0 1);t=0:0.001:5;s=exp(-t);subplot(325);plot(t,s);xlabel('模拟信号');ylabel('AM');axis(0 5 0 1);L=100;n=0:L-1;sa=exp(-n*Ts);subplot(326);stem(n,sa);xlabel('离散信号');ylabel('AM');axis(0 L-1

25、0 1);六、 实验报告要求1) 简述实验原理及目的。2) 结合实验中所得给定典型序列幅频特性曲线， 与理论结果比较， 并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。3) 总结实验所得主要结论。4) 简要回答思考题七、实验注意事项无八、思考题将实验内容的思考部分作详细阐述。实验四 设计IIR数字滤波器一、 实验目的1) 熟悉IIR数字滤波器的设计方法。2) 掌握模拟滤波器的matlab实现。3) 熟悉用脉冲响应不变法设计IIR数字滤波器的原理与方法。4) 熟悉用双线性变换法设计IIR数字滤波器的原理与方法。5) 掌握数字滤波器的计算机仿真方法。二、实验内容1）设计一个低通巴

26、特沃斯模拟滤波器，以满足 ：通带截止频率为100HZ，阻带截止频率为200HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为1000HZ）2）设计一个高通巴特沃斯模拟滤波器，以满足 ：通带范围为6001000HZ，阻带截止频率分别为400HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为2000HZ）3）设计一个带通巴特沃斯模拟滤波器，以满足 ：通带范围为100250HZ，阻带截止频率分别为50HZ、300HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为1000HZ）4）设计一个带阻巴特沃斯模拟滤波器，以满足 ：通带截止频率为100HZ、350HZ，阻带截止频率

27、分别为150HZ、300HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为1000HZ）5）用一个低通巴特沃模拟滤波器原型设计一低通数字滤波器（脉冲响应不变、双线性变换），以满足（T=1） ：三、实验仪器、设备及材料PC机，打印机，白纸四、 实验原理复习有关巴特沃斯模拟滤波器设计、切比雪夫模拟滤波器设计、用脉冲响应不变法设计IIR数字滤波器和用双线性变换法设计IIR数字滤波器的内容，按照设计思路计算出数字滤波器系统函数H(z)。例如假设最后求出的数字滤波器系统函数H(z)为：式中A=0.090 36 B1=1.2686,C1=-0.7051 B2=1.0106,C2=-0.3583

28、B3=0.9044,C3=-0.2155 由(1式)和(2式)可见， 滤波器H(z)由三个二阶滤波器H1(z), H2(z)和H3(z)级联组成， 如图4.1 所示。图4.1 滤波器H(z)的组成五、实验步骤1）设计一个低通巴特沃斯模拟滤波器，以满足 ：通带截止频率为100HZ，阻带截止频率为200HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为1000HZ）wp=100;ws=200;rp=3;rs=30;fs=1000;wp1=wp/fs*2;ws1=ws/fs*2;N,Wn=buttord(wp1,ws1,rp,rs);b,a=butter(N,Wn);freqz(b,a,5

29、12,fs);2）设计一个高通巴特沃斯模拟滤波器，以满足 ：通带范围为6001000HZ，阻带截止频率分别为400HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为2000HZ）wp=600;ws=400;rp=3;rs=30;fs=2000;wp1=wp/fs*2;ws1=ws/fs*2;N,Wn=buttord(wp1,ws1,rp,rs);b,a=butter(N,Wn，high);freqz(b,a,512,fs);3）设计一个带通巴特沃斯模拟滤波器，以满足 ：通带范围为100250HZ，阻带截止频率分别为50HZ、300HZ，通带波汶小于3dB,阻带不大于为30dB（采样频

30、率设为1000HZ）wp=100 250;ws=50 300;rp=3;rs=30;fs=1000;wp1=wp/fs*2;ws1=ws/fs*2;N,Wn=buttord(wp1,ws1,rp,rs);b,a=butter(N,Wn);freqz(b,a,512,fs);4）设计一个带阻巴特沃斯模拟滤波器，以满足 ：通带截止频率为100HZ、350HZ，阻带截止频率分别为150HZ、300HZ，通带波汶小于3dB,阻带不大于为30dB（采样频率设为1000HZ）wp=100 350;ws=150 300;rp=3;rs=30;fs=1000;wp1=wp/fs*2;ws1=ws/fs*2;N

31、,Wn=buttord(wp1,ws1,rp,rs);b,a=butter(N,Wn,'stop');freqz(b,a,512,fs);5）用一个低通巴特沃模拟滤波器原型设计一低通数字滤波器（脉冲响应不变、双线性变换），以满足（T=1） ：wp=0.2*pi;ws=0.3*pi;rp=1;rs=15;T=1;N,Wn=buttord(wp1,ws1,rp,rs);b,a=butter(N,Wn);bz1,az1=impinvar(b,a,T)figure(1);freqz(bz1,az1,512)figure(2);bz2,az2=bilinear(b,a,T)freqz(b

32、z2,az2,512)六、 实验报告要求1) 简述实验目的及原理。2) 由所打印的|H(ej)|特性曲线及设计过程简述双线性变换法的特点。3) 对比滤波前后的信号波形， 说明数字滤波器的滤波过程与滤波作用。七、实验注意事项无八、思考题脉冲响应不变、双线性变换的各自的应用范围是什么？并简单其原理。实验五 用窗函数法设计FIR数字滤波器一、 实验目的1）学习滤波器窗函数设计法2）比较用各种窗函数设计的线性相位滤波器的频率特性。二、实验内容1）设h(n)=-4,-1,-2.5,6.5,-2,-1,1,-4,求出振幅响应Ha（w）以及H(Z)的零点位置。并判断是线性相位几型FIR滤波器。2）根据下列技

33、术指标，用窗函数法设计一个数字FIR低通滤波器，以满足 ：三、实验仪器、设备及材料PC机，打印机，白纸四、 实验原理如果所希望的滤波器的理想频率响应函数为， 则其对应的单位脉冲响应为用窗函数w(n)将hd(n)截断， 并进行加权处理， 得到:h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列， 其频率响应函数为如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性， 可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。

34、五、实验步骤1）设h(n)=-4,-1,-2.5,6.5,-2,-1,1,-4,求出振幅响应Ha（w）以及H(Z)的零点位置。并判断是线性相位几型FIR滤波器。h=-4,1,-1,-2,5,6,5,-2,-1,1,-4;M=length(h);n1=0:M-1;L=(M-1)/2;a=h(L+1) 2*h(L:-1:1);n=0:1:L;w=0:1:500'*pi/500;Hr=cos(w*n)*a'amax=max(a)+1;amin=min(a)-1;subplot(221);stem(n1,h);axis(-1 2*L+1 amin amax);xlabel('n');ylabel('h(n)');title('脉冲响应');subplot(223);stem(0:L,a);axis(-1 2*L+1 amin amax);xlabel('n');ylabel('a(n)');title('a(n)系数');subplot(222);plot(w/pi,Hr);grid;xlabel('