登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

我的博客

电子信息学习 资料共享

 
 
 

日志

 
 

MATLAB语言在DSP设计中的应用  

2010-01-24 21:40:16|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

 

 MATLAB语言在DSP设计中的应用

随着计算机技术的高速发展,计算机语言也得到了迅速发展。我们熟知的BASIC、FORTRAN、C等广泛地应用于各种场合。但从工程计算和图形显示的角度,这些语言并不方便。1984年,美国Mathworks公司正式推出了 MATLAB语言。MATLAB是“矩阵实验室”(MATrix LABoratoy)的缩写,是一种科学计算软件,主要适用于控制和信息处理领域的分析设计。它是一种以矩阵运算为基础的交互式程序语言,能够满足工程计算和绘图的需求。与其它计算机语言相比,其特点是简洁和智能化,适应科技专业人员的思维方式和书写习惯,使得编程和调试效率大大提高,并且很容易由用户自行扩展。因此,当前已成为美国和其它发达国家大学教学和科学研究中必不可少的工具。

MATLAB语言自1988年推出3.x(DOS)版本,目前已出了4.x,5.x,6.x等(Windows)版本。随着版本的升级,内容不断扩充。

§1 MATLAB的工作环境

MATLAB的工作环境主要由命令窗(Command Windows)、文本编辑器(File Editor)、若干个图形窗(Figure Windows)及文件管理器组成。MATLAB视窗采用了WINDOWS视窗风格(如图5-1-1)。各视窗之间的切换可用快捷键Alt+Tab。

图5-1-1  MATLAB的命令窗、文本编辑窗和图形窗

使用MATLAB4.x以上的版本,可在WINDOWS主界面上直接点击MATLAB图标,进入MATLAB命令窗口。在MATLAB命令窗下键入一条命令,按Enter键,该指令就被立即执行并显示结果。   

如果一个程序稍复杂一些,则需要采用文件方式,把程序写成一个由多条语句构成的文件。这时就需要用到文本编辑器。建立一个新文件,应在MATLAB命令窗口下点击空白文档符号或在File菜单下点击New,将打开MATLAB文本编辑器窗口,显示一个空白的文档。对已经存在的文件,则点击打开文件或在File菜单下点击Open,会自动进入文件选择窗口,找到文件后点亮并打开即可进入MATLAB文本编辑器窗口。在MATLAB文本编辑器窗口中建立的文件默认为.m文件。如果要建立的文件是M函数文件,即希望被其它程序象MATLAB中的库函数那样被调用,则文件的第一句应是函数申明行,如:

function  [y,w]=XYZ(x,t)

式中,function为MATLAB关键字,[]放置输出宗量,()中放置输入宗量,XYZ为函数名。当其它程序调用该函数时,只需在程序中直接使用function关键字后面的部分。函数申明行是M函数文件必不可少的一部分。

程序执行的结果以图形方式显示时,将自动打开图形窗。在程序中,图形窗命令为figure。MATLAB允许打开多个图形窗。如果程序中对图形窗没有编号,将按程序执行的顺序自动给图形窗编号。

在MATLAB命令窗下,还具有许多文件管理的功能。例如,我们自己编写的文件放在一个专门的文件夹中,则需要将这个文件夹的路径存盘。否则,这个文件夹中的文件将不能在MATLAB环境下执行。在MATLAB命令窗口File下选set Path,将打开一个路径设置窗口。在这个窗口的Path菜单下选Add to Path,找到需要的文件夹,列入MATLAB路径,然后在File菜单下Save Path即可。

MATLAB提供了许多演示程序供使用者参考学习。在MATLAB命令窗下键入demo,将出现MATLAB演示图形窗。使用者可根据提示进行操作。通常画面的上半部是图形,下半部是相应的MATLAB程序语句。使用者可以在界面上直接修改其中的程序语句并执行,观察其结果。因此demo是一个很好的学习辅助手段。

MATLAB语言支持使用DOS命令。在MATLAB命令窗下执行DOS命令,只需在原DOS命令前加!(惊叹号),回车后将直接执行该命令。在用MATLAB语言编写的程序中也可以直接使用!加DOS命令,作为一条MATLAB程序来执行。

 

 

 

 

 

 

 

 

 

§2 MATLAB的基本语法

在MATLAB中,变量和常量的标识符最长允许19个字符,标识符中第一个字符必须是英文字母。MATLAB区分大小写,默认状态下,A和a被认为是两个不同的字符。

一、数组和矩阵

(一)数组的赋值

数组是指一组实数或复数排成的长方阵列。它可以是一维的“行”或“列”,可以是二维的“矩形”,也可以是三维的甚至更高的维数。在MATLAB中的变量和常量都代表数组,赋值语句的一般形式为

    变量=表达式(或数)

如键入a=[1 2 3; 4 5 6; 7 8 9]则将显示结果:

a=

    1    2    3

    4    5    6 

    7    8    9

如键入X=[-3.5 sin(6*pi) 8/5*(3+4) sqrt(2)]则将显示:

X =

    -3.5000   -0.0000   11.2000   1.4142

数组放置在[ ]中;数组元素用空格或逗号“,”分隔;数组行用分号“;”或“回车”隔离。

 

(二)复数

MATLAB中的每一个元素都可以是复数,实数是复数的特例。复数的虚部用i或j表示。

复数的赋值形式有两种:

z=[1+1i ,2+2i ;3+3i ,4+4i]

    z=[1,2;3,4]+[1,2;3,4]*i

得      z=1.000+1.000i    2.000+2.000i

      3.000+3.000i    4.000+4.000i

以上两式结果相同。注意,在第二式中“*”不能省略。

在复数运算中,有几个运算符是常用的。运算符“′”表示把矩阵作共轭转置,即把矩阵的行列互换,同时把各元素的虚部反号。函数conj表示只把各元素的虚部反号,即只取共轭。若想求转置而不要共轭,就把conj和“′”结合起来完成。例如键入

        w=z′,u=conj(z),v=conj(z)′

可得    w=1.000-1.000i    3.000-3.000i

  2.000-2.000i    4.000-4.000i

u=1.000-1.000i    2.000-2.000i

          3.000-3.000i    4.000-4.000i

v=1.000+1.000i    3.000+3.000i

          2.000+2.000i    4.000+4.000i

 

(三)数组寻访和赋值的格式

表M-1 常用子数组的寻访、赋值格式

子数组的寻访

   和赋值

                      使    用    说    明

    a(r,c)

  由a的“r指定行”和“c指定列”上的元素组成的子数组

    a(r,:)

  由a的“r指定行”和“全部列”上的元素组成的子数组

    a(:,c)

  由a的“全部行”和“c指定列”上的元素组成的子数组

   a(:)

  由a的各列按自左到右的次序,首尾相接而生成“一维长列”数组

    a(s)

  “单下标”寻访。生成“s指定的”一维数组。s若是“行数组”(或“列数组”),则a(s)就是长度相同的“行数组”(或“列数组”)

例:a=[1 2 3; 4 5 6; 7 8 9];

键入a(1,2)显示:

ans =

           2

键入a(2,:)显示:

ans =

           4     5     6

键入a(:,3)显示:

ans =

           3

           6

           9

其它情况读者可以自行上机观察使用,此处不再一一举例。

 

(四)执行数组运算的常用函数

表M-2   三角函数和双曲函数

  名  称

    含  义

  名  称

  含  义

  名  称

    含  义

 acos

  反余弦

  asinh

  反双曲正弦

  csch

  双曲余割

  acosh

  反双曲余弦

  atan

  反正切

  sec

  正割

  acot

  反余切

  atan2

 四象限反正切

  sech

  双曲正割

  acoth

  反双曲余切

  atanh

  反双曲正切

  sin

  正弦  

  acsc

  反余割

  cos

  余弦

  sinh

  双曲正弦

  acsch

  反双曲余割

  cosh

  双曲余弦

  tan

  正切

  asec

  反正割

  cot

  余切

  tanh

  双曲正切

  asech

  反双曲正割

  coth

  双曲余切

 

 

  asin

  反正弦

  csc

  余割   

 

 

 

表M-3  指数函数

  名 称

    含  义   

   名  称

含  义

  名  称

  含  义

   exp

  指数  

  log10

  常用对数

  pow2

  2的幂

   log

  自然对数

   log2

 以2为底的对数

  sqrt

  平方根

 

说明:表M-3、表M-4的使用形式与其它语言相似。如

X=tan(60),  Y=20*log(U/0.775),  Z=1-exp(-1.5*t)。

 

表M-4  复数函数

  名  称

含  义

   名  称

含  义

   名  称

  含  义

  abs

 模,或绝对值

  conj

  复数共轭

  real

  复数实部

  angle

 相角(弧度)

  imag

  复数虚部

 

 

例:已知h=a+jb,a=3,b=4,求h的模。

输入:a=3

b=4

h=a+b*j

abs(h)

将显示:

ans =

5

键入:angle(h)

将显示:

ans =

0.    9273

键入:real(h)

将显示:

ans =

3

键入:imag(h)

则显示:

ans =

4

 表M-5  取整函数和求余函数

     名  称

      含  义

      名  称

       含  义

   ceil

  向+∞舍入为整数

   rem(a,b)

  a整除b,求余数  

   fix

  向0 舍入为整数

   round

  四舍五入为整数

   floor

  向-∞舍入为整数

   sign

  符号函数

   mod(x,m)

  x整除m取正余数

 

 

 

例:键入ceil(1.45)

显示:

ans =

          2

键入:fix(1.45)

显示:

ans =

          1

键入:floor(-1.45)

显示:

ans =

     -2

    键入:round(1.45)

显示:

ans =

              1

键入:round(1.62)

显示:

ans =

2

    键入:mod(-55,7)

    显示:

ans =

              1

 键入:rem(-55,7)

 显示:

ans =

-6

 

(五)基本赋值数组

表M-6  常用基本数组和数组运算

                                        基本数组

zeros

全零数组(m×n阶)

logspace   

对数均分向量(1×n阶数组)

ones   

全么数组(m×n阶)

freqspace

频率特性的频率区间

rand

随机数数组(m×n阶)

meshgrid

画三阶曲面时的X,Y网格

randn

 正态随机数数组(m×n阶)

linspace

均分向量(1×n阶数组)

eye(n)

单位数组(方阵)

  将元素按列取出排成一列

                                     特殊变量和函数 

ans

   最近的答案

 Inf

 Infinity(无穷大)  

eps

浮点数相对精度

 NaN

 Not-a-Number(非数)

realmax

   最大浮点实数

 flops

  浮点运算次数 

realmin

最小浮点实数

 computer

 计算机类型 

   pi

3.14159235358579

 inputname *

  输入变量名

   i,j

虚数单位

 size

 多维数组的各维长度

   length 

   一维数组的长度

 

 

 

    为便于大量赋值,MATLAB提供了一些基本数组。举例说明:

A=ones(2,3),B=zeros(2,4),C=eye(3)

得  A=1   1   1         B=0  0  0  0         C=1  0  0

      1   1   1            0  0  0  0           0  1  0

                             0  0  1 

线性分割函数linespace(a,b,n)在a和b之间均匀地产生n个点值,形成1×n元向量。如

D=linspace(0,1,5)

得  D= 0  0.2500  0.5000  0.7500  1.0000

 

(六)数组运算和矩阵运算

    MATLAB中最基本的运算是矩阵运算。但是在MATLAB的运用中,大量使用的是数组运算。从外观形状和数据结构上看,二维数组和(数学中的)矩阵没有区别。但是,矩阵作为一种变换或映射算子的体现,其运算有着明确而严格的数学规则。而数组运算是MATLAB软件所定义的规则,其目的是为了数据管理方便、操作简单、指令形式自然简便以及执行计算有效。虽然数组运算尚缺乏严谨的数学推理,数组运算本身仍在完善和成熟中,但它的作用和影响正随着MATLAB的发展而扩大。

为更清晰地表述数组运算与矩阵运算的区别,我们以表M-7叙述各数组运算指令的意义。其中假定S=2,n=3,P=1.5。

A=[1 2 3; 4 5 6; 7 8 9],

B=[2 3 4; 5 6 7; 8 9 1]。

表M-7  举例说明数组运算指令的意义

  指 令

       含 义 

         运 算 结 果

  s+A

 标量s分别与A元素之和

     3     4     5

     6     7     8

     9    10    11

  A-s

 A分别与标量s的元素之差

-1     0     1

     2     3     4

     5     6     7

  s.*A

 标量s分别与A的元素之积

     2     4     6

     8    10    12

    14    16    18

s./A或

A.\s

 s分别被A的元素除

    2.0000    1.0000    0.6667

    0.5000    0.4000    0.3333

    0.2857    0.2500    0.2222

  A.^n

 A的每个元素自乘n次

     1     8    27

    64   125   216

   343   512   729

  p.^A

 以p为底,分别以A的元素为指数

求幂值

1.5000    2.2500    3.3750

    5.0625    7.5938   11.3906

   17.0859   25.6289   38.4434

  A+B

 对应元素相加

     3     5     7

     9    11    13

    15    17    10

  A-B

 对应元素相减

-1    -1    -1

    -1    -1    -1

    -1    -1     8

  A.*B

 对应元素相乘

     2     6    12

    20    30    42

    56    72     9

  A./B或 B.\A

  A的元素被B的对应元素除

    0.5000    0.6667    0.7500

    0.8000    0.8333    0.8571

    0.8750    0.8889    9.0000

  exp(A)

 以自然数e为底,分别以A的元素为指数,求幂

 1.0e+003 *

    0.0027    0.0074    0.0201

    0.0546    0.1484    0.4034

    1.0966    2.9810    8.1031

  log(A)

 对A的各元素求对数

          0    0.6931    1.0986

    1.3863    1.6094    1.7918

    1.9459    2.0794    2.1972

  sqrt(A)   

 对A的各元素求平方根

    1.0000    1.4142    1.7321

    2.0000    2.2361    2.4495

    2.6458    2.8284    3.0000

例:有一函数X(t)=tsin3t,在MATLAB程序中如何表示?

解:  X=t.*sin(3*t)

例:有一函数X(t)=sin3t/3t, 在MATLAB程序中如何表示?

解:  X=sin(3*t)./(3*t)

二、逻辑判断与流程控制

(一)关系运算

关系运算是指两个元素之间数值的比较,一共有六种可能。如表M-8所列。

关系运算的结果只有两种可能,即0或1。0表示该关系式为“假”,1表示该关系式为“真”。

例1:A=3+4==7,得  A=1。

例2:已知N=0,B=[N==0],得  B=1。

      若N=2,B=[N==0],得  B=0。

表M-8  关系运算符

    指  令

     含  义

     指  令

     含  义

      <

     小于

      >=

     大于等于

      <=

     小于等于

      ==

     等于

      >

     大于

      ~=

     不等于

(二)逻辑运算

逻辑量的基本运算为“与(&)”、“或(∣)”、“非(~)”三种,另外还可以用“异或(xor)”,如表M-9所示。

表M-9  逻辑运算符

 

  运  算

          A=0

             A=1

   B=0

    B=1

    B=0

    B=1

  A&B

      0

     0

     0

      1

  A|B

      0

     1

     1

      1

  ~A

      1

     1

     0

      0

  xor(A,B)

      0

     1

     1

      0

 

  (三)基本的流程控制语句

⑴ if条件执行语句

格式:  if  表达式 语句, end

       if  表达式1  语句组A, else   语句组B, end

if  表达式1  语句组A, elseif   表达式2  语句组B, else  语句组C, end

执行到该语句时,计算机先检验if后的逻辑表达式,为1则执行语句A;如为0则跳过A检验下一句程序,直到遇见end,执行end后面的一条语句。

例:if n<=2

       x=2;

    elseif n>3

       x=3;

    end

若n=5,则结果

   x=

3    

⑵ while循环语句

格式:  while 表达式  语句组A, end

执行到该语句时,计算机先检验while后的逻辑表达式,为1则执行语句A;到end处,它就跳回到while的入口,再检验表达式,如仍为1则再执行语句A,直到结果为0,就跳过语句组A,直接执行end后面的一条语句。

例:while k<=1000

         k=k+1;

    end

键入k将显示

    k=

      1001

⑶ for循环语句

格式:  for k=初值:增量:终值  语句组A, end

将语句组A重复执行N次,但每次执行时程序中k值不同。增量缺省值为1。

例:y=0;

for k=1:20

         y=y+k;

end

键入y将显示

y=

210

⑷ switch多分支语句

格式:  switch 表达式(标量或字符串)

         case  值1

           语句组A

         case  值2

           语句组B

         ……………

         otherwise

            语句组N

         end

当表达式的值与某case语句中的值相同时,它就执行该case语句后的语句组,然后直接跳到终点的end处。

三、基本绘图方法

(一) 二维图形函数

MATLAB语言支持二维和三维图形,这里我们主要介绍常用的二维图形函数。如表M-11所示。

表M-11  常用图形函数库

                                 基本X—Y图形  

plot

 线性X-Y座标绘图

 polar

 极座标绘图

 loglog

 双对数X-Y座标绘图

 plotyy

 用左、右两种Y座标画图

 semilogx

 半对数X座标绘图

 semilogy.

 半对数Y座标绘图

 stem

 绘制脉冲图

 stairs

 绘制阶梯图

 bar

 绘制条形图

 

 

                                  坐 标 控 制 

 axis

控制座标轴比例和外观

subplot  

 按平铺位置建立子图轴系

 hold

 保持当前图形

 

 

                                  图 形 注 释

 title

 标出图名(适用于三维图形)

 gtext

 用鼠标定位文字

 xlabel

 X轴标注(适用于三维图形)

 legend

 标注图例

 ylabel

 Y轴标注(适用于三维图形)

 grid

 图上加座标网格(适用于三维)

 text

 在图上标文字(适用于三维) 

 

 

                                     打 印 

 print

 打印图形或把图存为文件

 orient

 设定打印纸方向

 printopt

打印机默认选项

 

  

                             常用的三维曲线绘图命令

 Plot3

在三维空间画点和线

 mesh

 三维网格图

 fill3

在三维空间绘制填充多边形

 surf

 三维曲面图

最常用的命令使用说明:

plot(t,y)表示用线性X-Y座标绘图,X轴的变量为t,Y轴的变量为y。

subplot(2,2,1)建立2×2子图轴系,并选定图1。

axis([0  1  -0.1  1.2])表示建立一个座标,横座标的范围从0至1,纵座标的范围从-0.1至1.2。

title('X(n)曲线')在子图上端标注图名

作图时,线形、点形和颜色的选择可参考表M-12。

表M-12 线形、点形和颜色

标志符

  b

  c

  g

  k

  m

  r

  w

  y

 

颜  色

  蓝

  青

  绿

  黑

 品红

  红

  白

  黄

 

标志符

  ?

  。

  ×

  +

  -

  ﹡

  ﹕

  -.

 ┄

线、点

  点

 圆圈

 ×号

 +号

 实线

 星号

 点线

点划线

 虚线

 

(二)举例

以下举例说明二维图形函数在程序中的使用方法。

例1:作一条曲线 ,程序如下。

t=0:0.5:4*pi;                      %将t在0到4π间每间隔0.5取一点

y=exp(-0.1*t).*sin(t);

subplot(2,2,1),plot(t,y);       %建立2×2子图轴系,在图1处绘线性图

title('plot(t,y)');              %标注图名

subplot(2,2,2),stem(t,y);       %在2×2子图轴系图2处绘脉冲图

title('stem(t,y)');

subplot(2,2,3),stairs(t,y);     %在2×2子图轴系图3处绘阶梯图

title('stairs(t,y)');

subplot(2,2,4),bar(t,y);        %在2×2子图轴系图2处绘条形图

title('bar(t,y)');

例2:已知 , 。在同一座标系对两条曲线作图,用不同的颜色和线型区分。

方法一:将同时显示曲线的向量列入数组,t必须等长。显示的线型和颜色不能任意选择。

t=0:0.01:2;                  

y1=sin(2*pi*t);

y2=cos(4*pi*t);

plot(t,[y1;y2]);

图M-1 例2方法一

方法二:显示曲线的向量t不必等长,显示的线型和颜色能任意选择。作图时,先画第一条曲线保持住,再画第二条曲线。

 

  图M-2 例2方法二

t1=0:0.01:1;                  

y1=sin(2*pi*t1);

t2=0:0.01:2

y2=cos(4*pi*t2);

plot(t1,y1,'*m'),hold ;        %让第一条曲线保持住,再画第二条曲线

plot(t2,y2,'+b');      

 

 

 

 

 

 

 

 

 

 

 

 

 

 

§3 MATLAB在信号处理中常用的函数

MATLAB系统软件中具有专用的信号处理工具箱,对于我们学习信号与系统、数字信号处理等课程,进行通信、电子工程设计计算是一个非常有效的辅助手段。这里,我们仅列写出最常用的部分,供大家参考。

一、常用的信号及信号的波形

(一)常用的信号

    在 MATLAB中的信号处理工具箱中,主要提供的信号是离散信号。

由于MATLAB对下标的约定为从1开始递增,例如x=[5,4,3,2,1,0],表示x(1)=5,x(2)=4,X(3)=3…

因此要表示一个下标不由1开始的数组x(n),一般应采用两个矢量,如

        n=[-3,-2,-1,0,l,2,3,4,5];  

        x=[1,-l,3,2,0,4,5,2,1];

这表示了一个含9个采样点的矢量:X(n)={x(-3),x(-2),x(-1),x(0),x(1),x(2),x(3),x(4),x(5)}。

    1.单位取样序列   

                       

这一函数实现的方法有二:

方法一:可利用MATLAB的zeros函数。

        x=zeros(1,N);    %建立一个一行N列的全零数组

   x(1)=1;            %对X(1)赋1

    方法二:可借助于关系操作符实现

   n=1:N;             

   x=[n==1];  %n等于1时逻辑关系式结果为真,x=1;n不等于1时为假,x=0

如要产生              

则可采用MATLAB实现:

       n=n1:n2;

  x=[(n-n0)==0];%n=n0时逻辑关系式结果为真,x=1;n≠n0时为假,x=0

    2.单位阶跃序列

                     

    这一函数可利用MATLAB的ones函数实现:

        x=ones(1,N);

还可借助于关系操作符“>=”来实现。如要产生在n1≤n0≤n2上的单位阶跃序列

                      

则可采用MATLAB实现:

        n=n1:n2;

   x=[(n-n0)>=0]; %n-n0≥0为真,x=1;n-n0<0时为假,x=0

    3.实指数序列

采用MATLAB实现:

        n=0:N-l;

        x=a.^n; 

    4.复指数序列

采用MATLAB实现:

        n=0:N-1;

        x=exp((lu+j*w0)*n);

    5.正(余)弦序列

                                     

采用MATLAB实现:

        n=0:N-l;

        x=cos(w0*n+Q);

    6.随机序列

    MATLAB中提供了两类(伪)随机信号:

    rand(1,N)产生[0,1)上均匀分布的随机矢量;

    randn(1,N)产生均值为0,方差为1的高斯随机序列,也就是白噪声序列。其它分布的随机数可通过上述随机数的变换而产生。

    7.周期序列

                     

    例如,设t1表示T序列中一个周期的序列,要产生4个周期的T序列,用MATLAB实现:

T=[t1 t1 t1 t1];

 

(二)信号波形

    l. sawtooth

    功能:产生锯齿波或三角波。

    格式:x=sawtooth(t)

          x=sawtooth(t,width)

    说明:sawtooth(t)函数类似于sin(t),它产生周期为2π,幅值从-1到+1的锯齿波,在2π的整数倍处,其值为-1,并以1/π的斜率线性上升至+1。

    sawtooth(t,width)用于产生三角波,其width(O<width≤1的标量)用于确定最大值的位置,即从0到2π*width函数从-1上升到+1,然后在2π*width至2π之间又线性地从+1降至-1,周而复始。例如,当width=0.5时,可产生一对称的标准三角波,当width=1时,就产生锯齿波,即sawtooth(t,1)=sawtooth(t)。

   

2.square

    功能: 产生方波。

    格式: x=square(t)

          x=square(t,duty)

说明:square(t)类似于sin(t),产生周期是2π,幅值为士l的方波,square(t,duty)产生指定周期的方波,其duty用于指定正半周期的比例。

 

⒊ sinc

功能:产生sinc或 函数。

格式:y=sinc(x)

说明: MATLAB中的sinc子函数可用于计算下列函数

这个函数是宽度为2π,幅度为1的矩形脉冲的连续逆傅里叶变换,即

 

⒋ diric

功能: 产生Dirichlet或周期sinc函数。

格式: y=diric(x,n)

说明:在 y=diric(x,n)中,n必须为正整数,y为相应的x元素的Dirichlet函数,即

Dirichlet函数是周期信号,当n为奇数时,周期为2π;当n为偶数时,周期为4π。

 

二、多项式运算常用函数

多项式在近代信息和控制理论中具有十分重要的位置。MATLAB中提供了解决这些问题的工具,如表M-13所示。

多项式的直接表达方式:

    MATLAB约定,降幂多项式 其系数行向量表达式为:

表M-13  多项式运算子函数

                               多     项     式

roots

多项式求根

polyfit

用多项式曲线拟合数据

poly

按根组成多项式

polyder

多项式求导数

polyval

多项式求值

poly2str

由系数行向量求符号多项式

polyvalm

矩阵作变元的多项式求值

conv

多项式相乘,卷积

residue

部分分式展开(留数)

deconv

多项式相除,反卷积

 

在程序中可以按上式的格式直接输入。注意:多项式系数应以降幂次序排列。假如缺项,则认为该项系数为零。

若要表示 ,可建立

,再利用指令: 。多项式P是一个特征多项式, 的元素被认为是多项式P的根。

例1:已知多项式的根向量R,求其特征多项式P和以降幂形式表示的S的多项式。

R=[2 –2 3];               %建立根向量

P=poly(R)                  %求R的特征多项式

PR=poly2str(P,’S’)       %求以降幂形式表示的S多项式         

程序运行结果:

P =  1    -3    -4    12

PR =

  S^3 - 3 S^2 - 4 S + 12

注意:⒈ 特征多项式P的系数向量一定是(n+1)维的。

⒉ 特征多项式系数向量的第一个元素必是1。

例2:由给定的根向量R求X多项式。

R=[-4,–1+4*i,-1-4*i];          %建立根向量

P=poly(R)                      %求R的特征多项式

PR=real(P)                     %求P的实部

PPR=poly2str(PR,’X’)            %求以降幂形式表示的X多项式

程序运行结果:

P =

     1     6    25    68

PR =

     1     6    25    68

PPR =

   X^3 + 6 X^2 + 25 X + 68

注意:⒈ 实系数多项式要求根向量中的复数根必须共轭成对。

⒉ 含复数的根向量所生成的多项式系数向量(如P),其系数可能带有截断误差数量级的虚部。因而采用“real”指令取实部,将很小的虚部滤掉。

三、快速傅里叶变换函数

1.fft

    功能:一维快速傅里叶变换(FFT)。

    格式: y=fft(x)

           y=fft(x,n)

    说明:fft函数用于计算矢量或矩阵的离散傅里叶变换,这可通过

实现,其中 , 。

    y=fft(x)为利用FFT算法计算矢量x的离散傅里叶变换,当x为矩阵时,y为矩阵x每一列的FFT。当x的长度为2的幂次方时,则fft函数采用基2的FFT算法,否则采用稍慢的混合基算法。

y=fft(x,n)采用n点FFT。当x的长度小于n时,fft函数在x的尾部补零,以构成n点数据;当x的长度大于n时,fft函数会截断序列x。当x为矩阵时,fft函数按类似的方式处理列长度。

 

    2. ifft

功能:一维逆快速傅里叶变换(IFFT)。

格式: y=ifft(x)

          y=ifft(x,n)

    说明: ifft函数用于计算矢量或矩阵的逆傅里叶变换,即

                

其中 , 。

     y=ifft(x)用于计算矢量x的IFFT。当x为矩阵时,计算所得的y为矩阵x中每一列的IFFT。

     y=ifft(x,n)采用n点IFFT。当length(x)<n时,在x中补零;当length(x)>n时,将x截断,使length(x)=n。

ifft函数是fft函数的逆,其相应的M文件中采用的算法类似。

 

3.cftbyfft

功能:cftbyfft子函数采用FFT计算连续时间Fourier变换。

格式:[AW,f]=cftbyfft(wt,t,flag)

输出幅频谱数据对为[AW,f],输入量(W,t)为已经窗口化了的时间函数W(t),它们分别是长度为N的向量。flag取非0时(缺省使用),频率范围在[0,Fs];flag取0时,频率范围在[-Fs/2,Fs/2]。

四、系统分析与实现函数

1.  abs

功能:求绝对值(幅值)。

格式:y=abs(x)

说明:y=abs(x)用于计算x的绝对值。当x为复数时,得到的是复数模(幅值),即:

    

当x为字符串时,abs(x)得到字符串的各个字符的ASCⅡ码,例如x=‘123’,则abs(x) 得到49 50 51。

 

2.angle

功能:求相角。

格式:p=angle(h)

说明:p=angle(h)用于求取复矢量或复矩阵H的相角(以弧度为单位),相角介于-π和+π之间。例如,某复数h可用两种方法表示:

            

则幅值m和相角φ可由x+jy格式求出          

             m=abs(h)

             φ=angle(h)

当然,由m和φ也可求取x+jy格式

             h=m.*exp(j*φ)

             x=real(h)

             y=imag(h)

 

3.conv

功能:求卷积

格式:c=conv(a,b)

说明:conv(a,b)用于求取矢量a和b的卷积,即

              

其中N为矢量a和b的最大长度。例如,当a=[1 2 3],b=[ 4 5 6]时,则

               c=conv(a,b)

               c=

                  4 13 28 27 18

此函数可直接用于求两个有限长序列的卷积。设x(n)和h(n)的长度分别为M和N,则

             y=conv(x,h)

y的长度分别为N+M-1。

 

4.filter

功能:利用IIR滤波器和FIR滤波器对数据进行滤波。

格式:y=filter(b,a,x)

     [y,zf]=filter(b,a,x)

     y=filter(b,a,x,zi)

说明:filter采用数字滤波器对数据进行滤波,其实现采用移位直接Ⅱ型结构,因而适用于IIR和FIR滤波器。滤波器的系统函数为  

        

         即滤波器系数a=[a0 a1 a2 ...an],b=[b0 b1 ...bm],输入序列矢量为x。这里,标准形式为a0=1,如果输入矢量a时,a0≠1,则MATLAB将自动进行归一化系数的操作;如果a0=0,则给出出错信息。

y=filter(b,a,x)利用给定系数矢量a和b对x中的数据进行滤波,结果放入y矢量中,y的长度取max(N,M)。

y=filter(b,a,x,zi)可在zi中指定x的初始状态。

[y,zf]=filter(b,a,x)除得到矢量y外,还得到x的最终状态矢量zf。

 

5.freqs

功能:连续时间系统的频率响应。

格式:h=freqs(b,a,w)

[h,w]=freqs(b,a)

[h,w]=freqs(b,a,n)

      freqs(b,a)

说明:freqs用于计算由矢量a和b构成的连续时间系统

         

的复频响应H(jω)。其系统函数的系数a=[a0 a1 a2 ...an],b=[b0 b1 ...bm]。

h=freqs(b,a,w)用于计算连续时间系统的复频响应,其中实矢量w用于指定频率值。

[h,w]=freqs(b,a)自动设定200个频率点来计算频率响应,其200个频率值记录在w中。

[h,w]=freqs(b,a,n) 设定n个频率点计算频率响应。

不带输出变量的freqs函数,将在当前图形窗口中描绘幅频和相频曲线。

 

⒍ freqz  

功能:离散时间系统的频率响应。

格式:[h,w]=freqz(b,a,n)

      [h,f]=freqz(b,a,n,Fs)

      h=freqz(b,a,w)

      h=freqz(b,a,f,Fs)

      freqz(b,a,n)

说明:  freqz 用于计算数字滤波器H(Z)的频率响应函数H(ejω)。

[h,w]=freqz(b,a,n)可得到数字滤波器的n点复频响应值,这n个点均匀地分布在[0,π]上,并将这n个频点的频率记录在w中,相应的频响值记录在h中。要求n为大于零的整数,最好为2的整数次幂,以便采用FFT计算,提高速度。缺省时n =512。

    [h,f]=freqz(b,a,n,Fs)用于对H(ejω)在[0,Fs/2]上等间隔采样n点,采样点频率及相应频响值分别记录在f 和h中。由用户指定FS(以HZ为单位)值。

h=freqz(b,a,w)用于对H(ejω)在[0,2π]上进行采样,采样频率点由矢量w指定。

h=freqz(b,a,f,Fs) 用于对H(ejω)在[0,FS]上采样,采样频率点由矢量f指定。

freqz(b,a,n) 用于在当前图形窗口中绘制幅频和相频特性曲线。

 

7.freqz_m

功能:离散时间系统的幅值响应、相位响应及群迟延响应。

格式:[db,mag,pha,grd,w]= freqz_m(b,a)

说明:freqz_m函数是freqz函数的修正函数,可获得幅值响应(绝对和相对)、相位响应及群迟延响应。式中

db中记录了一组对应[0,π] 频率区域的相对幅值响应值;

mag中则记录了一组对应[0,π] 频率区域的绝对幅值响应值;

pha中记录了一组对应[0,π] 频率区域的相位响应值;

grd中记录了一组对应[0,π] 频率区域的群迟延响应值;

w中记录了对应[0,π] 频率区域的501个频点的频率值。

 

⒏ impulse

功能:求解连续系统的冲激响应。

    格式:impulse(b,a)

          impulse(b,a,t)

          y=impulse(b,a,t)

    说明:impulse用于计算由矢量a和b构成的连续时间系统

         

的冲激响应。其系统函数的系数a=[a0 a1 a2 ...an],b=[b0 b1 ...bm]。

impulse(b,a)计算并显示出连续系统的冲激响应h(t)的波形,其中t将自动选取。

impulse(b,a,t)可由用户指定t值。若t为一实数,将显示连续系统在0~t秒间的冲激响应波形;若t为数组例如[t1:dt:t2],则显示连续系统在指定时间t1~t2内的冲激响应波形,时间间隔为dt。

y=impulse(b,a,t)将结果存入输出变量y,不直接显示系统冲激响应波形。

 

⒐ impz

    功能:求解离散系统的冲激响应。

    格式:[h,t]=impz(b,a)

          [h,t]=impz(b,a, n)

          [h,t]=impz(b,a,n,Fs)

          impz(b,a)

    说明:由矢量a和b构成离散时间系统,即

         

    [h,t]=impz(b,a)计算出离散系统的冲激响应h,取样点数n由impz函数自动选取,并记录在矢量t中(t=[0:n-1]’)。

    [h,t]=impz(b,a, n)可由用户指定取样点或取样时刻。当n为标量时,t=[0:n-1]’,即在0~ n-1时刻计算冲激响应,0时刻表示离散系统的起始点;当n为矢量(其值应为整数),则表示t= n,即在这些指定的时刻计算冲激响应。

    [h,t]=impz(b,a,n,Fs) 表示取样间隔为1/Fs,在缺省Fs时,则取为1。

不带输出变量的impz将在当前图形窗口中利用stem(t,h)函数绘出冲激响应。

 

⒑ lsim

功能:对连续系统的响应进行仿真。

格式:lsim(b,a,x,t)

      y= lsim(b,a,x,t)

说明:a=[a0 a1 a2 ...an],b=[b0 b1 ...bm]是连续时间系统的传递函数的系数。

x和t是系统输入信号的行向量。例如

t=0:0.01:10;

x=sin(t);

定义输入信号为一正弦信号sin(t),且这个信号在0~10秒的时间内每间隔0.01秒选取一个取样点。

lsim(b,a,x,t)当将输入信号加在由a、b所定义的连续时间系统输入端时,将显示系统的零状态响应的时域仿真波形。

y= lsim(b,a,x,t)当将输入信号加在由a、b所定义的连续时间系统输入端时,不直接显示仿真波形,而是将求出的数值解存入输出变量y。

 

⒒ step

功能:求解连续系统的阶跃响应。

格式:step(b,a)

      step(b,a,t)

      y=step(b,a,t)

说明:step用于计算由矢量a和b构成的连续时间系统

         

的阶跃响应。其系统函数的系数a=[a0 a1 a2 ...an],b=[b0 b1 ...bm]。

step(b,a)计算并显示出连续系统的阶跃响应g(t)的波形,其中t将自动选取。

step(b,a,t)可由用户指定t值。若t为一实数,将显示连续系统在0~t秒间的阶跃响应波形;若t为数组例如[t1:dt:t2],则显示连续系统在指定时间t1~t2内的阶跃响应波形,时间间隔为dt。

y=step(b,a,t)将结果存入输出变量y,不直接显示系统阶跃响应波形。

五、IIR滤波器设计函数

⒈ butter

功能:Butterworth(巴特沃斯)模拟和数字滤波器设计。

格式:[b,a]=butter(n,Wn)

      [b,a]=butter(n,Wn,’ftype’)

      [b,a]=butter(n,Wn,’s’)

      [b,a]=butter(n,Wn, ’ftype’,’s’)

      [z,p,k]=butter(...)

      [A,B,C,D]=butter(...)

说明:butter函数可设计低通、带通、高通和带阻的数字和模拟滤波器,其特性为使通带内的幅度响应最大限度地平坦,这会损失截止频率处的下降斜率。因此,在期望通带平滑的情况下,可使用butter函数,但在期望下降斜率大的场合,应使用椭圆和Chebyshev滤波器。

⑴数字域

[b,a]=butter(n,Wn)可设计出截止频率为Wn的n阶低通Butterworth滤波器,其滤波器为 

   

截止频率是滤波器幅度响应下降至 处的频率,Wn∈[0,1],其中1相应于0.5Fs(取样频率,即Nyquist频率)。

    当Wn=[W1 W2] (W1<W2)时,butter函数产生一2 n阶的数字带通滤波器,其通带为W1<ω<W2。

[b,a]=butter(n,Wn,’ftype’)可设计出高通或带阻滤波器。

·当ftype=high时,可设计出截止频率为Wn的高通滤波器;

    ·当ftype=stop时,可设计出带阻滤波器,这时Wn=[W1 W2],且阻带为W1<ω<W2。

    利用输出变量个数的不同,可得到滤波器的另外两种表示:零极点增益和状态方程。

·[z,p,k]=butter(n,Wn)或[z,p,k]=butter(n,Wn,’ftype’) 可得到滤波器的零极点增益表示;

·[A,B,C,D]=butter(n,Wn)或[A,B,C,D]=butter(n,Wn,’ftype’) 可得到滤波器的状态空间表示。

⑵模拟域

    [b,a]=butter(n,Wn,’s’) 可设计出截止频率为Wn的n阶低通模拟Butterworth滤波器,

             

其中截止频率Wn>0。

模拟域的butter函数说明与数字域完全相同,也有六种形式,此处不再赘述。

   

    ⒉ cheby1

功能:Chebyshev(切比雪夫)Ⅰ型滤波器设计(通带等波纹)。

    格式:[b,a]= cheby1(n,Rp,Wn)

      [b,a]= cheby1(n,Rp,Wn, ’ftype’)

      [b,a]= cheby1(n,Rp,Wn, ’s’)

      [b,a]= cheby1(n,Rp,Wn, ’ftype’,’s’)

      [z,p,k]= cheby1(...)

      [A,B,C,D]= cheby1(...)

说明:cheby1函数可设计低通、带通、高通和带阻的数字和模拟ChebyshevⅠ型滤波器,其通带内为等波纹,阻带内为单调。ChebyshevⅠ型滤波器的下降斜率比Ⅱ型大,但其代价是在通带内波纹较大。

与butter函数类似,cheby1函数可设计数字域和模拟域的ChebyshevⅠ型滤波器。其通带内的波纹由Rp(分贝)确定。其它各公式的使用方法与butter函数相同。

 

⒊ cheby2

    功能:Chebyshev(切比雪夫)Ⅱ型滤波器设计(阻带等波纹)。

    格式:[b,a]= cheby2(n,As,Wn)

      [b,a]= cheby2(n,As,Wn, ’ftype’)

      [b,a]= cheby2(n,As,Wn, ’s’)

      [b,a]= cheby2(n,As,Wn, ’ftype’,’s’)

      [z,p,k]= cheby2(...)

      [A,B,C,D]= cheby2(...)

说明:cheby2函数可设计低通、带通、高通和带阻的数字和模拟ChebyshevⅡ型滤波器,与cheby1函数几乎一样,只不过cheby2函数其通带内为单调,阻带内为等波纹。因此,由As确定阻带内的波纹。

其它各公式的使用方法与butter函数相同,可参考相应公式。

六、窗函数

1. boxcar

功能:矩形窗

格式:w=boxcar(n)

说明:boxcar(n)函数可产生一长度为n的矩形窗函数。

 

2.triang

功能:三角窗

格式:w=triang(n)

说明:triang(n)函数可得到n点的三角窗函数。三角窗系数为

当n为奇数时: 

          

当n为偶数时:      

三角窗函数非常类似于Bartlett窗。Bartlett窗在取样点1和n上总是以零结束,而三角窗在这些点上并不为零。实际上,当n为奇数时,triang(n-2)的中心n-2个点等效于Bartlett(n)。

 

3.Bartlett

功能:Bartlett(巴特利特)窗。

格式:w=Bartlett(n)

说明:Bartlett(n)可得到n点的Bartlett窗函数。Bartlett窗函数系数为                     

Bartlett函数与三角窗窗非常类似,当n为奇数时,Bartlett(n)的中心n-2个点等效于triang(n-2)。

 

⒋ hamming

功能:Hamming(哈明)窗。

格式:w=hamming(n)

说明:hamming(n)可产生n点的Hamming窗,其系数为

    

⒌ Hanning

功能:hanning(汉宁)窗。

格式:w=hanning(n)

说明:hanning(n)可产生n点的Hanning窗,其系数为

          

 

⒍ blackman

功能:Blackman(布莱克曼)窗。

格式:w=blackman(n)

说明:blackman(n)可产生n点的Blackman窗,其系数为

    

与等长度的Hamming和Hanning窗相比,Blackman窗的主瓣稍宽,旁瓣稍低。

 

⒎ chebwin

功能:Chebyshev(切比雪夫)窗。

格式:w=chebwin(n,r)

说明:w=chebwin(n,r)可产生n点的Chebyshev窗函数,其傅里叶变换后的旁瓣波纹低于主瓣rdB。注意:当n为偶数时,窗函数的长度为n+1。

 

⒏ kaiser

功能:Kaiser(凯塞)窗。

格式:w=kaiser(n,beta)

说明:w=kaiser(n,beta)可产生n点的Kaiser窗函数,其中beta为影响窗函数旁瓣的β参数,其最小的旁瓣抑制α与β之间的关系为

增加β可使主瓣变宽,旁瓣的幅度降低。

七、 FIR数字滤波器设计函数

1. fir1

功能:基于窗函数的FIR数字滤波器设计——标准频率响应。

格式:b=fir1(n,Wn)

      b=fir1(n,Wn,'ftype')

      b=fir1(n,Wn,Window)

      b=fir1(n,Wn,'ftype',Window)

说明:fir1函数以经典方法实现加窗线性相位FIR滤波器设计,它可设计出标准的低通、带通、高通和带阻滤波器。

b=fir1(n,Wn)可得到n阶低通FIR滤波器,滤波器系数包含在b中,这可表示成:

 

这是一个截止频率为Wn的Hamming(汉明)加窗线性相位滤波器,0≤Wn≤1,Wn=1相应于0.5fs。

当Wn=[W1 W2]时,fir1函数可得到带通滤波器,其通带为W1<ω< W2。

b=fir1(n,Wn,'ftype')可设计高通和带阻滤波器,由ftype决定:

·当ftype=high时,设计高通FIR滤波器;

·当ftype=stop时,设计带阻FIR滤波器。

在设计高通和带阻滤波器时,fir1函数总是使用阶为偶数的结构,因此当输入的阶次为奇数时,fir1函数会自动加1。这是因为对奇数阶的滤波器,其在Nyquist频率处的频率响应为零,因此不适合于构成高通和带阻滤波器。

b=fir1(n,Wn,Window)则利用列矢量Window中指定的窗函数进行滤波器设计,Window长度为n+1。如果不指定Window参数,则fir1函数采用Hamming窗。

b=fir1(n,Wn,'ftype',Window)可利用ftype和Window参数,设计各种加窗的滤波器。

由fir1函数设计的FIR滤波器的群延迟为n/2。

 

2. fir2

功能:基于窗函数的FIR滤波器设计——任意频率响应。

格式:b=fir2(n,f,m)

      b=fir2(n,f,m,Windows)

      b=fir2(n,f,m,npt)

      b=fir2(n,f,m,npt,Windows)

说明:fir2函数用于设计任意频率响应的加窗数字FIR滤波器,对标准的低通、带通、高通和带阻滤波器的设计可采用fir1函数。

b=fir2(n,f,m)可设计出一n阶的FIR滤波器,其滤波器的频率特性由矢量f和m决定。

·f为频率点矢量,且f∈[0,1],当f=1时就相应于0.5fs。矢量f中按升序排列,且第一个必须为0,最后一个必须为1,并允许出现相同的频率值;

·矢量m中包含与f相对应的期望滤波器响应幅度;

·矢量f和m的长度必须相同。

b=fir2(n,f,m,Windows)可将列矢量Windows中指定的窗函数用于滤波器设计,如省略Windows,则自动选取Hamming窗。

b=fir2(n,f,m,npt)格式中,可利用参数npt指定fir2对频率响应进行内插的点数,对应的b=fir2(n,f,m,npt,Windows)格式中,可指定窗函数。

 

  评论这张
 
阅读(1667)| 评论(2)

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018