布莱克曼哈尔窗matlab,数字信号处理实验报告

数字信号处理实验报告

数字信号处理实验报告姓 名： 班 级： 13电信2 学 号： 2013302 2013302 2013302指导老师： 日期： 2016.6.62016.6.17 华南农业大学电子工程学院电子信息工程系实验一 常见离散信号的MATLAB产生和图形显示一、实验目的加深对常用离散信号的理解；二.实验原理1 单位抽样序列在MATLAB中可以利用zeros()函数实现。如果在时间轴上延迟了k个单位，得到即：2 单位阶越序列在MATLAB中可以利用ones()函数实现。3 正弦序列在MATLAB中4 复正弦序列在MATLAB中5 指数序列在MATLAB中6.卷积分析conv.m用来实现两个离散序列的线性卷积。其调用格式是：y=conv(x,h)若x的长度为N，h的长度为M，则y的长度L=N+M-1。三.实验内容1.画出信号x(n) = 1.5*d(n+1) - d(n-3)的波形。2.求序列x(n)和h(n)的线性卷积y(n)=x(n)*h(n)。x(n) = 3,-3,7,0,-1,5,2 , h(n) = 2,3,0,-5,2,1. 画出x(n),h(n),y(n)与n的离散序列图形四.实验要求1)画出信号x(n) = 1.5*d(n+1) - d(n-3)的波形。MATLAB程序如下：n3 = -3:3;x3 = (n3+1)=0;subplot(1,3,1);stem(n3,x3);n4 = -3:3;x4 = (n4-3)=0;subplot(1,3,2);stem(n4,x4);n5 = -3:3;x5 = 1.5*x3 - x4;subplot(1,3,3);stem(n5,x5);理论计算：x(n)=1.5 n=-1 -1 n=3 0 else 程序运行结果：图(1)从图(1)左侧起第一幅图是信号d(n+1)的波形，第二幅图是信号d(n-3)的波形，最后一幅图是信号x(n) = 1.5*d(n+1) - d(n-3)的波形。2)求序列x(n)和h(n)的线性卷积y(n)=x(n)*h(n)；x(n) = 3,-3,7,0,-1,5,2 , h(n) = 2,3,0,-5,2,1，画出x(n),h(n),y(n)与n的离散序列图形。MATLAB程序如下：n6 = 0:6;x6 = 3,-3,7,0,-1,5,2;subplot(1,3,1);stem(n6,x6);n7 = 0:5;x7 = 2,3,0,-5,2,1;subplot(1,3,2);stem(n7,x7);n8 = 0:11;x8 = conv(x6,x7);subplot(1,3,3);stem(n8,x8);理论分析：信号 x(n)的长度为 7，,h(n)的长度为 6，则线性卷积y(n)=x(n)*h(n)的长度为 7+6-1=12。y(n) = 6 ,3 ,5 ,6,-7 ,-25 ,30 ,21 ,-23 ,-1 ,9 ,2程序运行结果：图(2)从图(2)左侧起第一幅图是信号x(n)的波形，第二幅图是信号h(n)的波形，最后一幅图是线性卷积信号y(n)=x(n)*h(n)的波形。经过比较，理论与实验结果一致。实验二 离散系统的差分方程、冲激响应和卷积分析一、实验目的加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。二.实验原理离散系统其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号，。记系统单位冲激响应，则系统响应为如下的卷积计算式：当时，hn是有限长度的(n：0，M)，称系统为FIR系统；反之，称系统为IIR系统。1. filter可用来求一个离散系统的输出。调用格式：y=filter(b,a,x);2. impz可用来求一个离散系统的h(n)。调用格式：h=impz(b,a,N);h,t=impz(b,a,N);三.实验内容编制程序求解下列两个系统的输出、单位冲激响应，并绘出其图形。；四.实验要求给出理论计算结果和程序计算结果并讨论。1)理论结果：解：(1)对于计算可得系统函数为：H(z)=+由此可得h(n)=-5(-0.25)+6(-0.5)u(n)h(0)=1, h(1)=-1.75, h(2)=1.1875, h(3)=-0.828125程序计算结果：MATLAB程序如下：m = -30:30;b =1,-1;a = 1,0.75,0.125; x9 = (m-0)=1; h =filter(b,a,x9);n = (-10:50);subplot(1,2,1);stem(n,h);axis(-10,50,-1,1.5);title(Impluse Response);xlabel(n),ylabel(h(n);subplot(1,2,2);y1 = impz(b,a,n);stem(n,y1);程序运行结果：图3输出函数、单位冲激响应2) 理论结果：(2)对于由差分方程计算系统函数可得H(z)=0.25()则单位冲击响应为：h(n)=0.25(n-1)+(n-2)+(n-3)+(n-4),该波形与仿真结果一致。程序计算结果：MATLAB程序如下：m = -30:30;b =0.25,0.25,0.25,0.25;a =1; x9 = (m-0)=1; h =filter(b,a,x9);n = (-10:50);subplot(1,2,1);stem(n,h);axis(-10,50,-1,1.5);title(Impluse Response);xlabel(n),ylabel(h(n);subplot(1,2,2);y1 = impz(b,a,n);stem(n,y1);程序运行结果：图 4输出函数、单位冲激响应实验三、离散系统的零、极点分布及频率响应分析一、实验目的加深对离散系统的频率响应分析和零、极点分布的概念理解。二.实验原理离散系统的时域方程为freqz可用来求一个离散系统的频率响应。调用格式：H,w = freqz(B,A,N,whole);B和A分别为离散系统的(系统函数分子、分母多项式的系数向量)；N 表示选取单位圆的上半圆等间距的N个点作为频响输出；返回量H则包含了离散系统频响在0pi范围内N个频率等分点的值，向量w则包含范围内N个频率等分点。调用中若N默认，默认值为512。三.实验内容已知系统差分方程如下：y(n)-1.6y(n-1)+1.28y(n-2) =0.5x(n)+0.1x(n-1)四.实验要求(1)编程得到系统频响的幅度响应和相位响应，并画图。(2)编程得到系统的零极点分布图，分析系统的因果性和稳定性。(3)给出理论计算结果和程序计算结果并讨论。提示：幅度响应mag=abs(H)，plot(w，mag)相位响应ph=angle(H)，plot(w，ph)五实验结果分析处理y(n)-1.6y(n-1)+1.28y(n-2) =0.5x(n)+0.1x(n-1)1)编程得到系统频响的幅度响应和相位响应，并画图。MATLAB程序如下：n=-30:30;b = 0.5,0.1,a =1,-1.6,1.28;H,w = freqz(b,a,512,whole);mag = abs(H),subplot(1,2,1);plot(w,mag);ph = angle(H),subplot(1,2,2);plot(w,ph);MATLAB运行结果：图1 系统频响的幅度响应和相位响应(2)编程得到系统的零极点分布图，分析系统的因果性和稳定性。MATLAB程序如下：n=-30:30;b = 0.5,0.1,a =1,-1.6,1.28;z,p,K = tf2zp(b,a);zplane(z,p);MATLAB运行结果：图2 系统的零极点分布图理论计算分析：解：由理论计算可得H(z)=,解得：极点为：z1=0.8+0.8j, z2=0.8-0.8j。零点为:z=-0.2其中|z|=|z1|=|z2|=0.8=1.1312,当|z|=0.8时，收敛域包括无穷大，系统是因果系统；单位圆不在敛域内面，所以系统不稳定。当|z|0.8时，收敛域不包括无穷大，系统不是因果系统；单位圆在收敛域内面，所以系统稳定。实验四、DFT算法的应用一、实验目的利用DFT对信号(如由多个正弦信号组成的信号)进行频谱分析，并研究采样长度、截断(即加窗)和补零对分析频率的影响。二.实验原理1. 截断：通常情况下，信号都是无限长的。而在运用计算机进行模拟时，这是无法操作的。所以实际情况下，要把观测的信号限制在一定长的时间之内。为了从无限长的信号中得到有限长的数据，在时域乘一个窗函数，将信号截短，叫做加窗。2. 补零：为了增加频域抽样点数N，在不改变时域数据的情况下，在时域数据末端加一些零值点，叫做补零。3. 频率分辨率：指对两个最近的频谱峰值能够分辨的能力。利用DFT进行频谱分析，并研究不同数据长度，补零，加窗等对频率分辨率的影响。利用DFT计算频谱的目的在于，针对计算机只能计算有限个离散的点的取值这一特点，实现计算机对连续时间信号的频谱的模拟。所以我们比较关心的是模拟频谱和原信号频谱的拟合程度，我们希望拟合程度越高越好。这就需要增加频率分辨率，因为频率分辨率越高，根据公式，说明相同采样频率下，采样的长度就越长，也就是频谱采样的点数就越多，我们可以看到的模拟频谱图像就越清晰，这样与原信号的拟合程度就越好。MATLAB提供了2个内部函数用于计算DFT，它们分别是：fft(x,N)fft(x,N) 计算N点的DFT。三.实验内容 已知序列x=sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs)，fs=30(1)不同采样长度的影响：分别计算序列的32点和128点DFT，绘出幅度谱图形。F32=fft(x,32);F32mag=abs(F32);(2)补零对频率分辨率的影响：设定序列长度为32，采样长度64，对连续时间信号进行采样，首先对采样的信号时域图像后补若干个零，然后计算出采样长度情况下的频谱，并给出32点与64点的序列图与幅度谱图。m1=64;m=0:m1-1x2=x1,zeros(1,m1-n1)(3)加窗对频率分辨率的影响：一部分是设定采样长度N1，对连续时间信号进行采样，然后利用DFT计算出采样长度N1情况下的频谱，并给出时域和频域图像。一部分采用加窗的方法对时域图像进行采样，然后利用DFT计算出采样长度N1情况下的频谱，并给出时域和频域图像。这里要注意的是我们采用加窗的方法时，设定信号的长度为400，然后对长度为400的信号进行截短加窗，这时只有加窗处我们是对信号采样的，即得到的信号是有效的，其余位置信号都为0，所以这时我们仍然认为采样长度为128. 加窗：n=0:399;m=0:399;m2=128;fs=30;x= sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs);x1=x.*(heaviside(m)-heaviside(m-m2);F1=fft(x1);mag=abs(F1);stem(mag);四实验结果分析处理已知序列x=sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs)，fs=30(1)不同采样长度的影响：分别计算序列的32点和128点DFT，绘出幅度谱图形。MATLAB程序如下：fs = 30;n = 0:31;x1 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs);F32 = fft(x1,32);F32mag = abs(F32);subplot(1,2,1);stem(n,F32mag); m = 0:127;x2 = sin(5*2*pi*m/fs) + cos(3*2*pi*m/fs);F128 = fft(x2,128);F128mag = abs(F128);subplot(1,2,2);stem(m,F128mag);MATLAB运行结果：图1 32点和128点DFT幅度谱图形(2)补零对频率分辨率的影响：设定序列长度为32，采样长度64，对连续时间信号进行采样，首先对采样的信号时域图像后补若干个零，然后计算出采样长度情况下的频谱，并给出32点与64点的序列图与幅度谱图。MATLAB程序如下：fs = 30;n = 0:31;x1 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs);F32 = fft(x1,32);F32mag = abs(F32);subplot(1,2,1);stem(n,x1);subplot(1,2,2);stem(n,F32mag);m1= 64;n1 = 32;x2=x1,zeros(1,m1-n1);F64 = fft(x2,64);F64mag = abs(F64);subplot(1,2,1);stem(n2,x2);subplot(1,2,2);stem(n2,F64mag);MATLAB运行结果：图2补零后32点与64点的序列图与幅度谱图(3)加窗对频率分辨率的影响：一部分是设定采样长度N1，对连续时间信号进行采样，然后利用DFT计算出采样长度N1情况下的频谱，并给出时域和频域图像。一部分采用加窗的方法对时域图像进行采样，然后利用DFT计算出采样长度N1情况下的频谱，并给出时域和频域图像。MATLAB程序如下：fs = 30;n = 0:399;m = 0:399;m2=128;x3 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs);x5=ones(1,128),zeros(1,272);x4 = x3.*x5;F1 = fft(x4);mag = abs(F1);stem(mag);MATLAB运行结果：图3 400点DFT序列图分析：由结果图可知，加穿函数后信号的频谱图近似等于信号原谱图。实验五：IIR滤波器的设计一、实验目的1理解滤波器参数的意义；2掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的方法；3掌握利用Matlab设计其它各型IIR数字滤波器的方法；4掌握分析滤波器是否达到性能指标的方法。二、实验原理1.利用脉冲响应不变法，直接根据归一化的巴特沃斯低通模拟滤波器系统函数H(p)得到IIR数字低通滤波器方法是IIR数字滤波器的设计步骤：(1)按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标；模拟通带、阻带截止频率(2)根据模拟滤波器技术指标设计为响应的模拟低通滤波器；n,Wn=buttord(wp,ws,Ap,As,s)z,p,k=buttap(n);设计模拟滤波器b,a=zp2tf(z,p,k); 将系统函数的零极点转化为系统函数一般形式的系数bn,an=lp2lp(b,a,Wn);(3)根据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器；bz,az=bilinear(bn,an,1/Ts)bz,az= impinvar (bn,an,1/Ts)(4)如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到所要的滤波器。3.直接设计数字滤波器方法：n,Wn=buttord(wp/pi,ws/pi,Ap,As)b,a=butter(n,Wn,high)设计其它各型IIR数字滤波器的理论方法在这里不再给出，读者可参看有关内容。在Matlab中，设计滤波器的过程很简单，只要加上一些控制字符即可。控制字符省略或为“low”表示设计低通滤波器，控制字符为“high”表示设计高通滤波器，控制字符为“band”表示设计带通滤波器，控制字符为“stop”表示设计带阻滤波器。三、实验内容freqz(b,a);figure;y,t=impz(bz,az,101);stem(t,y)四、实验结果(1)脉冲响应不变法设计： MATLAB程序：Ts = 0.001;n1,Wn1 = buttord(0.2*pi/Ts,0.6*pi/Ts,2,15,s),z,p,k = buttap(n1), b,a = zp2tf(z,p,k);b1,a1 = lp2lp(b,a,Wn1);bz1,az1 = impinvar(b1,a1,1000);freqz(bz1,az1); figure;y1,t1 = impz(bz1,az1,101);stem(t1,y1); MATLAB仿真结果：图1 As=15dB频率响应截图图2 冲击响应图 理论计算分析：p = 0.2pi 时有 Ap = 2dB ；s = 0.6pi 时有 As 15dB，符合设计要求。(2)双线性变换法设计 MATLAB程序：Ts = 0.001;n1,Wn1 = buttord(0.2*pi/Ts,0.6*pi/Ts,2,15,s),z,p,k = buttap(n1),b,a = zp2tf(z,p,k);b1,a1 = lp2lp(b,a,Wn1);bz1,az1 = bilinear(b1,a1,1000);freqz(bz1,az1);figure;y1,t1 = impz(bz1,az1,101);stem(t1,y1); MATLAB仿真结果：图3 As=18dB频率响应截图图4 冲击响应图 理论分析：p = 0.2pi 时有 Ap = 2dB ；s = 0.6pi 时有 As 15dB，符合设计要求。(3)直接设计法设计高通滤波器 MATLAB程序：Ts = 0.001;n1,Wn1 = buttord(0.8,0.44,3,20),b,a = butter(n1,Wn1,high);freqz(b,a);figure;y1,t1 = impz(b,a,101);stem(t1,y1); MATLAB仿真结果：图6 As=20dB频率响应截图图7 冲击响应截图 理论分析：p = 0.8pi 时有 Ap = 3dB ；s = 0.44pi 时有 As 20dB，符合设计要求。实验六、FIR滤波器的设计一、实验目的1 掌握FIR数字滤波器的设计方法；2 熟悉利用MATLAB软件进行FIR数字滤波器设计，以及对所设计的滤波器进行分析；二、实验原理窗函数法的设计步骤窗函数设计FIR数字滤波器是从时域出发的，把理想滤波器的单位取样响应用于合适的窗函数截短成为有限长度的h(n),使得h(n)逼近理想的hd(n)。以实现所设计的滤波器的频率响应H(ejw)逼近与理想滤波器的频率响应Hd(ejw)。由过渡带宽度及阻带最小衰减要求，可选定窗形状，并估计窗口长度N。Matlab提供了fir1函数，以实现线性相位FIR滤波器。调用格式如下：b=fir1(N，wc，ftype，window(M)参数：N：阶数wc：归一化的数字频率，0wc1 wc=(wp+(ws-wp)/2)/piftype：滤波器类型，如高通、带阻等。window:应用的窗函数类型实验内容基于窗函数设计法，利用MATLAB软件设计满足设计要求的FIR数字低通滤波器。(1)要求利用窗函数设计法设计FIR数字低通滤波器，滤波器参数要求为： 。要求设计数字低通滤波器，同时要求给出滤波器的辐频特性曲线图以及幅度响应图和对数幅度响应图；参考函数：freqz(b,a)；h,w=freqz(b,a);plot(w,abs(h)plot(w,20*log10(abs(h)三、实验结果(1) 汉宁窗MATLAB程序：wp = 0.2*pi;ws = 0.3*pi;a = 1;wc = (wp+(ws-wp)/2)/pi;b = fir1(62,wc,low,hanning(63);freqz(b,a);figure;h,w = freqz(b,a);plot(w,abs(h);figure;plot(w,20*log10(abs(h);MATLAB运行结果：图1 辐频特性曲线图图2 幅度响应图图3 对数幅度响应图理论分析计算：当wp=0.2*pi时，Ap=0.25dB，当ws=0.3*pi，As=50dB。阻带衰减符合要求。(2) 海明窗MATLAB程序：wp = 0.2*pi;ws = 0.3*pi;a = 1;wc = (wp+(ws-wp)/2)/pi;b = fir1(66,wc,low,hamming(67);freqz(b,a);figure;h,w = freqz(b,a);plot(w,abs(h);figure;plot(w,20*log10(abs(h);MATLAB运行结果：图4 辐频特性曲线图图5 幅度响应图图6 对数幅度响应图理论分析计算：当wp=0.2*pi时，Ap=0.25dB，当ws=0.3*pi，As50dB。阻带衰减稍微不符合要求。(3) 布莱克曼窗MATLAB程序：wp = 0.2*pi;ws = 0.3*pi;a = 1;wc = (wp+(ws-wp)/2)/pi;b = fir1(110,wc,low,Blackman(111);freqz(b,a);figure;h,w = freqz(b,a);plot(w,abs(h);figure;plot(w,20*log10(abs(h);MATLAB运行结果：图7 辐频特性曲线图图8 幅度响应图图9 对数幅度响应图理论分析计算：当wp=0.2*pi时，Ap=0.25dB，当ws=0.3*pi，As=50dB。阻带衰减符合要求，甚至远大于要求衰减。四、实验分析窗函数的选择原则是保证阻带衰减的情况下选择主瓣窄的窗函数。由于滤波器的参数要求是，海明窗的阻带最小衰减为-53db，满足该滤波器阻带衰减要求，海宁窗为-44db，不满足要求，布拉克曼窗为-74db也满足要求，理论上说海明窗最适合设计该滤波器。由程序结果可知，海宁窗和海明窗的主瓣宽度一致，而布莱克曼窗主瓣宽度比他们大，但是阻带衰减最大。精选文档