常见离散信号产生和实现.doc

数字信号处理实验报告实验1 常见离散信号产生和实现一、实验目的1、加深对常用离散信号的理解；2、熟悉使用MATLAB在时域中产生一些基本的离散时间信号。 二、实验原理MATLAB 语言提供了一系列函数用来产生信号，如exp, sin, cos, square, sawtooth，ones, zeros等函数。1. 基本信号序列1) 单位抽样序列 在MATLAB中可以利用zeros()函数实现。x=1 zeros(1, n-1)示范程序:% Program P1_1% Generation of a Unit Sample Sequence clf;% Generate a vector from -10 to 20n = -10:20;% Generate the unit sample sequenceu = zeros(1,10) 1 zeros(1,20);% Plot the unit sample sequencestem(n,u);xlabel(Time index n);ylabel(Amplitude);title(Unit Sample Sequence);axis(-10 20 0 1.2); 如果在时间轴上延迟了k个单位，得到即： 2) 单位阶跃序列 在MATLAB中可以利用ones()函数实现。3) 实指数序列MATLAB实现：4) 复指数序列MATLAB实现：5) 随机序列 MATLAB提供了两种随机信号： rand(1,N) 产生0,1上均匀分布的随机矢量。 randn(1,N) 产生均值为0，方差为1的高斯随机序列，即白噪声序列。2. 基本周期波形1) 方波MATLAB工具箱函数square可以产生方波;t=0:0.1*pi:6*pi;y=square(t);axis(0 7*pi -1.5 1.5);plot(t,y);xlabel(时间 t);ylabel(幅度y);2) 正弦波在MATLAB中3) 锯齿波工具箱函数sawtooth函数可以产生锯齿波Fs=10000;t=0:1/Fs:1.5; %抽样长度1.5s, 抽样频率为10kHzx=sawtooth(2*pi*50*t); %信号频率为50Hzplot(t, x);axis()0 0.2 -1 1; %画出0.2秒的波形3. 基本非周期波形工具箱函数chirp能产生一种扫射频率信号，其特点是信号的瞬时频率随时间按照一定规律变化 t=0:1/1000:2 %抽样频率1kHz, 抽样时间2s. x=chirp(t, 0.1, 150) %0时刻为DC信号，1s时频率为150Hz。 specgram(x, 256, 1000, 256, 250);4. sinc信号MATLAB实现： t=linspace(-5,5); x=sinc(t); plot(t,x);5. 序列的操作1) 信号加 x(n)=x1(n)+x2(n)MATLAB实现： x=x1+x2;注意：x1和x2序列应该具有相同的长度，位置对应，才能相加。2) 信号乘 x(n)=x1(n)*x2(n)MATLAB实现： x=x1.*x2; %数组乘法3) 改变比例 y(n)=k*x(n)MATLAB实现： y=k*x;4) 折叠 y(n)=x(-n);MATLAB实现： y=fliplr(x);5) 抽样和 MATLAB实现： y=sum(x(n1:n2);6) 抽样积 MATLAB实现： y=prod(x(n1:n2);7) 信号能量 MATLAB实现： Ex=sum(abs(x).2);8) 信号功率 MATLAB实现：Px=sum(abs(x).2)/ N;实验和程序：1)单位抽样序列 function uss(n)%构造函数 N=0:n-1 u = 1, zeros(1,n-1);%构造矩阵 stem(N,u); xlabel(Time index n);ylabel(Amplitude); title(Unit Sample Sequence); axis(-10 n 0 1.2); 右图为取n为10的图形2)单位阶跃序列function jieyue(n)%构造函数 N = -10:n-1;u = zeros(1,10) ones(1,n);stem(N,u);xlabel(Time index n);ylabel(Amplitude);title(Unit jump Sequence);axis(-10 n 0 1.2); 右图为取n为10的图形3)实指数序列function index(z1,N)%构造函数% n1=0:N;n=0:N/2;x=z1.n;stem(n,x);xlabel(Time index n);ylabel(Amplitude);右图为取z1=2,N=30的图形4)复指数序列function findex(z1,z2)%构造函数N=10;n=0:N/2;x=exp(z1+j*z2).*n);stem(n,abs(x);xlabel(Time index n);ylabel(Amplitude);右图为z1=2，z2=3的图形5)方波function square1(g)%G为方波的占空比 t=0:0.1:6*pi;y=square(t,g);plot(t,y);xlabel(时间 t);ylabel(幅度y);axis(0 7*pi -1.5 1.5); %要放在最后才可以正确显示右图分别为G=70，G=30的图形6）正弦波function sin1(N,A,Fs,f,Q)%N为显示 的范围，A为幅度， %FS为抽样频率，f为模拟信号的频率，Q为初相位n=0:N/2;x=A*sin(2*pi*f*n/Fs+Q)stem(n,x); xlabel(Time index n);ylabel(Amplitude);axis(0 N/2 -2.5 2.5);右图为sin1(60,2,1000,50,10)的图形7）锯齿波function saw(f,F)%f为锯齿波的频率，F为间隔 %其中F不能取太大，不易观 察图形 t=0:1/F:2;y=sawtooth(2*pi*f*t);plot(t,y);axis(0 2 -1.2 1.2);右图为saw(10,100)的图形8）基本非周期波形 %利用chirp来实现 t=0:1/500:2 %抽样频率500Hz, 抽样时间2s. x=chirp(t, 0.1, 150) %0时刻为DC信号，1s时频率为150Hz。 specgram(x, 256, 500, 256, 250);右图分别为抽样频率为500与10000的图形9）sinc信号 function sinc1(t1,t2,A)%t1,t2为定义域A为幅度 t=linspace(t1,t2); x=A*sinc(t); plot(t,x); 右图为t1,t2=-10，10，A=2的图形 10）序列相加 function sequenceplus(x1,x2) n1=length(x1); n2=length(x2); n=-n1/2:(n1/2-1);subplot(221) stem(n,x1) subplot(222) stem(n,x2) if(n1=n2) % x1，x2 的长度必须相同 x=x1+x2; subplot(223) stem(n,x); axis(-n1 n1 0 2) else error(x1与x2不一样长) end右图为x3=ones(1,20),ones(1,10),ones(1,10);x4=zeros(1,20),ones(1,10),zeros(1,10); sequenceplus(x3,x4)的图形如果输入x3=ones(1,10),ones(1,10),ones(1,10);x4=zeros(1,20),ones(1,10),zeros(1,10); 则会显示x1与x2不一样长11）序列相 乘function sequencenmultiply(x1,x2) n1=length(x1); n2=length(x2); n=-n1/2:(n1/2-1);subplot(221) stem(n,x1) subplot(222) stem(n,x2) if(n1=n2) % x1，x2的长度必须相同 x=x1.*x2; subplot(223) stem(n,x); axis(-n1 n1 0 2) else disp(error x1与x2不一样长) end右图为x3=ones(1,20),ones(1,10),ones(1,10);x4=zeros(1,20),ones(1,10),zeros(1,10); sequenceplus(x3,x4)的图形如果输入x3=ones(1,10),ones(1,10),ones(1,10);x4=zeros(1,20),ones(1,10),zeros(1,10); 则会显示x1与x2不一样长12）改变比例 function changerate(x,k) n1=length(x); %x为需要改变的序列 n=-n1/2:(n1/2-1); ，k为改变的比例 x1=k.*x; stem(n,x1);axis(-n1 n1 0 k+1);右图为x=ones(1,20),ones(1,10),ones(1,10);K=1.5的图形13)折叠 function pucker(x) n1=length(x); n=-n1/2:(n1/2-1);%这里的n必须关于0两边对称 x1=fliplr(x); ，不然就会出错 subplot(211) stem(n,x)%画出原序列 axis(-n1 n1 0 1.2); subplot(212) stem(n,x1);%画出折叠后的序列 axis(-n1 n1 0 1.2); 14)抽样和函数 n=0:24; x1=zeros(1,5),ones(1,3),zeros(1,5),ones(1,6),zeros(1,6);x=sum(x1)计算结果为x = 915)序列积 n=0:25; x1=rand(1,26)x=prod(x1)计算结果为x1 = Columns 1 through 6 0.9501 0.2311 0.6068 0.4860 0.8913 0.7621 Columns 7 through 12 0.4565 0.0185 0.8214 0.4447 0.6154 0.7919 Columns 13 through 18 0.9218 0.7382 0.1763 0.4057 0.9355 0.9169 Columns 19 through 24 0.4103 0.8936 0.0579 0.3529 0.8132 0.0099 Columns 25 through 26 0.1389 0.2028x = 4.6704e-01216) 序列能量 x3=ones(1,15),ones(1,20),ones(1,16); p=sum(abs(x3).2)计算结果为p = 5117)信号功率 x1=zeros(1,5),-2*ones(1,3),zeros(1,5),5*ones(1,6),zeros(1,6); N=length(x1) p=sum(abs(x1).2)/N计算结果为N = 25p =6.4800四、讨论复指数序列的性质。绘出信号，当、时、时的信号实部和虚部图；当时呢？此时信号周期为多少？调用复指数函数，绘出上面的z的图形如下所示（1）（2）（3）（4）（5）讨论：由上面的图形可以得出当z的实部为负数时，经过复指数运算后的实部和虚部都是振幅按指数衰减的正弦振荡。而z的实部为正数时，经过复指数运算后的实部和虚部都是振幅按指数增长的正弦振荡。对于实数经过复指数运算后只有实部，而当z为纯虚数时，经过复指数运算后的实部与虚部都是正弦振荡的，且振荡周期为虚部除以PI再2倍，上题中时，周期为12。实验2 离散系统的时域分析一、实验目的1、熟悉并掌握离散系统的差分方程表示法；2、加深对冲激响应和卷积分析方法的理解。 二、实验原理在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号,记系统单位冲激响应，则系统响应为如下的卷积计算式：当时，hn是有限长度的（），称系统为FIR系统；反之，称系统为IIR系统四、实验内容1、用MATLAB求系统响应1） 卷积的实现线性移不变系统可由它的单位脉冲响应来表征。若已知了单位脉冲响应和系统激励就可通过卷积运算来求取系统响应，即程序：function conv1(x,h) y=conv(x,h); N=length(y); n=0:1:N-1; disp(output sequence=); disp(y); stem(n,y); xlabel(Time index n);ylabel(Amplitude);输入为： x=-2 0 1 -1 3 h=1 2 0 -1function conv1(x,h)输出为：output sequence= -2 -4 1 3 1 5 1 -32） 单位脉冲响应的求取线性时不变因果系统可用MATLAB的函数filter来仿真 y=filter(b,a,x); 其中，x和y是长度相等的两个矢量。矢量x表示激励，矢量a，b表示系统函数形式滤波器的分母和分子系数，得到的响应为矢量y。例如计算以下系统的单位脉冲响应y(n)+0.7y(n-1)-0.45y(y-2)-0.6y(y-3)=0.8x(n)-0.44x(n-1)+0.36x(n-2)+0.02x(n-3) 程序实现： function filter1(N,b,a) % a=input(Type in the vector a=); N=input(Desired impuse response x=1 zeros(1,N-1); length=);%b=input(Type in the vector b=); y=filter(b,a,x); disp(output sequence) disp(y) k=0:1:N-1; stem(k,y); xlabel(Time index n); ylabel(Amplitude);输入数据：N=41b=0.8 -0.44 0.36 0.02a=1 0.7 -0.45 -0.6输出为：output sequence Columns 1 through 6 0.8000 -1.0000 1.4200 -0.9440 0.6998 -0.0627 Columns 7 through 12 -0.2076 0.5370 -0.5069 0.4719 -0.2363 0.0736 Columns 13 through 18 0.1253 -0.1964 0.2380 -0.1798 0.1151 -0.0187 Columns 19 through 24 -0.0430 0.0908 -0.0941 0.0809 -0.0445 0.0111 Columns 25 through 30 0.0207 -0.0362 0.0414 -0.0328 0.0198 -0.0038 Columns 31 through 36 -0.0081 0.0158 -0.0170 0.0142 -0.0081 0.0018 Columns 37 through 41 0.0036 -0.0065 0.0073 -0.0059 0.00352、以下程序中分别使用conv和filter函数计算h和x的卷积y和y1，运行程序，并分析y和y1是否有差别，为什么要使用xn补零后的x1来产生y1；具体分析当hn有i个值，xn有j个值，使用filter完成卷积功能，需要如何补零？clf;h = 3 2 1 -2 1 0 -4 0 3; %impulse responsex = 1 -2 3 -4 3 2 1; %input sequencey = conv(h,x);n = 0:14;subplot(2,1,1);stem(n,y);xlabel(Time index n); ylabel(Amplitude);title(Output Obtained by Convolution); grid; n1=length(h);x1 = x zeros(1,n1-1) %要补零,因为输入x1与输出y1要等长的,y1 = filter(h,1,x1)subplot(2,1,2);stem(n,y1);xlabel(Time index n); ylabel(Amplitude);title(Output Generated by Filtering); grid; max(y-y1) %判断两种方法求出的有没有差别以上程序段为有判断y与y1的差别的,用 conv 与filter求的结果一样,只是两个的原理不同.一个是卷积,一个是循环卷积.而对于循环卷积来说,它的两个对象必须长度相等,所以就需要进行补零当hn有i个值，xn有j个值，使用filter完成卷积功能，需要补零,首先判断哪个序列的长度更大,再将短的一个序列后面补零,使其与长序列相等长度.3、编制程序求解下列两个系统的单位冲激响应，分别用filter 和 impz实现，并绘出其图形。给出理论计算结果和程序计算结果并讨论。下面为程序段:function zzz%(1)yn+0.75yn-1+0.125yn-2=xn-xn-1%用filter求t1=1 0.75 0.125;%系统函数的分母t2=1 -1;%系统函数的分子subplot(2,1,1)filter(40,t2,t1)title(filter)%用impz求t1=1 0.75 0.125;%系统函数的分母t2=1 -1;%系统函数的分子subplot(212)impz(20,t2,t1)title(impz)%(2)yn=0.5xn-1+xn-2+xn-3+xn-4%用filter求t1=1%系统函数的分母t2=0.5 0.5 0.5 0.5%系统函数的分子subplot(211)filter1(40,t2,t1) title(filter)%用impz求t1=1t2=0.5 0.5 0.5 0.5subplot(212)filter1(40,t2,t1)title(impz)function filter1(N,b,a)% a=input(Type in the vector a=); N=input(Desired impuse response length=);%b=input(Type in the vector b=); x=1 zeros(1,N-1); y=filter(b,a,x); disp(output sequence) ; disp(y); k=0:1:N-1; stem(k,y); xlabel(Time index n); ylabel(Amplitude);% function impz1(N,b,a)% a=input(Type in the vector a=); N=input(Desired impuse response length=);%b=input(Type in the vector b=); % x=1 zeros(1,N-1); y1=impz(b,a,N); disp(output sequence) disp(y1); k=0:1:N-1; stem(k,y1);xlabel(Time index n); ylabel(Amplitude);第一,二题的图形如下图所示:理论计算的结果为(1) h(n)=-5*(-0.25)n*u(n)+6*(-0.5)n*u(n)(2) h(n)=0.25*d(n-1)+d(n-2)+d(n-3)+d(n-4) d代表单位冲激响应下面是用理论计算的结果会出的第一题的图形,与程序的结果大致相同.分析以上结果可以知道,用两种函数都可以求出差分方程的单位冲激响应.所求结果也大致相同.实验3 离散系统的变换域分析一、实验目的1、熟悉对离散系统的频率响应分析方法；2、加深对零、极点分布的概念理解。 二、实验原理离散系统的时域方程为其变换域分析方法如下：频域:系统的频率响应为:Z域:系统函数为:分解因式:，其中和称为零、极点。四、实验内容求系统的零、极点和幅度频率响应和相位响应。程序段如下;num=0.0528 0.0797 0.1295 0.1295 0.797 0.0528;%分子矩阵den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;%分母矩阵z,p,k=tf2zp(num,den);disp(零点);disp(z);disp(极点);disp(p);k=5000;%采样频率为5000Hzw=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,abs(h);grid;title(幅度谱,FS=5000);xlabel(omega/pi);ylabel(幅值);subplot(2,2,2);plot(w/pi,angle(h);grid; %w/pi为归一化的形式title(相位谱,FS=5000);xlabel(omega/pi);ylabel(弧度); k=100;%采样频率为100Hzw=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,3);plot(w/pi,abs(h);grid;title(幅度谱,FS=100);xlabel(omega/pi);ylabel(幅值);subplot(2,2,4);plot(w/pi,angle(h);grid; %w/pi为归一化的形式title(相位谱,FS=100);xlabel(omega/pi);ylabel(弧度); 运行结果如下: 零点 -1.0000 -0.3018 + 0.9534i -0.3018 - 0.9534i 0.0471 + 0.9989i 0.0471 - 0.9989i极点 0.2788 + 0.8973i 0.2788 - 0.8973i 0.3811 + 0.6274i 0.3811 - 0.6274i 0.4910 分样实验教材中的练习可以得出:频率响应的幅度特性与相频特性与H(Z)的零极点有关,极点位置主要影响频率响应的峰值位置及尖锐程度.零点主要影响频率响应的谷点位置与形状,对于一个因果稳定系统,其极点全部位于单位圆内.实验4 FFT算法的应用一、实验目的1、加深对离散信号的DFT的理解；2、在MATLAB中实现FFT算法。二、实验原理N点序列的DFT和IDFT变换定义式如下：,利用旋转因子具有周期性，可以得到快速算法（FFT）。在MATLAB中，可以用函数和计算N点序列的DFT正、反变换。四、实验内容1、2N点实数序列N=64。用一个64点的复数FFT程序，一次算出，并绘出 的图形。2、已知某序列在单位圆上的N=64等分样点的Z变换为:。用N点IFFT程序计算出。实验内容的程序段如下:(1) format compactN=64;x=zeros(1,128)for n=1:2*N;x(n)=cos(2*pi*7*(n-1)/N)+0.5*(2*pi*19*(n-1)/N);endxs=zeros(1,64);%定义两个数组用来存放X的奇偶项xx=zeros(1,64);for n1=1:N %X分为奇偶项 xs(n1)=x(2*n1-1); xx(n1)=x(2*n1);endx1=xs+j*xx; %定义一个新的复序列,长度为64X1=fft(x1);X11=conj(X1); %取共轭for k=1:N Xep(k)=0.5*(X1(k)+X11(N+1-k) %共轭对称项 Xop(k)=0.5*(X1(k)-X11(N+1-k) %共轭反对称项 end Xs=Xep;%偶数项对应的 Xx=-j*Xop;%奇数项对应的for n3=1:N X(n3)=Xs(n3)+exp(-j*2*pi*n3/(2*N)*Xx(n3); X(n3+N)=Xs(n3)-exp(-j*2*pi*n3/(2*N)*Xx(n3); endn4=0:127; stem(n4,abs(X) 注：本图为实验程序所得图形2）注：本图为直接采用128点FFT算法得出的图形比较上面两图形，发现采用基2FFT算法求出的IDFT在N=64左右会出现误差，即出现频谱泄漏这是由于将原序列分为两个64点的序列再进行复数的FFT运算造成的，也就是截断效应还有混叠现象造成的。format compactN=64; n1=0:N-1;for n=1:N X3=1./(1-0.8*exp(-j*2*pi*(n-1)/N);end x3=ifft(X3,N)stem(n1,x3)在序列进行FFT算法时，在截取序列时应该要2的整数次冥方。不够的要补零。还有就是要关于谱分析范围与频率分分辨率，这就要求采样频率FS2fc,而分辨率F与采样间隔有关，故在选择FS时要使N2fc/F,NT=1/F;这样就会使分析的谱更理想。在分析谱时，如果不知道序列的周期，可以先截取M点进行DFT，求出一个X1（K），再将截取长度扩大一倍，再进行DFT，求出X2（K）；如果二者的主谱判别满足分析误差要求，则以其中一个近似表示x(n)的频谱。否则继续取长序列，直到前后两次分析的谱满足误差要求。实验5 无限冲激响应数字滤波器设计一、实验目的1、掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法;2、熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。二、实验原理在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1)利用buttord和cheb1ord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；如：求阶数N, Wn = cheb1ord (Wp, Ws, Rp, Rs,s)选择项说明：high类别。缺省为low； s模/数，缺省为数Rp即 ap, Rs即as ； WnChebyshev自然频率（3dB频率），数字设计： Wpwp/p Wsws/p 。2)num,den=butter（N,Wn）(巴特沃斯)和num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；3）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；4)使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；5)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。四、实验内容利用MATLAB编程，用脉冲响应不变法和双线性变换法设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：；阻带边缘频率：，最小阻带衰减： 。程序段如下:%采用切比雪夫%脉冲响应format compactfs=1000;%采样频率wp1=0.45*pi*fs;wp2=0.65*pi*fs;ws1=0.3*pi*fs;ws2=0.8*pi*fs;N,wn=cheb1ord(wp1 wp2,ws1 ws2,1,40,s);B,A=cheby1(N,1,wn,s);num,den=impinvar(B,A,fs);h1,w=freqz(num,den); %双线性法wp3=2*fs*tan(pi*0.45/2); wp4=2*fs*tan(pi*0.65/2);ws3=2*fs*tan(pi*0.3/2);ws4=2*fs*tan(pi*0.8/2);N,wn=cheb1ord(wp1 wp2,ws1 ws2,1,40,s);B,A=cheby1(N,1,wn,s);num,den=bilinear(B,A,fs);h2,w=freqz(num,den);f=w/pi*fs/2;plot(f,20*log10(abs(h1),-.,f,20*log10(abs(h2),-);axis(0,600,-80,10);grid;title(采用切比雪夫)xlabel(f)ylabel(幅度/dB)%采用巴特沃斯format compactfs=1000;%采样频率wp1=0.45*pi*fs;wp2=0.65*pi*fs;ws1=0.3*pi*fs;ws2=0.8*pi*fs;N,wn=buttord(wp1 wp2,ws1 ws2,1,40,s);B,A=butter(N,wn,s);num,den=impinvar(B,A,fs);h1,w=freqz(num,den); %双线性法wp3=2*fs*tan(pi*0.45/2); wp4=2*fs*tan(pi*0.65/2);ws3=2*fs*tan(pi*0.3/2);ws4=2*fs*tan(pi*0.8/2);N,wn=buttord(wp1 wp2,ws1 ws2,1,40,s);B,A=butter(N,wn,s);num,den=bilinear(B,A,fs);h2,w=freqz(num,den);f=w/pi*fs/2;plot(f,20*log10(abs(h1),-.,f,20*log10(abs(h2),-);axis(0,600,-80,10);grid;title(采用巴特沃斯)xlabel(f)ylabel(幅度/dB) 采用切比雪夫与采用巴特沃斯产生的滤波器频谱如