数字信号处理MATLAB仿真设计.docx

数字信号处理MATLAB仿真设计.docx

《数字信号处理MATLAB仿真设计.docx》由会员分享，可在线阅读，更多相关《数字信号处理MATLAB仿真设计.docx（32页珍藏版）》请在冰豆网上搜索。

数字信号处理MATLAB仿真设计

实验一数字信号处理的Matlab仿真

一、实验目的

1、掌握连续信号及其MATLAB实现方法；

2、掌握离散信号及其MATLAB实现方法

3、掌握离散信号的基本运算方法，以及MATLAB实现

4、了解离散傅里叶变换的MATLAB实现

5、了解IIR数字滤波器设计

6、了解FIR数字滤波器设计1

二、实验设备

计算机，Matlab软件

三、实验容

（一）、连续信号及其MATLAB实现

1、单位冲击信号

例1.1：

t=1/A=50时，单位脉冲序列的MATLAB实现程序如下：

clearall;

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=[t1t2t3];

u=[u1u2u3];

plot（t,u）

axis（[-0.510A+2]）

2、任意函数

例1.2：

用MATLAB画出如下表达式的脉冲序列

clearall;

t=-2:

1:

3;

N=length（t）;

x=zeros（1,N）;

x

（1）=0.4;

x

（2）=0.8

x（3）=1.2;

x（4）=1.5;

x（5）=1.0;

x（6）=0.7;

stem（t,x）;

axis（[-2.23.201.7]）

3、单位阶跃函数

例1.3：

用MATLAB实现单位阶跃函数

clearall;

t=-0.5:

0.001:

1;

t0=0;

u=stepfun（t,t0）;

plot（t,u）

axis（[-0.51-0.21.2]）

4、斜坡函数

例1.4：

用MATLAB实现g（t）=3（t-1）

clearall;

t=0:

0.01:

3;

B=3;

t0=1;

u=stepfun（t,t0）;

n=length（t）;

fori=1:

n

u（i）=B*u（i）*（t（i）-t0）;

end

plot（t,u）

axis（[-0.23.1-0.26.2]）

5、实指数函数

例1.5：

用MATLAB实现

clearall;

t=0:

0.001:

3;

A=3;

a=0.5;

u=A*exp（a*t）;

plot（t,u）

axis（[-0.23.1-0.214]）

6、正弦函数

例1.6：

用MATLAB实现正弦函数f（t）=3cos（10πt+1）

clearall;

t=-0.5:

0.001:

1;

A=3;

f=5;

fai=1;

u=A*sin（2*pi*f*t+fai）;

plot（t,u）

axis（[-0.51-3.23.2]）

（二）、离散信号及其MATLAB实现

1、单位冲激序列

例2.1：

用MATLAB产生64点的单位冲激序列

clearall;

N=64;

x=zeros（1,N）;

x

（1）=1;

xn=0:

N-1;

stem（xn,x）

axis（[-16501.1]）

2、任意序列

例2.2：

用MATLAB画出如下表达式的脉冲序列

clearall;

N=8;

x=zeros（1,N）;

x

（1）=8.0;

x

（2）=3.4

x（3）=1.8;

x（4）=5.6;

x（5）=2.9;

x（6）=0.7;

xn=0:

N-1;

stem（xn,x）

axis（[-1808.2]）

3、单位阶跃序列

例2.3：

用MATLAB实现单位阶跃函数

clearall;

N=32;

x=ones（1,N）;

xn=0:

N-1;

stem（xn,x）

axis（[-13201.1]）

4、斜坡序列

例2.4：

用MATLAB实现g（n）=3（n-4）点数为32的斜坡序列

clearall;

N=32;

k=4

B=3;

t0=1;

x=[zeros（1,k）ones（1,N-k）];

fori=1:

N

x（i）=B*x（i）*（i-k）;

end

xn=0:

N-1;

stem（xn,x）

axis（[-132090]）

5、正弦序列

例2.5：

用MATLAB实现幅度A=3，频率f=100，初始相位Φ=1.2，点数为32的正弦信号

clearall;

N=32;

A=3;

f=100;

fai=1.2;

xn=0:

N-1;

x=A*sin（2*pi*f*（xn/N）+fai）;

stem（xn,x）

axis（[-132-3.23.2]）

6、实指数序列

例2.6：

用MATLAB实现

，点数为32的实指数序列

clearall;

N=32;

A=3;

a=0.7;

xn=0:

N-1;

x=A*a.^xn;

stem（xn,x）

7、复指数序列

例2.7：

用MATLAB实现幅度A=3，a=0.7，角频率ω=314，点数为32的实指数序列

clearall;

N=32;

A=3;

a=0.7;

w=314;

xn=0:

N-1;

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的均匀分布的随机序列与高斯随机序列

clearall;

N=32;

x_rand=rand（1,N）;

x_randn=randn（1,N）;

xn=0:

N-1;

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程序如下：

clearall;

N=32;

w=100;

k=3;

x1=zeros（1,k）;

xn=0:

N-1;

x2=sin（100*xn）;

figure

（1）

stem（xn,x2）

x=[x1x2];

axis（[-1N-1.11.1]）

N=N+k;

xn=0:

N-1;

figure

（2）

stem（xn,x）

axis（[-1N-1.11.1]）

2、信号相加

若信号

，值得注意的是当序列

和

的长度不相等或者位置不对应时，首先应该使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加

例3.2：

用MATLAB实现两序列相加

clearall;

n1=0:

3

x1=[20.50.91];

figure

（1）

stem（n1,x1）

axis（[-1802.1]）

n2=0:

7

x2=[00.10.20.30.40.50.60.7];

figure

（2）

stem（n2,x2）

axis（[-1800.8]）

n=0:

7;

x1=[x1zeros（1,8-length（n1））];

x2=[zeros（1,8-length（n2））,x2];

x=x1+x2;

figure（3）

stem（n,x）

axis（[-1802.1]）

3、信号相乘

信号序列

和

相乘所得信号

的表达式为：

这是样本与样本之间的点乘运算，在MATLAB中可采用“.*”来实现，但是在信号序列相乘之前，应对其做与相加运算一样的操作。

例3.3：

用MATLAB实现上例中两序列相乘

clearall;

n1=0:

3

x1=[20.50.91];

figure

（1）

stem（n1,x1）

axis（[-1802.1]）

n2=0:

7

x2=[00.10.20.30.40.50.60.7];

figure

（2）

stem（n2,x2）

axis（[-1800.8]）

n=0:

7;

x1=[x1zeros（1,8-length（n1））];

x2=[zeros（1,8-length（n2））,x2];

x=x1.*x2;

figure（3）

stem（n,x）

axis（[-1800.35]）

4、信号翻转

信号翻转的表达式为：

y（n）=x（-n），在MATLAB中可以用fliplr函数实现此操作

例3.4：

用MATLAB实现“信号相加”中的

序列翻转

clearall;

n=0:

3

x1=[20.50.91];

x=fliplr（x1）;

stem（n,x）

axis（[-1402.1]）

5、信号和

对于N点信号

，其和的定义为：

例3.5：

用MATLAB实现“信号相加”中的

序列和

clearall;

n=0:

3

x1=[20.50.91];

x=sum（x1）

6、信号积

对于N点信号

，其积的定义为：

例3.5：

用MATLAB实现“信号相加”中的

序列积

clearall;

n=0:

3

x1=[20.50.91];

x=prod（x1）

（四）、离散傅里叶变换的MATLAB实现

例４：

若

是一个N=32的有限序列，利用MATLAB计算它的DFT并画出图形。

N=32;

n=0:

N-1;

xn=cos（pi*n/6）;

k=0:

N-1;

WN=exp（-j*2*pi/N）;

nk=n’*k;

WNnk=WN.^nk;

Xk=xn*WNnk;

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，绘出该滤波器的幅频于相频特性，以及其冲击响应图

clearall;

N=10;

Wn=[100200]/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=1000Hz

clearall;

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：

设计一个切比雪夫Ⅱ型数字带通滤波器，要求带通围100-250Hz，带阻上限为300Hz，下限为50Hz，通带纹波小于3dB，阻带纹波为30dB，抽样频率为1000Hz，并利用最小的阶次实现。

clearall;

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）矩形窗（RectangleWindow）

调用格式：

w=boxcar（n），根据长度n产生一个矩形窗w。

（2）三角窗（TriangularWindow）

调用格式：

w=triang（n），根据长度n产生一个三角窗w。

（3）汉宁窗（HanningWindow）

调用格式：

w=hanning（n），根据长度n产生一个汉宁窗w。

（4）海明窗（HammingWindow）

调用格式：

w=hamming（n），根据长度n产生一个海明窗w。

（5）布拉克曼窗（BlackmanWindow）

调用格式：

w=blackman（n），根据长度n产生一个布拉克曼窗w。

（6）恺撒窗（KaiserWindow）

调用格式：

w=kaiser（n,beta），根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w。

2、基于窗函数的FIR滤波器设计

利用MATLAB提供的函数firl来实现

调用格式：

firl（n,Wn,’ftype’,Window），n为阶数、Wn是截止频率（如果输入是形如[W1W2]

的矢量时，本函数将设计带通滤波器，其通带为W1<ω

通-省略该参数、高通-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.5

π

Window=blackman（16）;

b=fir1（15,[0.30.5],Window）;

freqz（b,1）

例6.3：

MATLAB中的chirp.mat文件中存储信号ｙ的数据，该信号的大部分号能量集中在Fs/4（或二分之一奈奎斯特）以上，试设计一个34阶的FIR高通滤波器，滤除频率低于Fs/4的信号成分，其中滤波器的截止频率为0.48，阻带衰减为30dB，滤波器窗采用切比雪夫窗

clearall;

loadchirp

window=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,10*log10（Py）,’.’,fyfit,10*log10（Pyfit））;

gridon

ylabel（‘幅度（dB）’）

xlabel（‘频率（Hz’）

legend（‘滤波前的线性调频信号’,‘滤波后的线性调频信号’）

实验一：

MATLAB语言上机操作实践

1.简单的数组赋值方法

（1）

a=[123;456;789]%给定数组a

a（4,2）=11%给数组a的第四行第二列元素赋11

a（5,:

）=[-13-14-15]%给数组a的第五行元素赋值

a（4,3）=abs（a（5,1））%数组a的第四行第三列元素为第五行第一列元素的模

a（[2,5],:

）=[]%数组a的第二行和第五行取空

a/2%对数组a各元素对应除以2

a（4,:

）=[sqrt（3）（4+5）/6*2-7]%给数组a的第四行元素赋值

（2）

B=[1+2i,3+4i;5+6i,7+8i]

C=[1,3;5,7]+[2,4;6,8]*i

C=[1,3;5,7]+[2,4;6,8]i

D=sqrt（2+3i）%对复数2+3i求根号

D*D%对D求平方

E=C'%对C转置

F=conj（C）

G=conj（C）'

（3）

H1=ones（3,2）

H2=zeros（2,3）

H3=eye（4）

2.数组的基本运算

（1）

A=[135]

B=[246]

C=A+B

D=A-2

E=B-A

（2）

F1=A*3

F2=A.*B

F3=A./B

F4=A.\B

F5=B.\A

F6=B.^A

F7=2./B

F8=B.\2

（3）

Z1=A*B'

Z2=B'*A

3.常用函数及相应的信号波形显示

例1-1：

①

t=0:

0.05:

3;%给建立时间数组

f=2*sin（2*pi*t）;%生成函数f

plot（t,f）;%用plot作连续信号的曲线

title（'f（t）-t曲线'）;%在图上端标注图名

xlabel（'t'）;%标注横坐标

ylabel（'f（t）'）;%标注纵坐标

②

t=0:

0.05:

3;%给建立时间数组

f=2*sin（2*pi*t）;%生成函数f

subplot（2,2,1）,plot（t,f）;%建立2*2子图轴系，在图1处绘线性图

title（'plot（t,f）'）;

subplot（2,2,2）,stem（t,f）;%建立2*2子图轴系，在图2处绘脉冲图

title（'stem（t,f）'）;

subplot（2,2,3）,stairs（t,f）;%建立2*2子图轴系，在图3处绘阶梯图

title（'stairs（t,f）'）;

subplot（2,2,4）,bar（t,f）;%建立2*2子图轴系，在图4处绘条形图

title（'bar（t,f）'）;

练习题：

（1）

t=0:

0.05:

4;%给建立时间数组

f=4*exp（（-2）*t）;%生成函数f

subplot（2,2,1）,plot（t,f）;%建立2*2子图轴系，在图1处绘线性图

title（'plot（t,f）'）;

subplot（2,2,2）,stem（t,f）;%建立2*2子图轴系，在图2处绘脉冲图

title（'stem（t,f）'）;

subplot（2,2,3）,stairs（t,f）;%建立2*2子图轴系，在图3处绘阶梯图

title（'stairs（t,f）'）;

subplot（2,2,4）,bar（t,f）;%建立2*2子图轴系，在图4处绘条形图

title（'bar（t,f）'）;

（2）

t=0:

0.05:

3;%给建立时间数组

f=exp（（-1）*t）.*cos（2*pi*t）;%生成函数f

subplot（2,2,1）,plot（t,f）;%建立2*2子图轴系，在图1处绘线性图

title（'plot（t,f）'）;

subplot（2,2,2）,stem（t,f）;%建立2*2子图轴系，在图2处绘脉冲图

title（'stem（t,f）'）;

subplot（2,2,3）,stairs（t,f）;%建立2*2子图轴系，在图3处绘阶梯图

title（'stairs（t,f）'）;

subplot（2,2,4）,bar（t,f）;%建立2*2子图轴系，在图4处绘条形图

title（'bar（t,f）'）;

（3）

k=0:

0.05:

10;%给建立时间数组

f=k;%生成函数f

subplot（2,2,1）,plot（k,f）;%建立2*2子图轴系，在图1处绘线性图

title（'plot（k,f）'）;

subplot（2,2,2）,stem（k,f）;%建立2*2子图轴系，在图2处绘脉冲图

title（'stem（k,f）'）;

subplot（2,2,3）,stairs（k,f）;%建立2*2子图轴系，在图3处绘阶梯图

title（'stairs（k,f）'）;

subplot（2,2,4）,bar（k,f）;%建立2*2子图轴系，在图4处绘条形图

title（'bar（k,f）'）;

（4）

k=-20:

0.05:

20;

f=k.*sin（k）;

subplot（2,2,1）,plot（k,f）;%建立2*2子图轴系，在图1处绘线性图

title（'plot（k,f）'）;

subplot（2,2,2）,stem（k,f）;%建立2*2子图轴系，在图2处绘脉冲图

title（'stem（k,f）'）;

subplot（2,2,3）,stairs（k,f）;%建立2*2子图轴系，在图3处绘阶梯图

title（'stairs（k,f）'）;

subplot（2,2,4）,bar（k,f）;%建立2*2子图轴系，在图4处绘条形图

title（'bar（k,f）'）;

4.简单的流程控制编程

例1-2

X=0;%给变量X赋初值0

forn=1:

32%确定循环变量及其取值围

X=X+n^2;%循环体

end

练习题：

（1）

X=0;%给变量X赋初值0

forn=1:

20%确定循环变量及其取值围

X=X+（2*n-1）^2;%循环体

end

（2）

X=0;%给变量X赋初值0

forn=1:

99%确定循环变量及其取值围

X=X+n*（n+1）;

end

（3）

a=[];%定义数组a为空

a

（1）=1;%给数组第一个元素赋值1

a

（2）=1;%给数组第二个元素赋值1

fork=1:

18%确定循环变量及其取值围

a（k+2）=a（k）+a（k+1）;%循环体

end

实验二：

时域离散信号的产生

（2）

①

n1=-3;n2=4;n0=0;%在起点为n1、终点为n2的围，于n0处产生冲激

n=n1:

n2;%生成离散信号的时间序列

f=[n==n0];%建立信号f

stem（n,f,'filled'）;%绘制脉冲杆图，且在圆心处用实心圆表示

title（'单位脉冲序列'）

axis（[n1n201.1*max（f）]）;%限定横坐标和纵坐标的显示围

ylabel（'f（n）'）;

②

n1=-5;n2=5;n0=0;

n=n1:

n2;%生成离散信号的时间序列

f=[n>=n0];%生成离散序列f

stem（n,f,'filled'）;

title（'单位阶跃序列'）

axis（[n1n201.1*max（f）]）;%限定横坐标和纵坐标的显示围

ylabel（'幅度f（n）'）;

③

n=0:

16;

f=exp（（0.1+1.6*pi*j）*n）;%建立信号f

stem（n,f,'fill'）,gridon

xlabel（'n'）,title（'f（n）'）;

④

n=0:

20;

f=3*sin（pi/4*n）;%建立信号f

stem（n,f,'fill'）,xlabel（'n'）,gridon

title（'f（n）'）

axis（[0,20,-3.5,3.5]）;%限定横坐标和纵坐标的显示围

⑤

n=-20:

20;

f=sinc（n/5）;%建立信号f

stem（n