matlab 基础及数字信号处理实验.docx
《matlab 基础及数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《matlab 基础及数字信号处理实验.docx（26页珍藏版）》请在冰点文库上搜索。
matlab基础及数字信号处理实验
MATLAB基础
一、MATLAB的工作环境
1．命令窗口：
是用户和MATLAB系统交互的主要窗口。
在该窗口中，用户可以运行函数，执行MATLAB的基本操作命令以及对MATLAB系统的参数设置等操作。
在命令提示符>>后输入命令
如：
t=[1,2,3;4,5,6;7,8,9];
完成对t的赋值
2．帮助的使用：
直接在命令行输入help
在菜单栏选择help—matlabhelp
直接按F1
3．图形窗（Figure）：
用于显示绘出的图形。
通常只要执行了任意一种绘图命令，图形窗都会自动产生。
绘图都是在这个图形窗中进行。
如果要在建一个图形窗，则可在命令窗中输入figure命令，MATALB就回新建一个图形窗口，并自动给它排出序号。
4．文本编辑窗：
其作用是用来创建、编辑和调试MATLAB的相关文件（或称程序，即.M文件），它与一般的编辑调试器有相似的功能。
通常，MATLAB命令编辑有行命令方式和文件两种。
行命令方式，即在命令窗口中一行一行的输入命令，计算机对每一行命令做出反应。
这种方式，只能用于编辑简单的程序，当程序比较复杂的时候，把程序写成一个由多行语句组成的文件（.M文件），让MATLAB来执行这个程序中的全部语句，MATLAB文本编辑器的功能就是完成编写，修改和调试这种程序。
其进入方式：
1）file—new/open—M-file
2）直接在命令窗口输入：
edit/edit文件名
二、MATLAB的基本语法
1．变量及其赋值
1）赋值要求
在MATLAB中，变量和常量的标识符最长允许19个字符。
MATLAB内部只有一种数据格式，就是双精度类型，对应于64位二进制。
赋值就是把数赋予代表常量或变量的标识符。
基本格式为：
变量=表达式（或数）
在MATLAB中，变量都代表矩阵。
列矢量可被当做只有一行的矩阵，行矢量也可被当做只有一个行的矩阵，标量（或常数）应该做是1×1的矩阵。
输入矩阵时，应遵循一下规则：
整个矩阵的值都应放在方括号中；同一行的个元素以逗号或空格分开；不同行的元素以分号隔开。
例如：
t=[1,2,3;4,5,6;7,8,9]
若不希望有处理结果，则在语句的后面加上分号
例如：
t=[1,2,3;4,5,6;7,8,9];
若想看t的值，接着输入>>t即可
2）变量元素的赋值
单独给矩阵的某个元素（用圆括号表示）赋值，如：
>>t（2,3）=10;t（1,2）=1.5;
如果赋值元素的下标超出了原有矩阵的大小，矩阵行列会自动扩展，如
>>t（4,2）=11;
变量的阶数可以用命令size来获取：
>>size（t）
ans为MATLAB自动给出的一个临时变量。
3）赋值技巧
在MATLAB中，为变量的赋值提供了一些简便快捷的方法
1利用冒号，给全行元素赋值
>>t（4,:
）=[5,3,2]
2利用特殊矩阵和数组赋值
zeros：
生成一个元素全部为0的矩阵或数组
ones：
生成一个元素全部为1的矩阵或数组
rand：
生成随机矩阵或数组，元素是在（0,1）之间服从均匀分布
eye：
生成一个单位矩阵或数组
linspace：
生成一个线性间隔的行矢量
randn：
生成随机矩阵或数组，元素服从均值为0，方差为1的正态分布
具体函数的用法请查help
4）复数的赋值方式
复数的虚部部分用i或j表示，这是MATLAB启动时自动设定的：
>>c=3+2.5i
>>z=[1+2j,3+4j;5+6j,7+8j]
>>z=[1+2*j,3+4*j;5+6*j,7+8*j]
5）MATLAB内部特殊变量和常数
①变量ans：
临时变量，通常表示当前的答案
3常数eps：
表示浮点相对精度，按IEEE标准，eps=2-52近似为2.2204e-016
4常数pi：
表示圆周率
5常数Inf：
表示正无穷大
6虚数单位i,j：
表示复数的虚部单位
7NaN：
表示非数值，如Inf-Inf
2．运算符、复数运算及流程控制
算术运算、关系运算、逻辑操作（对复数都有效）
1）>>f=sqrt（1+2i）
>>f*f
2）复数的共轭可利用函数conj来完成
>>x=conj（f）
3.基本的数学函数：
>>helpelfun查看该库函数
4．基本绘图命令与画图方法
1）基本绘图命令
二维图形（graph2d子目录）和三维图形（graphics子目录），特殊绘图命令在specgraph子目录。
1figure：
新建一个绘图窗口，并自动给它排序号
2subplot：
将窗口分成n×m个子图形窗口，并选择第p个子图形窗口作为当前窗口，供绘图函数使用。
格式：
subplot（m,n,p）或subplot（nmp）
3二维曲线绘制函数plot：
格式一：
plot（y）——输入一个数组的情况
功能：
如果y是一个数组，plot（y）给出线性直角坐标的二维图。
以y中元素的下标作为x坐标，y中元素的值作为y坐标，并将各点连线
>>y=5*（rand（1,10）-0.5）;
>>plot（y）
---------------------------------------------
>>xlable（'n'）;ylable（'Y'）;%标注x轴,y轴
>>xlabel（'n'）;ylabel（'Y'）;%标注x轴,y轴
>>plot（y）
>>gridon;%设置网格
格式二：
plot（x,y）——输入2个数组的情况
功能：
绘出以x元素为横坐标，y元素为纵坐标的曲线，x和y必须具有相同的长度
closeall;%关闭所有的图形窗口
>>clear;%清除工作空间中的所有变量
>>t=0:
0.02:
2;%利用冒号，对t从[0,2]赋值，步长为0.02
>>y=3*sin（t）.*exp（-2*t）;%对y赋值，.*为矩阵元素相乘
>>ploy（t,y）
格式三：
plot（y,’字符串’）或plot（x,y,’字符串’）
用户只需要在plot变量后加一个单引号，单引号内部放入线型、颜色以及标记符号的标识符即可完成。
plot（t,y,'*r-'）
颜色标识符
符号标识符
线型标识符
b
蓝
.
点
-
实线
g
绿
o
圆圈
：
点线
r
红
x
叉
-.
点划线
c
青
+
+号
--
虚线
m
品红
*
星号
y
黄
s
平方号
k
黑
d
钻石符号
例：
利用subplot命令将图形窗口划分为2×1的子窗口，并在不同的子窗口用不同的颜色线型和标记符号绘制图形：
t=0:
0.5:
4*pi;%对t变量进行赋值
y1=exp（-0.1*t）.*sin（t）;
y2=exp（-0.1*t）.*sin（t+1）;
figure
（1）;
subplot（2,1,1）;
plot（t,y1,'-bx'）;
title（'y1=exp（-0.1t）*sin（t）'）;
xlabel（'t'）;ylabel（'y1'）;
subplot（2,1,2）;
plot（t,y2,'--k*'）;
title（'y2=exp（-0.1t）*sin（t+1）'）;
xlabel（'t'）;ylabel（'y2'）;
格式四：
plot（x1,y1,x2,y2,…,xn,yn）——输入多对数组的情况
功能：
x1,y1,x2,y2,…,xn,yn分别为向量对，每一对x-y可以绘出一条图线，这样就可以在一张图上画出多图线，每一组向量的长度可以不同，在其后都可以加线型标记符号。
t1=0:
0.5:
4*pi;y1=exp（-0.1*t1）.*sin（t1）;
t2=0:
0.2:
2*pi;y2=exp（-0.5*t2）.*sin（5*t2+1）;
plot（t1,y1,'+k-',t2,y2,':
b'）;
xlabel（'时间'）;ylabel（'Y'）;
legend（'\函数y1','\函数y2'）;%标注图例
三、MATLAB在数字信号处理中的应用
1.序列的图形表示stem
格式一：
stem（x,’string’）
功能：
绘制序列x（n）的棒状图，横坐标为序列的下标序号，棒的末端用string指定的标记符号表示，如前表所示。
格式二：
stem（x,y,’string’）
功能：
绘制序列的棒状图，横坐标由矢量x（n）指定。
例：
绘制离散时间信号x[-1]=-1，x[0]=1,x[1]=2,x[2]=1,x[3]=0,x[4]=-1，其他时间x[n]=0的棒状图。
x=[0,0,-1,1,2,1,-1,0,0];
stem（x）;
line（[0,9],[0,0]）;%绘制横轴线
xlabel（'n'）;ylabel（'x[n]'）;
gridon;
出现问题：
横坐标不对，向量下标从1开始，不能取0或者负值，用格式二，用另一个向量来表示横坐标
n=-3:
5;%定位时间变量
x=[0,0,-1,1,2,1,-1,0,0];
stem（n,x,'.'）;
line（[-3,5],[0,0]）;%绘制横轴线
xlabel（'n'）;ylabel（'x[n]'）;
gridon;
2．信号的基本运算
1）序列的反褶与序列的累加
2）两序列的卷积
例：
计算下列卷积，并图示各序列及其卷积结果
x（n）=0.9nR20（n）,h（n）=R10（n）
其MATLAB程序如下：
xn=0:
19;x=（0.9）.^xn;%产生x（n）
hn=0:
9;h=ones（1,10）;%利用ones产生h（n）
y=conv（x,h）;%计算y（n）=x（n）*h（n）
yn=0:
（20+10-2）;%计算卷积序列y（n）的时间序列
subplot（3,1,1）,stem（xn,x,'.'）,xlabel（'n'）,ylabel（'x（n）'）;
subplot（3,1,2）,stem（hn,h,'.'）,xlabel（'n'）,ylabel（'h（n）'）;
subplot（3,1,3）,stem（yn,y,'.'）,xlabel（'n'）,ylabel（'y（n）'）;
3）两序列的相关运算：
4）离散信号的能量和功率
matlab实现为：
或
matlab实现为：
或
3．LTI系统描述
1）LTI系统的描述模型
1常系数线性差分方程
是n时刻的输出，
是n时刻的输入，
和
均为常数，N为系统的阶数。
2系统的传递函数（系统函数）
在matlab中，系统函数由分子，分母两个多项式的系数来表示，系数为降幂排列，若某阶次项没有，则需用0补上。
例如：
可表示为：
分子：
num=[1,0.2,3]或A=[1,0.2,3]
分母：
den=[1,0.5,0,2]或B=[1,0.5,0,2]
3零极点增益模型
在matlab中，增益系数、零点向量、极点向量分别用k,z,p表示，向量为列向量
4二次分式模型
是零极点增益模型的变形
5状态空间模型
设x为状态变量，u为输入，离散LTI系统的状态方程为
在matlab中用矩阵（或向量）ABCD表示系统的状态空间模型
2）模型函数之间的转换
在matlab中，分别用ss，sos，tf，zp符号表示状态空间模型、二次分式模型、传递函数模型和零极点增益模型，并同时提供了一组不同线性系统之间转换的函数，如表所示：
函数名
功能说明
tf2zp
传递函数模型转换为零极点增益模型
tf2sos
传递函数模型转换为二次分式模型
zp2tf
零极点增益模型转换为传递函数模型
zp2sos
零极点增益模型转换为二次分式模型
sos2tf
二次分式模型转换为传递函数模型
sos2zp
二次分式模型转换为零极点增益模型
ss2tf
状态空间模型转换为传递函数模型
ss2zp
状态空间模型转换为零极点增益模型
ss2sos
状态空间模型转换为二次分式模型
zp2ss
零极点增益模型转换为状态空间模型
sos2ss
二次分式模型转换为状态空间模型
例如：
tf2zp函数的基本格式为
[z,p,k]=tf2zp（num,den）
功能：
求系统传递函数的零点向量z、极点向量p和增益系数k。
用于离散系统时分子和分母多项式的长度必须相同，否则补零。
例：
求离散时间系统
的零极点向量和增益系数
num=[2,3];
den=[1,0.4,1];
[num,den]=eqtflength（num,den）;%使长度相等
[z,p,k]=tf2zp（num,den）
--------------------------------
z=0
-1.5000
p=-0.2000+0.9798i
-0.2000-0.9798i
k=2
表示零点为0，-1.5；极点共轭对；增益系数2
4．离散LTI系统的时域响应
1）已知单位冲激响应时系统对任意输入的响应
2）已知系统函数时对任意输入的响应
格式：
功能：
计算出传递函数
形式下系统的响应，x为输入信号，当dlsim函数不带输出变量时，可在当前图形窗口中直接绘出系统的输出响应曲线；当带有输出变量调用函数时，可得到系统输出响应曲线的数据，而不直接绘出曲线。
例：
有二阶系统
，求系统对100点随机噪声的响应曲线。
num=[2,-3.4,5.5];
den=[1,-1.2,0.8];
u=randn（1,100）;%产生随机数
dlsim（num,den,u）;
title（'随机噪声响应'）;
3）已知系统函数求LTI系统的冲激响应
利用系统函数
求单位冲激响应
，matlab提供了一条专用的函数命令：
格式一：
h=dimpulse（num,den）
功能：
返回多项式传递函数
表示的冲激响应，其中num和den按z的降幂排列的多项式系数。
格式二：
dimpulse（num,den,n）
功能：
可利用用户指定的样点n来求出系统的单位冲激响应。
注意：
当dimpulse函数不带输出变量时，可在当前图形窗口中直接绘出系统的输出响应曲线；当带有输出变量调用函数时，可得到系统输出响应曲线的数据，而不直接绘出曲线。
例：
有二阶系统
，求系统的单位冲激响应。
num=[2,-3.5,1.5];
den=[1,-1.7,0.3];
dimpulse（num,den）;
title（'离散LTI系统的冲激响应'）;
5．线性时不变系统的频率响应
1）求数字滤波器
的频率响应函数freqz
格式一：
H=freqz（B,A,W）
功能：
计算由向量W（rad）指定的数字频率点上（通常在[0,π]范围的频率）数字滤波器
的频率响应
，结果存于H向量中。
向量B和A分别为
的分子和分母多项式系数。
格式二：
[H,W]=freqz（B,A,N）
功能：
计算出N个频率点上的频率响应存放在H向量中，N个频率存放在向量W中。
freqz函数自动将这N个频率点均匀设置在频率范围[0,π]上。
默认W和N时，freqz函数将自动选取512个频率点计算。
格式三：
[H,W]=freqz（B,A,N,’whole’）
功能：
计算出N个频率点上的频率响应存放在H向量中，N个频率存放在向量W中。
freqz函数自动将这N个频率点均匀设置在频率范围[0,2π]上。
格式四：
[H,F]=freqz（B,A,N,Fs）或[H,F]=freqz（B,A,N,Fs,’whole’）
功能：
按指的抽样频率Fs（Hz）以及N个频率点，返回频率响应向量H和频率向量F（Hz）
格式五：
H=freqz（B,A,F,Fs）
功能：
计算出指定频率向量F（Hz）和抽样频率Fs（Hz）的频率响应
格式六：
[H,W,S]=freqz（…）或[H,F,S]=freqz（…）
功能：
返回有关函数freqzplot的相关信息。
说明：
不带输出向量时freqz函数将自动绘出频率响应的幅频和相频曲线
例：
已知某滤波器的系统函数为
，求该滤波器的频率响应
B=[1,0,0,0,0,0,0,0,-1];
A=1;
freqz（B,A）;
2）滤波函数filter
格式：
y=filter（B,A,x）
功能：
对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。
B和A分别为数字滤波器系统函数
的分子和分母多项式系数。
要求a
（1）=1，否则就应归一化。
例：
设系统的差分方程为
，求该系统对信号
的响应。
B=1;
A=[1,-0.8];
n=0:
31;x=0.8.^n;
y=filter（B,A,x）;
subplot（2,1,1）;stem（x）
subplot（2,1,2）;stem（y）
6．快速傅里叶变换FFT算法
1）一维快速傅里叶正变换函数fft
格式：
X=fft（x,N）
功能：
采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，按x的长度计算，当N为2的整数次幂时，按基-2算法，否则用混合算法。
2）一维快速傅里叶反变换函数ifft
格式：
X=ifft（x,N）
功能：
采用FFT算法计算序列向量x的N点IDFT变换
3）线性调频z变换czt
格式：
y=czt（x,m,w,a）
功能：
计算由z=a*w.^（0:
m-1）定义的z平面螺旋线上各点的z变换。
a规定了起点，w规定了相邻点的比例，m规定了变换长度。
后三个变元的默认值是a=1,w=exp（j*2*pi/m）以及m=length（x）。
因此y=czt（x）就等于y=fft（x）。
4）基于FFT重叠相加法FIR滤波器实现函数fftfilt
格式：
y=fftfilt（h,x）
功能：
采用重叠相加法FFT实现对信号向量x的快速滤波，得到输出序列向量y，向量h为FIR滤波器的单位脉冲响应。
例：
用快速傅里叶变换FFT计算下面两个序列的卷积
并测试直接卷积计算和快速卷积的时间。
xn=sin（0.4*[1:
15]）;
hn=0.9.^（1:
20）;
tic,%计时开始
y=conv（xn,hn）;
toc,%计时结束
M=length（xn）;N=length（hn）;
nx=1:
M;nh=1:
N;
L=pow2（nextpow2（M+N-1））;%取L为大于等于且最接近M+N-1的2的正次幂
tic,%计时开始
Xk=fft（xn,L）;
Hk=fft（hn,L）;
Yk=Xk.*Hk;
yn=ifft（Yk,L）;
toc,%计时结束
subplot（2,2,1）,stem（nx,xn,'.'）;ylabel（'x（n）'）;
subplot（2,2,2）,stem（nh,hn,'.'）;ylabel（'h（n）'）;
subplot（2,1,2）;ny=1:
L;stem（ny,real（yn）,'.'）;ylabel（'y（n）'）;
Elapsedtimeis0.003000seconds.
Elapsedtimeis0.000000seconds.
7．IIR数字滤波器的设计方法
1）冲激响应不变法impinvar
格式：
[BZ,AZ]=impinvar（B,A,Fs）
功能：
把具有[B,A]模拟滤波器传递函数模型转换为采样频率为Fs的数字滤波器的传递函数模型[BZ,AZ]，Fs默认值为1。
例：
一个4阶的Butterworth模拟低通滤波器的系统函数如下：
试用冲激响应不变法求出Butterworth模拟低通数字滤波器的系统函数。
num=1;
den=[1,sqrt（5）,2,sqrt
（2）,1];
[num1,den1]=impinvar（num,den）
---------------------------------------------------
num1=-0.00000.09420.21580.0311
den1=1.0000-2.00321.9982-0.76120.1069
2）双线性变换法bilinear
格式一：
[Zd,Pd,Kd]=bilinear（Z,P,K,Fs）
功能：
把模拟滤波器的零极点模型转换成数字滤波器的零极点模型，Fs是采样频率
格式二：
[numd,dend]=bilinear（num.den,Fs）
功能：
把模拟滤波器的传递函数模型转换为数字滤波器的传递函数模型。
例：
一个三阶的模拟Butterworth模拟低通滤波器的系统函数如下：
，试用双线性变换法求出数字Butterworth数字低通滤波器的系统函数。
num=1;
den=[1,sqrt（3）,sqrt
（2）,1];
[num1,den1]=bilinear（num,den,1）
----------------------------------------------
num1=0.05330.15990.15990.0533
den1=1.0000-1.33820.9193-0.1546
3）IIR数字滤波器的频率变换实现
步骤：
1按一定的规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标
2根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc
3利用最小阶数N产生模拟低通原型
4利用截止频率Wc把模拟低通滤波器原型转换为模拟低通、高通、带通、带阻滤波器
5利用冲激响应不变法或双线性变换法把模拟滤波器转换为数字滤波器
表一IIR滤波器阶次估计
函数名
功能说明
buttord
计算Butterworth滤波器的阶次及截止频率
cleb1ord
计算ChebyshevⅠ滤波器的阶次
cheb2ord
计算ChebyshevⅡ滤波器的阶次
ellipord
计算椭圆滤波器的最小阶次
表二模拟低通滤波器原型设计
函数名
功能说明
buttap
Butterworth模拟低通滤波器原型设计
cleb1ap
ChebyshevⅠ模拟低通滤波器原型设计
cheb2ap
ChebyshevⅡ模拟低通滤波器原型设计
ellipap
椭圆模拟低通滤波器原型设计
表三模拟滤波器变换函数
函数名
功能说明
lp2bp
模拟低通转换为带通
lp2bs
模拟低通转换为带阻
lp2hp
模拟低通转换为高通
lp2
改变模拟低通的截止频率
例：
设计一个数字信号处理系统，它的采样频率为Fs=100Hz，希望在该系统中设计一个Butterworth型高通滤波器，使其通带中允许的最大衰减为0.5dB，阻带内最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。
wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;
[N,Wc]=buttord（wp,ws,rp,rs,'s'）;
[Z,P,K]=buttap（N）;%创建巴特沃斯低通滤波器模型
[A,B,C,D]=zp2ss（Z,P,K）;%零极点增益模型——>状态空间模型
[AT,BT,CT,DT]=lp2hp（A,B,C,D,Wc）;%低通——>高通
[num1,den1]=ss2tf（AT,BT,CT,
升级会员 