数字信号处理MATLAB仿真.doc

数字信号处理MATLAB仿真 作者： 日期：13 个人收集整理 勿做商业用途实验一 数字信号处理的Matlab仿真一、实验目的1、掌握连续信号及其MATLAB实现方法；2、掌握离散信号及其MATLAB实现方法3、掌握离散信号的基本运算方法，以及MATLAB实现4、了解离散傅里叶变换的MATLAB实现5、了解IIR数字滤波器设计6、了解FIR数字滤波器设计1二、实验设备计算机，Matlab软件三、实验内容(一）、 连续信号及其MATLAB实现1、 单位冲击信号例1.1:t=1/A=50时，单位脉冲序列的MATLAB实现程序如下：clear all;t1=0.5:0。001：0；A=50；A1=1/A；n1=length(t1）；u1=zeros(1，n1);t2=0：0。001：A1；t0=0；u2=A*stepfun（t2，t0);t3=A1：0。001：1；n3=length（t3)；u3=zeros(1,n3)；t=t1 t2 t3;u=u1 u2 u3；plot(t，u)axis(0。5 1 0 A+2）2、 任意函数例1.2：用MATLAB画出如下表达式的脉冲序列clear all;t=2：1:3；N=length（t)；x=zeros（1，N);x（1）=0。4；x(2)=0。8x（3）=1。2;x(4）=1.5；x(5）=1.0；x(6)=0.7；stem（t,x);axis(-2。2 3。2 0 1.7）3、 单位阶跃函数例1。3：用MATLAB实现单位阶跃函数clear all;t=0。5：0.001：1;t0=0；u=stepfun（t,t0）；plot(t，u）axis(-0.5 1 0.2 1。2）4、 斜坡函数例1.4:用MATLAB实现g(t）=3（t1)clear all;t=0:0。01:3;B=3；t0=1；u=stepfun(t，t0）;n=length（t）;for i=1:n u（i）=B*u（i)(t(i）-t0）;endplot（t,u)axis（0。2 3。1 0。2 6.2）5、 实指数函数例1.5:用MATLAB实现clear all;t=0:0.001:3；A=3;a=0.5；u=A*exp（at);plot(t，u)axis(-0。2 3.1 -0。2 14)6、 正弦函数例1。6:用MATLAB实现正弦函数f（t)=3cos（10t+1）clear all；t=-0。5：0.001：1；A=3；f=5;fai=1；u=Asin（2*pi*ft+fai）;plot（t,u）axis(0.5 1 -3。2 3.2）（二）、离散信号及其MATLAB实现1、 单位冲激序列例2.1：用MATLAB产生64点的单位冲激序列clear all;N=64;x=zeros（1,N)；x(1)=1;xn=0：N1；stem(xn,x)axis(1 65 0 1。1)2、 任意序列例2.2：用MATLAB画出如下表达式的脉冲序列clear all;N=8；x=zeros(1,N）；x（1）=8。0;x（2）=3.4x（3）=1。8;x(4)=5.6;x(5)=2。9;x（6）=0.7;xn=0：N1;stem（xn，x)axis(1 8 0 8.2)3、 单位阶跃序列例2.3:用MATLAB实现单位阶跃函数clear all；N=32；x=ones(1，N);xn=0:N-1；stem（xn，x)axis（1 32 0 1.1)4、 斜坡序列例2。4：用MATLAB实现g（n)=3（n-4)点数为32的斜坡序列clear all;N=32；k=4B=3；t0=1；x=zeros(1,k) ones(1,N-k)；for i=1：N x（i)=B*x（i）（i-k）；endxn=0：N-1;stem（xn，x)axis(1 32 0 90）5、 正弦序列例2。5:用MATLAB实现幅度A=3，频率f=100，初始相位=1.2,点数为32的正弦信号clear all；N=32；A=3;f=100；fai=1.2;xn=0：N-1;x=Asin（2pif*(xn/N)+fai）；stem(xn,x）axis(-1 32 3.2 3.2）6、 实指数序列例2。6：用MATLAB实现,点数为32的实指数序列clear all；N=32;A=3；a=0。7；xn=0:N1；x=Aa.xn；stem（xn,x）7、 复指数序列例2.7:用MATLAB实现幅度A=3,a=0。7，角频率=314，点数为32的实指数序列clear all；N=32；A=3；a=0。7;w=314；xn=0：N1；x=A*exp(a+j*w)*xn);stem(xn，x）8、 随机序列利用MATLAB产生两种随机信号：rand(1,N）在区间上产生N点均匀分布的随机序列randn（1，N）产生均值为0，方差为1的高斯随机序列，即白噪声序列例2。8:用MATLAB产生点数为32的均匀分布的随机序列与高斯随机序列clear all;N=32;x_rand=rand(1,N）；x_randn=randn(1,N）;xn=0：N1；figure（1)；stem（xn,x_rand)figure(2）;stem(xn，x_randn）（三）、离散信号的基本运算1、 信号的延迟给定离散信号x(n），若信号y（n）定义为：y(n）=x(n-k)，那么y（n)是信号x（n)在时间轴上右移k个抽样周期得到的新序列。例3.1：正弦序列y(n)=sin(100n)右移3个抽样周期后所得的序列，MATLAB程序如下：clear all;N=32;w=100;k=3；x1=zeros(1,k);xn=0：N1;x2=sin(100*xn);figure（1)stem（xn，x2)x=x1 x2；axis(1 N -1。1 1.1)N=N+k;xn=0：N-1;figure（2）stem（xn,x)axis（-1 N -1.1 1.1)2、 信号相加若信号，值得注意的是当序列和的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加例3。2:用MATLAB实现两序列相加clear all；n1=0:3x1=2 0.5 0。9 1；figure（1)stem(n1,x1）axis(-1 8 0 2。1 ）n2=0：7x2= 0 0。1 0。2 0。3 0.4 0.5 0。6 0。7；figure(2）stem(n2，x2）axis（-1 8 0 0.8 ）n=0:7;x1=x1 zeros（1,8-length(n1）)；x2= zeros（1,8length（n2）,x2;x=x1+x2;figure(3)stem（n，x）axis(1 8 0 2.1）3、 信号相乘信号序列和相乘所得信号的表达式为：这是样本与样本之间的点乘运算，在MATLAB中可采用“。来实现,但是在信号序列相乘之前,应对其做与相加运算一样的操作。例3.3：用MATLAB实现上例中两序列相乘clear all；n1=0:3x1=2 0.5 0.9 1；figure（1）stem(n1,x1)axis（-1 8 0 2。1 ）n2=0：7x2= 0 0。1 0.2 0.3 0.4 0。5 0.6 0.7；figure(2）stem(n2,x2）axis(-1 8 0 0。8 ）n=0：7;x1=x1 zeros（1,8-length(n1)；x2= zeros(1,8-length(n2)），x2;x=x1.x2;figure（3)stem(n，x）axis（1 8 0 0.35)4、 信号翻转信号翻转的表达式为：y（n)=x（n）,在MATLAB中可以用fliplr函数实现此操作例3.4：用MATLAB实现“信号相加”中的序列翻转clear all;n=0:3x1=2 0。5 0。9 1；x=fliplr(x1);stem（n，x）axis(1 4 0 2。1 ）5、 信号和对于N点信号，其和的定义为：例3.5:用MATLAB实现“信号相加”中的序列和clear all；n=0:3x1=2 0.5 0。9 1；x=sum（x1)6、 信号积对于N点信号,其积的定义为：例3。5：用MATLAB实现“信号相加”中的序列积clear all；n=0：3x1=2 0。5 0。9 1;x=prod（x1)（四）、离散傅里叶变换的MATLAB实现例:若是一个N=32的有限序列,利用MATLAB计算它的DFT并画出图形。N=32;n=0：N1;xn=cos（pi*n/6）;k=0：N-1；WN=exp（j2*pi/N)；nk=n*k;WNnk=WN.nk；Xk=xnWNnk;figure（1)stem（n，xn）figure(2)stem（k,abs（Xk）)在MATLAB中，可以直接利用内部函数fft来实现FFT算法，该函数是机器语言，而不是MATLAB指令写成的，执行速度很快。常用格式为：y=fft（x）y=fft（x,N)（五）、IIR数字滤波器设计1、 基于巴特沃斯法直接设计IIR数字滤波器例。1：设计一个10阶的带通巴特沃斯数字滤波器，带通频率为100Hz到200Hz，采样频率为1000Hz，绘出该滤波器的幅频于相频特性，以及其冲击响应图clear all；N=10;Wn=100 200/500;b,a=butter（N，Wn，bandpass）；freqz(b,a，128，1000）figure（2)y，t=impz(b,a，101)；stem(t，y)2、 基于切比雪夫法直接设计IIR数字滤波器例5。2：设计一个切比雪夫型数字低通滤波器，要求：Ws=200Hz,Wp=100Hz,Rp=3dB,Rs=30dB，Fs=1000Hzclear all；Wp=100；Rp=3；Ws=200;Rs=30；Fs=1000；N,Wn=cheb1ord（Wp/(Fs/2）,Ws/（Fs/2）,Rp,Rs）;b，a=cheby1(N，Rp,Wn);freqz（b，a,512,1000);例5。3：设计一个切比雪夫型数字带通滤波器，要求带通范围100250Hz，带阻上限为300Hz，下限为50Hz,通带内纹波小于3dB，阻带纹波为30 dB，抽样频率为1000 Hz,并利用最小的阶次实现。clear all;Wpl=100;Wph=250；Wp=Wpl，Wph;Rp=3;Wsl=50;Wsh=300；Ws=Wsl，Wsh；Rs=30；Fs=1000；N，Wn=cheb2ord(Wp/(Fs/2)，Ws/（Fs/2),Rp，Rs）；b，a=cheby2（N,Rp,Wn);freqz（b,a,512,1000);（六)、FIR数字滤波器设计1、在MATLAB 中产生窗函数十分简单：（1）矩形窗（Rectangle Window）调用格式:w=boxcar（n)，根据长度n 产生一个矩形窗w.（2)三角窗(Triangular Window)调用格式:w=triang（n） ，根据长度n 产生一个三角窗w.（3)汉宁窗(Hanning Window）调用格式：w=hanning(n) ，根据长度n 产生一个汉宁窗w.(4）海明窗（Hamming Window）调用格式：w=hamming（n) ,根据长度n 产生一个海明窗w.（5)布拉克曼窗(Blackman Window）调用格式:w=blackman(n) ，根据长度n 产生一个布拉克曼窗w。（6)恺撒窗（Kaiser Window）调用格式：w=kaiser（n，beta) ，根据长度n 和影响窗函数旁瓣的参数产生一个恺撒窗w。2、基于窗函数的FIR 滤波器设计利用MATLAB 提供的函数firl 来实现调用格式:firl （n,Wn,ftype,Window),n 为阶数、Wn 是截止频率（如果输入是形如W1 W2的矢量时，本函数将设计带通滤波器，其通带为W1W2）、ftype 是滤波器的类型(低通-省略该参数、高通-ftype=high、带阻-ftype=stop）、Window 是窗函数.例6.1： 设计一个长度为8 的线性相位FIR 滤波器.其理想幅频特性满足Window=boxcar（8)；b=fir1(7，0.4，Window）；freqz(b,1）例6。2:设计线性相位带通滤波器，其长度N=15，上下边带截止频率分别为W1= 0。3，w2=0。5Window=blackman（16）；b=fir1（15，0。3 0。5，Window）；freqz(b，1） 例6.3：MATLAB中的chirp。mat文件中存储信号的数据，该信号的大部分号能量集中在Fs/4（或二分之一奈奎斯特）以上，试设计一个34阶的FIR高通滤波器,滤除频率低于Fs/4的信号成分,其中滤波器的截止频率为0.48，阻带衰减为30dB，滤波器窗采用切比雪夫窗clear all；load chirpwindow=chebwin(35，30）；b=fir1（34，0.48,high，window）；yfit=filter(b,1,y）;Py，fy=pburg(y,10,512,Fs）;Pyfit，fyfit=pburg（yfit,10,512，Fs);plot（fy,10log10（Py）,.，fyfit， 10*log10(Pyfit);grid onylabel(幅度（dB)xlabel（频率(Hz）legend(滤波前的线性调频信号， 滤波后的线性调频信号）