数字信号处理实验指导书(广东石油化工学院).docx

《数字信号处理实验指导书(广东石油化工学院).docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书(广东石油化工学院).docx（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 数字信号处理实验指导书广东石油化工学院电子信息工程系 实验 1 信号产生的基本方法.1实验 2 离散时间系统的时域分析.6实验 3 FFT 算法与离散时间信号的频域分析.1实验 4 离散系统的变换域分析.9实验 5 无限冲激响应数字滤波器的设计.5实验 6 有限冲激响应数字滤波器设计.1 实验 1 信号产生的基本方法一、实验目的1、加深对常用离散信号的理解；2、熟悉使用 MATLAB信号处理的基本概念。二、实验原理1、信号或序列的表示方法信号在 Matlab中用数据阵列来表示，数据阵列可以是实数或复数，也可以是一维或多维。一维抽样信号或序列用 1n或 n1的的矢量（阵列）来表示，例如：x=2

2、,3,1,4,-9，若矢量x转置，x=x，则x=2314-9多维信号用多维阵列来表示，例如，三通道信号，可以是这样：y=x 2*x x/piy =2.00003.00001.00004.0000-9.00004.00006.00000.63660.95490.31831.2732-2.86482.00008.0000-18.00002、离散时间信号的表示方法在 Matlab中，大多数工具箱函数需要用一个矢量表示时间。例如，考虑产生 1000Hz的抽样频率的数据，时间可以是t = (0:0.001:1);或 t = 0:0.001:1。时间 t是 1001个元素的矢量，时间是按间隔 1ms从 0

3、到 1抽样。给定时间 t，抽样频率 1000Hz，我们可以产生一个 70Hz和 130Hz的持续时间 1s的正弦信号的和：y = 3*sin(2*pi*70*t) + 4*sin(2*pi*130*t);离散时间信号 y有 1001个元素。3、单位抽样序列、单位阶跃序列、斜坡序列的表示方法 1 n = 0(n) =d单位取样序列：0 n 0，在 Matlab中实现过程如下：imp= 1; zeros(99,1);% 单位取样信号1 n 00 n 0单位阶跃序列： u(n) =，在 Matlab中实现过程如下：unit_step = ones(100,1); % 从 0开始单位阶跃信号t t 0

4、斜坡信号： ( ) =r t，在 Matlab中实现过程如下：0 t 0t = (0:0.001:1);ramp_sig= t;% 斜坡信号4、复指数序列、指数序列的表示方法pj n3(n) = re ，例如，序列 x(n) = 3ewj nx复制数序列：(n = 0,1,.,9) ，在 MATLAB 中实现过程如下：w=pi/3;n=0:9;x=r*exp(j*w*n);x(n) = a u(n)，例如， x(n) = 0.7 u(n) ，指数序列：在 MATLAB中实现过程如下：nnN=10;a=0.7;n=0:N-1;x=a.n;5、模拟时间信号的绘图模拟时间信号的绘图采用 plot命令

5、。例如，产生一个 时间为 1.5 s ，频率 为 50 Hz 的锯齿波，抽样频率为 1kHz，画出 0.2 s 的波形。fs = 1000;t = 0:1/fs:1.5;x = sawtooth(2*pi*50*t);plot(t,x), axis(0 0.2 -1 1)程序运行结果如下：10.50-0.5-100.050.10.150.2 又例如：例如，产生一个 时间为 1.5 s ，信号为 y(t) = 3sin(140pt) + 4*sin(260pt)，抽样频率为 10kHz，画出 0.2 s 的波形。在 MATLAB中实现过程如下：fs = 10000;t = 0:1/fs:1.5;

6、x = 3*sin(2*pi*70*t) + 4*sin(2*pi*130*t);plot(t,x), axis(0 0.2 -8 8)程序运行结果如下：86420-2-4-6-800.050.10.150.2sin(pt)再例如：产生 sinc信号 x(t) =如下：，绘制从 t=-5到 t=5的波形图,在 MATLAB中实现过程ptx = linspace(-5,5);y = sinc(x);plot(x,y)程序运行结果如下：10.50-0.5-5056、离散时间信号的绘图离散序列的绘图采用 stem命令。例如，绘制 x(n) = e-0.07n cos(n) 的波形，在 MATLAB中

7、实现过程如下：n = 0:25; x = exp(-.07*n).*cos(n);stem(n,x)程序运行结果如下：10.50-0.5-10510152025sin(pn / 4)又例如：产生 sinc序列 x(n) =程如下：，绘制从 n=-5到 n=5的波形图,在 MATLAB中实现过p / 4nn= -25:25;x = sinc(n/4);stem(n,x,.)程序运行结果如下：10.50-0.5-30-20-100102030三、预习要求1、预先阅读实验原理（MATLAB基础介绍）；2、讨论正弦序列、复指数序列的性质。1p1p1ppA绘出离散信号 x(n) = e ，当 = -zn

8、+、 =+、 =、 = 1.7 +j和 =z jzjzjzz1271271277时的信号实部和虚部图？周期信号的周期为多少？B绘出信号 x(n) = 0.7sin(2 p 0.3n) 的数字频率是多少？周期是多少？产生一个数字频率为0.7的正弦序列，并显示该信号，说明其周期。 3、使用 Matlab 帮助功能学习 square(方波),sawtooth(锯齿波)和 sinc 函数，并绘图。四、实验内容1、编制程序产生复指数信号、指数衰减信号、square(方波)、sawtooth(锯齿波)和 sinc函数的模拟信号，时间和抽样频率可输入确定，函数需要的参数也可输入确定，并绘出其图形。2、编制程

9、序产生复指数序列、指数衰减序列、square(方波)、sawtooth(锯齿波)和 sinc函数的离散信号，长度可输入确定，函数需要的参数可输入确定，并绘出其图形。 实验 2 离散时间系统的时域分析一、实验目的1、熟悉并掌握离散时间系统的差分方程表示法；2、加深对冲激响应和卷积分析方法的理解。二、实验原理在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:x(n)离散时间系统其输入、输出关系可用以下差分方程描述：MNy(n) = b x(n - k) - a y(n - k)kkk=0k=1a = 0,k = 1,2,., Nk当时，hn是有限长

10、度的，称系统为 FIR 系统；反之，称系统为 IIR系统。1 n = 0(n) =d当输入信号分单位取样序列：0 n 0，系统输出为单位取样响应 h(n) 。若深入信号 x(n) ，则系统响应为如下的卷积计算式：y(n) = x(n)* h(n) =x(k)h(n - k)k=-三、预习要求 y = filter(A,B, x)1、在 MATLAB中，熟悉利用函数实现差分方程的仿真；y = cov(h, x)y = impz(A,B, N)求系统冲激响应2、在 MATLAB中，熟悉用函数计算卷积，用的过程。四、实验内容1、分别使用 conv和 filter函数计算序列 h(n)和 x(n)的卷

11、积，运行下面程序，并分析两个结果是否有差别，在使用 filter完成卷积功能时，需要如何补零？例程如下：clf;h = 3 2 1 -2 1 0 -4 0 3; %脉冲响应x = 1 -2 3 -4 3 2 1; %输入序列y = conv(h,x);n = 0:14;subplot(2,1,1);stem(n,y);xlabel( n); ylabel(幅度);title(卷积); grid;x1 = x zeros(1,8);y1 = filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel(n); ylabel(幅度);title(filter的结果

12、); grid;2、编制程序求解下列两个系统的单位取样响应，并绘出其图形。要求分别用 filter、conv、impz三种函数完成。y(n) + 0.67y(n -1)+ 0.126y(n - 2) = x(n) - x(n -1)1y(n) = x(n) + x(n -1)+ x(n - 2) + x(n - 3) + x(n - 4)，5给出理论计算结果和程序计算结果并讨论。 实验 3 离散系统的变换域分析一、实验目的1、熟悉对离散系统的频率响应分析方法；2、加深对零、极点分布的概念理解。二、实验原理离散系统的时域方程为：MNy(n) = b x(n - k) - a y(n - k)kk

13、k=0k=1其变换域分析方法如下：y(n) = x(n)* h(n) Y (e w ) = X (e w )H(e w )频域:Z 域:。jjjy(n) = x(n)* h(n) Y (z) = X (z)H(z)，系统的转移函数为:B(z) b + b z + + b z-1-NH(z) =01NA(z) a + a z + + a z-1-M01MB(z)A(z)的根为系统的极点。其中的根称为系统的零点、三、预习要求1. 在 MATLAB 中，熟悉函数 tf2zp、zplane、freqz、residuez、zp2sos 的使用，其中：z，p，K=tf2zp（b，a）求得有理分式形式的系统

14、转移函数的零、极点；zplane（z，p）绘制零、极点分布图；h=freqz(b,a,w)求系统的单位频率响应；r，p，k=residuez（b，a）完成部分分式展开计算；sos=zp2sos（z，p，K）完成将高阶系统分解为 2 阶系统的串联。2. 阅读实验内容的参考程序，学习频率分析法在 MATLAB 中的实现；3. 编程实现系统参数输入，绘出幅度频率响应和相位响应曲线和零、极点分布图。 四、实验内容0.0528 + 0.0797z + 0.1295z + 0.1295z + 0.797z + 0.0528z-1-2-3-4-51、求系统 H(z) =1-1.8107z + 2.4947z

15、 -1.8801z + 0.9537z - 0.2336z-1-2-3-4-5的零、极点和幅度频率响应和相位响应。2、求下列直接型系统函数的零、极点，并将它转换成二阶节形式1- 0.1z - 0.3z - 0.3z - 0.2z-1-2-3-4-4H(z) =1+ 0.1z + 0.2z + 0.2z + 0.5z-1-2-3用 MATLAB参考程序如下：b=1 -0.1 -0.3 -0.3 -0.2;a=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(b,a);m=abs(p);disp(零点);disp(z);disp(极点);disp(p);disp(增益系数);disp(k

16、);sos=zp2sos(z,p,k);disp(二阶节);disp(real(sos);zplane(b,a)3. 差分方程：y(n) + 0.7y(n -1) - 0.45y(n - 2) - 0.6y(n - 3)= 0.8x(n) - 0.44x(n -1) + 0.36x(n - 2) + 0.02x(n - 3)所对应的系统的频率响应。MATLAB计算的程序参考如下：k=256;b=0.8 -0.44 0.36 0.02;a=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(b,a,w);subplot(2,2,1);plot(w/pi,real(h);gr

17、idtitle(实部)xlabel(omega/pi);ylabel(幅度)subplot(2,2,2);plot(w/pi,imag(h);grid title(虚部)xlabel(omega/pi);ylabel(Amplitude)subplot(2,2,3);plot(w/pi,abs(h);gridtitle(幅度谱)xlabel(omega/pi);ylabel(幅值)subplot(2,2,4);plot(w/pi,angle(h);gridtitle(相位谱)xlabel(omega/pi);ylabel(弧度) 实验 4 FFT 算法与离散时间信号的频域分析一、实验目的1、加

18、深对离散信号的 DFT 的理解；2、熟练使用 FFT 函数并对信号进行幅频和相频分析。二、实验原理在 MATLAB 中，可以用函数 X=fft(x,N)和 x=ifft(X,N)计算 N 点序列的 DFT 正、反变换。N点序列的 DFT 和 IDFT 变换定义式如下：N -1X (k) = x(n)W,nkNn=01 NN -1x(n) =X (k)W-nk ,Nn=0p- j2W = e利用旋转因子具有周期性，可以得到快速算法（FFT）。nkNNX(k) = X(k) eX(k)arg(X (k)，相位为 。在频域中，jarg(X (k) ，幅度为三、预习要求1、在 MATLAB 中，熟悉函

19、数 fft、ifft 的使用；2、阅读扩展练习中的实例，学习在 MATLAB 中的实现 FFT 算法的实现；3、利用 MATLAB 编程完成计算，绘出相应图形。四、实验内容2p12px(n) = cos( 7n) + cos( 7n), n = 0,1, , N -11、实序列。分别用一个 64 点和 128N2NX(k)arg(X (k)，并绘出的图形。点的 fft 函数计算幅度为，相位为 2、已知某序列 x(n)在单位圆上的 N=64点傅立叶变换为:1X (k) =,k = 0,1,2, , N -11- 0.7ep- j2 k/Nx(n) = IDFTX (k)用 N点 ifft程序计算

20、。3. 在噪声信号中寻找频率信号，考虑抽样频率 1000Hz，信号为频率 50 Hz幅度 0.7和频率 120Hz幅度 1的两个正弦波之和，它被零均值的随机噪声污染，程序如下：Fs = 1000;% 抽样频率% 抽样间隔% 信号的长度% 时间矢量T = 1/Fs;L = 1000;t = (0:L-1)*T;% 50 Hz 和 120 Hz正弦波相加x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);y = x + 2*randn(size(t);plot(Fs*t(1:50),y(1:50)title(被噪声污染的信号)xlabel(时间 (ms)% 加噪声被噪

21、声污染的信号6420-2-4-601020304050时间 (ms)NFFT = 500;Y = fft(y,NFFT)/L;f = Fs/2*linspace(0,1,NFFT/2+1);% 绘制单边幅度谱plot(f,2*abs(Y(1:NFFT/2+1)title(信号 y(t)的单边幅度谱)xlabel(频率(Hz)ylabel(|Y(f)|) 信号y(t) 的单边幅度谱0.50.40.30.20.10|)f(Y|0100200300400500频率(Hz)f = 8 fs4. 对连续的单一频率周期信号按采样频率=16，观察其 DFT结果的幅度谱。此时离散序列采样，截取长度 分别选 N

22、 =20和NNax(n) = sin(2 nf / f ) = sin(2pn / 8)p，即 k=8。as用 MATLAB计算并作图，函数 fft用于计算离散傅里叶变换 DFT，程序如下：k=8;n1=0:1:19;xa1=sin(2*pi*n1/k);subplot(2,2,1)plot(n1,xa1)xlabel(t/T);ylabel(x(n);xk1=fft(xa1);xk1=abs(xk1);subplot(2,2,2)stem(n1,xk1)xlabel(k);ylabel(X(k);n2=0:1:15;xa2=sin(2*pi*n2/k);subplot(2,2,3)plot(

23、n2,xa2)xlabel(t/T);ylabel(x(n);xk2=fft(xa2);xk2=abs(xk2);subplot(2,2,4)stem(n2,xk2)xlabel(k);ylabel(X(k);程序运行结果： 10.5086)Xnk 4(x-0.5-120010t/T20010k2010.501086)Xnk(x420-0.5-1051015051015t/Tk计算结果示于图 1，(a)和(b)分别是 N =20 时的截取信号和 DFT 结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是 N=16 时的截取信号和 DFT 结果，由于截取了两个整周期，得到单一谱线

24、的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。 实验 5 无限冲激响应数字滤波器的设计一、实验目的1、掌握 IIR数字滤波器的具体设计方法;2、熟悉设计低通、高通、带通和带阻 IIR数字滤波器的几种方法。二、实验原理1、数字滤波器设计的参数指标：w，阻带截止频率w ，通带增大衰减a设计数字滤波器的技术指标：通带截止频率和阻psp带的最大a。sw和阻带截止频率w为二元向量。在带通滤波器的设计中，通带截止频率ps2、滤波器的设计中，数字频率与模拟频率的关系：数字频率 ，模拟频率为 f或 ，抽样频率 fs，归一化频率 f=f/fs或 = / s， = /2 ， 因为 f、

25、 和 ， 都是无量纲， 刻度是一样的， 将 f、 、 、 f、 、 的定标值对应关系。 fs 1fs- f 2f 20fss0.500.51f -WWssss 10.50.51Ww-2-02w 10.500.513、 数字滤波器 H(z)B(z) b + b z + + b z-1-NH(z) =01NA(z) a + a z + + a z-1-M01M三、预习要求 1. 熟悉MATLAB中巴特沃斯滤波器的设计函数：N,wc=buttord(wp,ws,Rp,As); N,wc=buttord(wp,ws,Rp,As,s); B,A=butter(N,wc,ftype)；B,A=butter

26、(N,wc,ftype,s)。2. 熟悉MATLAB中切比雪夫型滤波器的设计函数：N,wpo= cheb1ord(wp,ws,Rp,As); N,wpo= cheb1ord(wp,ws,Rp,As,s); B,A=cheby1(N,Rp, wpo,ftype); B,A= cheby1(N,Rp, wpo,ftype,s) 。3. 熟悉MATLAB中切比雪夫型滤波器的设计函数：N,wpo= cheb2ord(wp,ws,Rp,As); N,wpo= cheb2ord(wp,ws,Rp,As,s); B,A=cheby2(N,Rp, wpo,ftype); B,A= cheby2(N,Rp, w

27、po,ftype,s) 。4. 熟悉MATLAB中椭圆滤波器的设计函数：N,wpo= cheb2ord(wp,ws,Rp,As); N,wpo= cheb2ord(wp,ws,Rp,As,s); B,A=cheby2(N,Rp, wpo,ftype); B,A= cheby2(N,Rp, wpo,ftype,s) 。四、实验内容1. 利用设计一个椭圆数字带通滤波器，指标要求如下：w 0.45p w0.65p1= dB，通带峰值起伏：a；通带边缘频率：阻带边缘频率：，p1p2pw 0.3p w=0.8pa 40 dB。=，最小阻带衰减：s1s2s2. 设采样周期T=250s（采样频率fs =4k

28、Hz），设计一个三阶巴特沃斯波器，其3dB边界频率为fc =1kHz。3、设计一切比雪夫型数字高通滤波器，它的通带为400500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。4、 设计切比雪夫型带通滤波器，其dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3 =80kHz和f4=120kHz处的最小衰减大于dB，采样频率fs=400kHz。5、一数字滤波器采样频率fs=1kHz，要求滤除100Hz的干扰，其dB的边界频率为95Hz和105Hz，原型归一化低通滤波器的参考程序为：w1=95/500;w2=105/5

29、00;B,A=butter(1,w1, w2,stop);h,w=freqz(B,A);f=w/pi*500;plot(f,20*log10(abs(h); axis(50,150,-30,10);grid;xlabel(频率/Hz)ylabel(幅度/dB) 实验 6 有限冲激响应数字滤波器设计一、实验目的：1、加深对数字滤波器的常用指标理解。2、学习数字滤波器的设计方法。二、实验原理：1、数字滤波器设计的参数指标：w，阻带截止频率w ，通带增大衰减a设计数字滤波器的技术指标：通带截止频率和阻psp带的最大a。sw和阻带截止频率w为二元向量。在带通滤波器的设计中，通带截止频率ps2、滤波器的

30、设计中，数字频率与模拟频率的关系：数字频率 ，模拟频率为 f或 ，抽样频率 fs，归一化频率 f=f/fs或 = / s， = /2 ， 因为 f、 和 ， 都是无量纲， 刻度是一样的， 将 f、 、 、 f、 、 的定标值对应关系。 fs 1fsf 20fss0.500.51f -WWssss 10.50.51Ww-2-02w 10.500.513、 FIR数字滤波器 H(z)H(z) = b + b z + + b z-1-N01N 三、预习要求1. 在 MATLAB中，熟悉函数 fir1、kaiserord 、remezord、remez的使用；B = fir1(n,Wn,high,no

31、scale)设计滤波器；n,Wn,beta,ftype = kaiserord(f,a,dev)估计滤波器阶数；n,fo,ao,w = remezord (f,a,dev,fs)计算等波纹滤波器阶数 n和加权函数 w()；B=remez(n,f,a)进行等波纹滤波器的设计。2. 阅读扩展练习中的实例，学习 FIR滤波器的设计方法及其在 MATLAB中的实现；3. 给出 FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。4. 数字滤波器有 IIR和 FIR两种类型，它们的特点和设计方法不同。四、实验内容：1、利用 MATLAB编程,分别用窗函数法和等波纹滤波

32、器法设计两种 FIR数字滤波器，指标要求如下：w 0.45p w0.65p1= dB，通带峰值起伏：a；通带边缘频率：阻带边缘频率：，p1p2pw 0.3p w=0.8pa 40 dB。=，最小阻带衰减：s1s2sw 0.3pw 0.5p=，2、=用凯塞窗设计一 FIR低通滤波器，通带边界频率，阻带边界频率ps阻带衰减不小于 50dB。凯塞窗设计对应的 MATLAB参考程序为：wn=kaiser(30,4.55);nn=0:1:29;alfa=(30-1)/2;hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa);h=hd.*wn;h1,w1=freqz(h,1);或者

33、：b = fir1(29,0.4,kaiser(30,4.55);h1,w1=freqz(b,1);plot(w1/pi,20*log10(abs(h1);axis(0,1,-80,10); grid;xlabel(归一化频率/p);ylabel(幅度/dB);还可以使用n,Wn,beta,ftype=kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器：fcuts = 0.3 0.5; %归一化频率 omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs);

34、 %计算出凯塞窗 N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),noscale);freqz(hh);实际中，一般调用 MATLAB信号处理工具箱函数 remezord来计算等波纹滤波器阶数 N和加权函数 W()，调用函数 remez可进行等波纹滤波器的设计，直接求出滤波器系数。函数 remezord中的数组 fedge为通带和阻带边界频率，数组 mval是两个边界处的幅值，而数组 dev是通带和阻带的波动，fs是采样频率单位为 Hz。3、利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通 FIR数字滤波器，其指标为：通带边界频率 fc=80

35、0Hz，阻带边界 fr=1000Hz，通带波动 0.5dB，阻带最小衰减为 40dB，采样频率 fs=4000Hz。= 1-10 a = 0.0593=10 a = 0.01， 在 MATLAB中可以用 remezordp提示：d，d-/20-/20s12和 remez两个函数设计， MATLAB参考程序如下：fedge=800 1000;mval=1 0;dev=0.0593 0.01;fs=4000;N,fpts,mag,wt=remezord(fedge,mval,dev,fs);b=remez(N,fpts,mag,wt);h,w=freqz(b,1,256);plot(w*2000/pi,20*log10(abs(h);grid;xlabel(频率/Hz);ylabel(幅度/dB);