第6章 信号处理

null第6章 信号处理 第6章 信号处理 6.1 MATLAB的数据可视化达 6.2 MATLAB仿真中的信号处理 6.3 信号特征参数的计算 6.4 信号的频谱分析 6.5 线性系统 6.6 数字滤波器设计与实现6.1 MATLAB的数据可视化表达 6.1 MATLAB的数据可视化表达 6.1.1 二维图形的绘制 1.常用的二维图形绘图函数 基本的二维绘图函数有： ·plot：绘制二维曲线； ·title：给图形加标题； ·grid：显示网格线； ·xlabel：给x轴加标记； ·ylabel：给y轴加标记； ·text：在坐标图中加入文字注释。 null ［例1］画出函数x=sin2πx的曲线图，其中x从0到2π,步进为π/100。程序如下： 程序6-1   %ch6plot.m X=0:pi/100:2*pi; Y=sin(X); plot(X,Y);%作图 gridon;%网格线显示，若改为gridoff,则不显示网格 ylabel(′y=sin2＼pix′);%Y轴标注，可以有汉字 xlabel(′x′);%X轴标注，可以有汉字null title(′functionploty=sin2＼pix′);%图标题 text(0.5,sin(0.5),′＼leftarrowsin2＼pi0.5′); %指令text可以在指定坐标处写文字标注 text(2.3,sin(2.3),′＼leftarrowsin2＼pi2.3′);%所有标注中均可使用汉字 %对于特殊符号，如希腊字母、箭头等需要采用LaTeX格式 结果如图6-1所示。  null图6-1 基本的二维绘图函数用法 null 2.图形的线型和颜色控制 在命令“plot”的高级用法中，可以设置作图的线型，标记类型，线和标记的颜色、粗细等特征。用命令“docLineSpec”和“docplot”可以查询详细的帮助文档。作图时常用的线型、标记以及颜色的定义参见表6-1。 null表6-1 作图时常用的线型、标记以及颜色的定义null 利用线型属性和标记属性可以随心所欲地设计作图图式。常用的线型属性和标记属性参见表6-2。 null表6-2 常用的线型属性和标记属性null 程序6-2 t=0:pi/20:2*pi; plot(t,sin(t),′-.r*′)%用红色点画线和星号作图 holdon%保持当前图形不被擦除 plot(sin(t-pi/2),′--bp′)%用蓝色虚线和五角星标记作图 plot(sin(t-pi),′:ks′)%用黑色虚点线和方框标记作图 holdoff 结果如图6-2所示(黑白印刷时颜色表现不出来)。 null图6-2 二维绘图线型和标记的特征控制(1)null ［例3］ 程序6-3   t=0:pi/20:2*pi; plot(t,sin(2*t),′-mo′,...%线型:实线,洋红色,小圆标记 ′LineWidth′,2,...%线宽为2 ′MarkerEdgeColor′,′k′,...%标记边缘颜色:黑色 ′MarkerFaceColor′,［.491.63］,...%标记面颜色:淡绿 ′MarkerSize′,12)%标记大小:12  结果如图6-3所示。null图6-3 二维绘图线型和标记的特征控制(2)null 3.图形的标注 图形的标注可以用“text”函数。其用法是： text(x,y,′字符串′); text(...′PropertyName′,PropertyValue...);   其中，字符串中若有特殊符号，如希腊字母、箭头等，需要采用LaTeX格式表示。′PropertyName′为字符属性名称，PropertyValue为相应的属性取值。用“doctext-props”null 可以得到详细的用法手册。对于常用的一些，举例如下： ·改变字符大小，属性为′FontSize′，取值为：10，12，16等。 ·改变字符字体，属性为′FontName′，取值为：′Courier′，′宋体′，′黑体′等。 ·改变字体背景颜色，属性为′BackgroundColor′，取值为：［R,G,B］和′r′，′b′，′k′，′w′(红，蓝，黑，白)等,参见手册中的“docColorSpec”。 命令“xlabel”，“ylabel”，“title”等也可用类似方法修改字体属性，详见帮助文档。 null 修改上例为： 程序6-4 X=0:pi/100:2*pi; Y=sin(X); plot(X,Y,′--r′);%用红色虚线作图 gridon; ylabel(′y=sin2＼pix′,′FontSize′,14);%Y轴标注,设置了字号 xlabel(′x′,′FontSize′,14);%X轴标注,设置了字号 title(′functionploty=sin2＼pix′,′FontSize′,14);%图标题,设置了字号 text(0.5,sin(0.5),′＼leftarrowsin2＼pi0.5′,...%...为续行号 ′FontSize′,18,...null %字号为18号 ′BackgroundColor′,′w′);%背景为白,字符部分将盖住网格线,使字符更清晰 text(2.3,sin(2.3),′＼leftarrowsin2＼pi2.3′,′BackgroundColor′,［0.80.80.8］); %′BackgroundColor′,［0.70.70.7］使得背景为灰色 结果如图6-4所示。 null图6-4 修改作图字体属性和线型属性null 4.坐标轴的控制方法 MATLAB中对作图坐标轴的设置十分灵活，功能十分强大。用命令“docaxes”可以获得完整的帮助文档。与坐标轴设置相关的命令有axis，axex，get，set，gca等等。 下面对常用的坐标轴设置作简要介绍。 1)作图坐标范围设置(axis命令)   axis(［xminxmaxyminymax］)%用于设置x,y坐标作图范围 axisoff%用于不显示坐标 axison%用于显示坐标(默认) null 2)网格的控制(grid命令等)  gridon%用于显示网格 gridoff%用于不显示网格 set(gca,′XGrid′,′on′)%用于只显示X方向网格 set(gca,′YGrid′,′on′)%用于只显示Y方向网格 set(gca,′GridLineStyle′,′:′)%用于设置网格的线型,具体如下: %′-′表示实线;′--′表示虚线;′:′表示虚点线(默认);′-.′表示点画线 null 3)坐标轴线型的控制 set(gca,′LineWidth′,2)%控制坐标轴线宽度为2(默认为1) 4)坐标形式的控制   set(gca,′box′,′on′)%用于显示封闭形式的坐标(默认) set(gca,′box′,′off′)%用于显示开放形式的坐标 null 5)坐标刻度方向控制   set(gca,′TickDir′,′in′)%坐标刻度朝内(默认) set(gca,′TickDir′,′out′)%坐标刻度朝外 6)坐标颜色控制 set(gca,′Color′,′y′)%坐标面背景颜色设置,本例为:黄 set(gca,′XColor′,′k′)%设置横坐标轴,刻度,字符的颜色 set(gca,′YColor′,′r′)%设置纵坐标轴,刻度,字符的颜色null 7)坐标刻度字形的控制 set(gca,′FontSize′,14)%控制字体大小 set(gca,′FontWeight′,′bold′)%设置字体粗细 %有{normal}|bold|light|demi四种 8)坐标位置和方向控制 set(gca,′XAxisLocation′,′top′)%横坐标轴位于下方(bottom默认) 或上方(top)set(gca,′YAxisLocation′,′right′)%纵坐标轴位于左方(left默认)或右方(right)set(gca,′XDir′,′reverse′)%横坐标反方向(由右到左为增)set(gca,′YDir′,′reverse′)%纵坐标反方向(由右到左为增)null 9)坐标刻度线性/对数标度的设置 set(gca,′XScale′,′log′)%横坐标轴位作对数标度 set(gca,′YScale′,′log′)%纵坐标轴位作对数标度 %默认为线性标度′linear′ %用semilogx或semilogy,loglog可直接得到对数标度的作图 10)坐标刻度数的控制   set(gca,′XTick′,［］)%横坐标不标度 set(gca,′XTick′,［051629］)%在横坐标值为0,5,16,29处标度 null set(gca,′XTick′,［20:10:100］)%标度从20开始,间隔10标度,直到100 set(gca,′YTick′,［］)%对纵坐标的标度设置,同上 set(gca,′YTick′,［0:5:20］) set(gca,′XTickLabel′,{′One′;′Two′;′Three′;′Four′}) %将开始的4个刻度依次标记为字符One,Two,Three,Four %然后循环利用这4个标记将其余刻度全部标完 set(gca,′YTickLabel′,{′One′;′Two′;′Three′;′Four′}) %同上,对纵坐标作标记null ［例4］二元信号的误码概率曲线的计算为 现在要求用MATLAB作出这两条曲线，曲线使用宽度为2的粗实线，颜色为黑；手工确定作图坐标范围并手工作出坐标刻度；进行标注，注意标注中需要写入公式等特殊字符，需用LaTeX格式;网格线需要设定为细实线；横坐标表示每比特SNR，用分贝表示,纵坐标是对数刻度的。null 作图程序如下： 程序6-5  Q=inline(′0.5.*erfc(x./sqrt(2))′,′x′);%Q函数定义 gama-b-dB=0:0.5:14;%横坐标范围(分贝) gama-b=10.^(gama-b-dB./10);%横坐标范围 Pb1=Q(sqrt(gama-b));%曲线1计算 Pb2=Q(sqrt(2*gama-b));%曲线2计算 plot(gama-b-dB,Pb1,′-k′,gama-b-dB,Pb2,′-k′,′LineWidth′,2); %作图,线型为黑实线,宽度2像素null set(gca,′YScale′,′log′)%纵坐标轴位作对数标度 axis(［01410e-710e-1］);%手工设置作图范围 xlabel(′SNRperbit,＼gamma-b(dB)′,′FontSize′,12); %横轴标注,并设定标注字号 ylabel(′Probabilityoferror,P-b′,′FontSize′,12); %纵轴标注,并设定标注字号 set(gca,′GridLineStyle′,′-′)%用于设置网格的线型为实线 gridon;%开启网格线 set(gca,′MinorGridLineStyle′,′none′);%将对数分格的虚线去掉set(gca,′XTick′,［0:2:14］)%在横坐标值为0,2,4，…处标度null %下面是在图中写字 text(2.2,5e-3,′＼rho-r=-1′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(2,2e-3,′Antipodal′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(2,1e-3,′signals′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(2,0.4e-3,′P-b=Q(＼surd2＼gamma-b)′,...null ′FontSize′,12,... ′BackgroundColor′,′w′); text(10,2e-2,′＼rho-r=0′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(10,9e-3,′Orthogonal′,... ′FontSize′,12,... ′BackgroundColor′,′w′);null text(10,4e-3,′signals′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(10,1.5e-3,′P-b=Q(＼surd＼gamma-b)′,... ′FontSize′,12,... ′BackgroundColor′,′w′);   结果如图6-5所示。 null图6-5 二元信号的误码概率曲线作图实例null ［例5］学习连续信号及其采样后的离散信号的表示方法。 以取样函数f(x)=sinx/x为例，作出该函数在x∈(-10,10)内的波形。程序如下： 程序6-6  f=inline(′sin(x)./x′,′x′);%定义波形函数 x=-10:0.1:10;%x的计算范围,步进0.1 y=f(x+1e-16);%计算波形,为避免0/0,x加一微小值 plot(x,y,′--k′);%用黑色虚线作图(看一看作图结果) axis(［-1010-0.31.1］);%至此作图坐标有何变化?null holdon;%保持前图 boxoff;%坐标盒子打开(看一看坐标有何变化?) sample-time=-10:1:10;%设定离散信号的取样间隔为1 y-sample=f(sample-time+1e-16);%计算离散信号样值 h=stem(sample-time,y-sample,′fill′,′-k′); %stem的用法与plot相同,专门用于画离散信号的′火柴杆′图 %stem的用法详见docstem帮助 set(h(1),′Marker′,′o′,′MarkerSize′,5); %对离散信号作图形式进行修正,修正标记符,标记大小 null set(h(2),′LineWidth′,2); %修正到标记的竖线线型 set(h(3),′Color′,′r′); %修正基线的颜色属性 %用get(h(1)),get(h(2)),get(h(3))可以获得可修改的属性和取值 结果如图6-6所示。null图6-6 连续信号与离散信号在同一图中作出(注意boxoff的坐标形式)null 此外，我们还可以进一步对坐标轴标度进行手工设定。例如，将横坐标标度的字符进行任意设置，接上例，如果继续执行以下两句指令:   set(gca,′XTick′,［-10:2.5:0,4:4:10］)%设定标度位置 set(gca,′XTickLabel′,{′-10Ts′;′-7.5Ts′;′-5Ts′;′-2.5Ts′;... ′0′;′4Ts′;′8Ts′});%设定标度的符号 set(gca,′FontSize′,14);%设定坐标标注字号 则获得的结果如图6-7所示。 null图6-7 连续信号与离散信号在同一图中作出 (对坐标标度进行了修改)null ［例6］学习其他常用的特殊二维图形的绘制方法。利用 “bar”可以作出二维条形图，利用“stairs”可以作二维阶梯图，其用法与“plot”类似。作图程序如下： 程序6-7 t=0:1/pi:2*pi; plot(t,sin(t),′k′,′LineWidth′,2); holdon; stairs(t,sin(t),′r′);%阶梯图,注意与plot所得图的区别 %阶梯图常用来表现取样后零阶保持器的输出波形 bar(t,0.5*sin(t),′m′);%条形图,注意正弦波幅度减小了 axis(［02*pi-1.11.1］);%坐标范围 则获得的结果如图6-8所示。 null图6-8 正弦波的plot，stairs，bar作图表达的比较null ［例7］特殊函数和参数方程所表达的函数曲线的绘制问题。请在同一图上绘制以下3个方程的函数曲线： x2+3y2=5 (6-1) 和 x2+y2=1 (6-2) 以及 x=sin3tcost y=sin3tsint t∈(0,π)  (6-3)null 分析：显然，前两个方程是隐函数表示，不易得到函数的显式表达。第一个方程表示一个椭圆，而第二个方程是一个单位圆，最后一个方程则是由参数方程所表达的曲线。将前两个方程化为参数方程，即采用极坐标表示，则可以方便地进行作图。为了利用极坐标，设 x=rcost y=rsint  代入第一个方程，可得 r2(1+2sin2t)=4 (6-5)(6-4)r∈［0,+∞),t∈［0,2π)null 即 所以第一个方程对应的参数方程为 (6-6) (6-7) null 类似地推导可得，第二个方程(即单位圆)对应的参数方程为 x=cost y=sint t∈［0,2π) 下面根据这些结果编程作图。 程序6-8   t=0:0.01:2*pi；%参数t的变化范围和步进 r=sqrt(4./(1+2*(sin(t)).^2));%第一个方程的参数方程 x=r.*cos(t);y=r.*sin(t);%null plot(x,y,′-k′,′LineWidth′,3);%以3点宽黑实线对第一个方程曲线作图 holdon;%保持图不被擦除 plot(cos(t),sin(t),′-k′,′LineWidth′,1);%对单位圆作图 axis(［-33-22］);%设定作图的坐标范围 axisequal%使得纵、横坐标比例相同,保证图不失真 null %请比较一下该句执行前后的图形变化 set(gca,′FontSize′,20);%将坐标标度字体变大 set(gca,′LineWidth′,3);%将坐标轴线变粗 t=0:0.01:pi;%根据题设修改参数变化范围 x2=sin(3*t).*cos(t);y2=sin(3*t).*sin(t);%对第三个参数方程计算 plot(x2,y2);%并以默认方式作图   程序运行结果如图6-9所示。 null图6-9 隐函数和参数方程曲线的作图null 利用极坐标作图命令“polar”可以绘制极坐标表达的函数曲线。例如，绘制方程  r=0.001θ2和方程ρ=sin2tcos2t： 程序6-9 theta=0:0.1:10*pi; r=0.001*theta.^2; polar(theta,r);holdon; t=0:0.01:2*pi; polar(t,sin(2*t).*cos(2*t),′-k′); 运行程序，得到极坐标图如图6-10所示。 null图6-10 极坐标作图的例子null MATLAB提供了丰富的二维绘图指令，它们的用法与基本绘图指令“plot”类似，读者可参考在线帮助。 6.1.2 三维图形的绘制 1.三维曲线的绘制 用命令“plot3”可以进行三维空间曲线的绘制。“plot3”的常用格式是： plot3(X1,Y1,Z1,...)或 plot3(...,′PropertyName′,PropertyValue,...) 详细用法参见“docplot3”显示的帮助文档。例如，绘制三维曲线 x=(10π-t) sint y=(10π-t) cost z=t t∈(0,15π)null 的程序如下： 程序6-10   t=0:0.1:15*pi; x=(10*pi-t).*sin(t);y=(10*pi-t).*cos(t);z=t; plot3(x,y,z,′-k′,′LineWidth′,3);%作图,设定线型 gridon;%看一看,曲线像不像沙发的弹簧?   程序运行结果如图6-11所示。 null图6-11 三维曲线作图指令plot3的例子null 采用命令“stem3”可以作出三维火柴杆图。例如： 程序6-11 t=0:0.1:2*pi; x=(10*pi-t).*sin(t);y=(10*pi-t).*cos(t);z=t; stem3(x,y,z);%作图,设定线型 gridon; 运行结果如图6-12所示。 null图6-12 三维曲线作图指令stem3的例子null 2.三维曲面的绘制 常用的绘制三维曲面的命令有： (1)mesh(x,y,z)：绘制三维表面网格。 (2)meshc(x,y,z)：绘制三维表面网格,带等高线图。 (3)meshz(x,y,z)：绘制三维表面网格,给出零基准平面。 (4)surf(x,y,z)：绘制三维表面图。 (5)surfc(x,y,z)：绘制三维表面图，带等高线图。 (6)surfl(x,y,z)：绘制三维表面图，带光照效果。 null (7)contour(x,y,z,n)：绘制平面等高线图，等高线条数为n。 (8)contour3(x,y,z,n)：在三维空间绘制等高线图，等高线条数为n。 (9)waterfall(x,y,z)：绘制三维瀑布图。 这些命令的详细用法请参考帮助文档。下面举例加以说明。 考虑绘制一个二元函数 z=f(x,y)=(x2-2x) e-x2-y2-xy (6-10) null 其中，绘制范围为x∈(-3,3); y∈(-2,2)。首先用“meshgrid”函数产生一个x和y的网格矩阵，即产生一个x轴坐标起始于-3，终止于3，步进为1, y轴坐标起始于-2，终止于2，步进为1的网格分割。程序如下： 程序6-12 ［x,y］=meshgrid(-3:0.1:3,-2:0.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); mesh(x,y,z); 绘图结果如图6-13所示。 null图6-13 三维曲面作图指令mesh的例子null 如果使用指令“meshc、meshz”代替“mesh”指令，即   figure(1);meshc(x,y,z);figure(2);meshz(x,y,z);   则分别作出带等高线图的以及给出零基准平面的三维网格图，如图6-14(a)、(b)所示。 null 图6-14带等高线图的以及给出零基准平面的三维网格图 (a）三维曲面作图指令meshc的例子； (b）三维曲面作图指令meshz的例子 null 如果使用指令“surf、surfc、surfl”代替“mesh”指令，即 figure(1);surf(x,y,z);figure(2);surfc(x,y,z);figure(3);surfl(x,y,z); 则分别作出三维表面图、带等高线的三维表面图以及具有光照效果的三维表面图，如图6-15 (a)、(b)、(c)所示。 null图6-15三维表面图、带等高线的以及具有光照效果的三维表面图 (a）指令surf的例子；(b）指令surfc的例子；(c）指令surfl的例子 null 另外，使用“contour”指令可以绘制平面等高线图，使用“contour3”指令能够在三维空间绘制出等高线图来，而使用“waterfall”则可绘制出三维瀑布图，即 figure(1);contour(x,y,z,30);figure(2);contour3(x,y,z,30);figure(3);waterfall(x,y,z);   绘图结果如图6-16所示。 null图6-16平面等高线图、三维等高线图以及三维瀑布图 (a）指令contour的例子；(b）指令contour3的例子； (c）指令waterfall的例子null 6.1.3 子图的绘制 我们可以用“subplot”命令将整个图形窗口分割成若干子图，每个子图可以用不同绘图语句独立绘制出所需图形。指令“subplot(m，n，p)”用于创建m行×n列子图中的第p个子图,现举例说明。以下程序的作图结果如图6-17所示。 程序6-13 t=0:pi/10:4*pi; subplot(2,2,1);plot(t,sin(t));title(′subfig1′) subplot(2,2,2);plot(t,sin(2*t));title(′subfig2′) null subplot(2,2,3);stem(t,sin(t));title(′subfig3′) ［x,y］=meshgrid(-3:0.1:3,-2:0.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); subplot(2,2,4);mesh(x,y,z);title(′subfig4′) null图6-17 用subplot指令创建的子图null 6.1.4 创建图形窗口 如果在一个程序中需要创建多个图形窗口来作多幅图，那么就需要创建不同的图形窗口。创建图形窗口的指令是“figure(n)”。例如，下面程序创建了三个图形窗口，其中第三个图形窗口中又分割为上、下两个子图，如图6-18所示。 null图6-18 用figure指令创建多个作图窗口null 程序6-14 t=0:0.5:6*pi; figure(1);%创建图1 plot(t,t.*sin(t)); figure(2);%创建图2 stem(t,cos(t)); figure(3);%创建图3 subplot(2,1,1);plot(t,cos(t));%图3中上图 subplot(2,1,2);stairs(t,cos(t));%图3中下图 null 6.1.5 图像表现 MATLAB可以通过“imread”命令从图像文件中读取图像数据矩阵，并通过“imshow”命令显示出来，还可以由“imwrite”将矩阵数据写为图像格式，为图像信号处理提供方便。在Toolbox＼images＼imdemos子下有多幅图片和图像处理的例子供测试用，例如图片cameraman.tif等。如下程序将演示如何读取灰度图，如何进行对比度和亮度设置、图像二值化、灰度直方图绘制以及格式转换、存盘等等，运行结果如图6-19所示。 null图6-19 图像处理的输出结果null 程序6-15  I=imread(′cameraman.tif′);%读取tif格式的图片256×256，灰度图 I=double(I);%将图像数据由uint8类型转换为double类型，便于数值运算 Ia=fliplr(I);%图像左右翻转(上下翻转可用flipud） Ib=70+I*0.6;%演示提高亮度，减小对比度 Ic=255*(I>128);%演示门限为128的二值化 I=uint8(I);%将计算结果由double类型转换为uint8类型，以便显示 null Ia=uint8(Ia); Ib=uint8(Ib); Ic=uint8(Ic); subplot(2,2,1);imshow(I);%显示原图片 subplot(2,2,2);imshow(Ia); subplot(2,2,3);imshow(Ib); subplot(2,2,4);imshow(Ic); imwrite(I,′cameraman.jpg′,′jpg′); %将图像矩阵数据存为jpg格式，当前工作路径下名为cameraman.jpg null 6.1.6 动画 利用MATLAB中的动态图形指令可以方便地得到图形或曲线的动态变化效果。 1.影片动画 使用指令“moviein、getframe”可以制作动画，然后用指令“movie”播出动画。首先，使用“moviein(n)”指令预先分配一个n列矩阵，每列将存储一帧动画数据；然后,使用“getframe”指令将当前图形窗口的像素以列的形式存储到分配的矩阵中；最后,通过“movie” 指令按照一定的速度播放这些存储的帧，于是就形成了动画。 null 例如： 程序6-16 M=moviein(15);%预先分配一个能够存储15帧的矩阵  Z=peaks;surf(Z);%第0帧作图，并通过以下两句来固定作图坐标系 axistight;%或使用axismanual set(gca,′nextplot′,′replacechildren′);%保持作图的帧使用相同坐标系 fork=1:15%循环，形成15帧图形数据并记录null surf(sin(2*pi*k/20)*Z,Z);%作图 M(k)=getframe;%录制作图的每一帧，每次循环得到一帧 end movie(M,20,10)%反复播放20次，播放速度每秒10帧 使用指令“movie2avi”可以将动画帧矩阵保存为avi格式的影片。另外，用“aviread(′filename′)”可以将avi影片文件读入，然后使用“movie”指令播放。例如：   %接上例 movie2avi(M,′myavi′);%将动画帧矩阵M保存为avi格式的影片myavi.avinull clear;%清空内存变量 F=aviread(′myavi′);%读入avi文件数据 movie(F,10,5);%反复播放10次，播放速度每秒5帧 2.彗星轨线 采用彗星轨线指令“comet”、“comet3”可以分别在二维、三维上绘制出质点的运动轨迹。其常用用法是： comet(x,y,p) comet3(x,y,z,p) 其中，x，y，z分别是轨迹的坐标点序列，与“plot”、“plot3”指令的要求相同。 参数p是可选项，默认值为0.1，代表作图的慧尾长度参数。p可取0到1之间的数。  null 例如： 程序6-17 t=-10*pi:pi/250:10*pi; comet3((cos(2*t).^2).*sin(t),(sin(2*t).^2).*cos(t),t); 得到的动画绘图结果如图6-20所示。(动画的速度与计算机配置有关。） null图6-20 用comet3指令动画绘制三维曲线的瞬间6.2 MATLAB仿真中的信号处理6.2 MATLAB仿真中的信号处理 6.2.1 数据采集 1.MATLAB中信号的表示 连续信号通常指的是时间连续、幅度连续的信号，也称为模拟信号，在数学上表示为一个时间连续函数f(t) 。数字信号是指时间离散而且幅度也离散的信号，在数学上可以表示为有限精度的一个离散时间序列{x1,x2,…,xn,…}。 对模拟信号进行时间域取样，就获得了时间域离散的取样信号。若取样速率满足奈奎斯特取样定理，即取样速率大于等于模拟信号的最高频率的2倍，那么模拟信号可以由其取样值序列构成的时间离散信号无失真地表达。null 对时间离散信号再进行量化，就获得了幅度离散的样值序列，也就是数字信号。量化将引入量化误差，形成量化噪声成分。当量化足够精细时，量化噪声可以忽略。这时我们就认为，所获得的数字信号就无失真地代表了原来的模拟信号。 在数字计算机中所能够存储的是数字序列，也即模拟信号必须通过取样和量化后，变成相应的数字信号，才能被计算机存储和处理。最后，可以将处理结果序列以取样速率送入理想低通滤波器，恢复为所需的模拟信号。 对于音频信号来说，实现模拟音频信号与数字音频信号之间的转化的模块就是声卡，MATLAB可以方便地对声卡进行诸如采样率等输入/输出参数的配置，从而实现录音和放音，而且其声卡操作函数是跨操作平台的，也是硬件无关的。 null 2.音频信号的获取与输出 对于配置了声卡的多媒体计算机，MATLAB中可以采用命令“wavrecord”来录音。“wavrecord”的调用格式是: y=wavrecord(n,Fs,ch,dtype); %n为总的取样点数，Fs为取样速率(样点/秒) %标准取样速率可设为8000，11025(默认)，2250以及44100(样点/秒) %用户也可以设定其他取样速率值，如Fs=10000等，但需满足取样定理要求 %ch为录音声道数，默认ch=1，为单声道录音，若ch=2，为立体声录音 %dtype为记录的数据格式，有double(默认)，single，int16，int8类型null %详细用法参见“docwavrecord”帮助 %常用用法是: y=wavrecord(n，Fs);%回车后立即开始录音，直到录音结束才返回提示符 %该语句以双精度记录n点采样速率为Fs的单声道数字音频。y的值域在-1到1之间特别需要强调的是，录音输出序列y已经是一个n×1的数字信号序列。采用均匀量化规则，对于“double”(默认)，“single”，“int16”数据类型，每个样值的量化精度为16bit，这对工程应用来说是足够的。一般情况下，我们可以认为这样的量化过程的量化噪声几乎可以忽略。 当我们需要研究8bitA律PCM的语音质量时，就可以将上述的录音结果视为量化之前的取样样值序列。 null 使用命令“sound”将数字序列以设定的采样速率输出到声卡，通过声卡转化为模拟音频信号。“sound”的用法是:   sound(y，Fs); %y为取值在［-1，1］区间的n行1列的数字序列(单声道输出) %对于立体声输出，y为取值在［-1，1］区间的n行2列的数字序列 %Fs为设定的取样速率，一般声卡支持Fs的范围是5000(Hz)到441000(Hz)   MATLAB也可以将记录的音频信号直接保存为“*.wav”格式或“*.au”格式。在Windows环境中，“*.wav”格式是最为常用的；在UNIX/Linux系统中，“*.au”格式是通行的。null 利用命令“wavwrite(y，Fs，′filename′);”就可以将向量y存储为取样率为Fs的wav音频文件。wav文件可以采用Windows中的多媒体播放器进行播放。如果要保存为“*.au”格式，则采用命令“auwrite(y，Fs，′aufile′);”即可。注意，数字音频信号的向量y是n行1列的(对于单声道而言)或n行2列的(对于双声道立体声而言)。另外，y的取值范围要在［-1,1］区间。 null MATLAB可以直接读取“*.wav”格式和“*.au”格式的音频文件，其常用的调用命令分别是:   ［y,Fs］=wavread(′filename′);%读取*.wav文件,返回y为样值序列,Fs为采样率 ［y,Fs］=wavread(′filename′,［N1N2］); %读取*.wav文件中从第N1采样点到N2采样点之间的样值序列 %返回y为样值序列,Fs为采样率 ［y,Fs］=auread(′aufile′);%读取*.au文件 ［y,Fs］=auread(′aufile′,［N1,N2］); %读取*.au文件中从第N1采样点到N2采样点之间的样值序列 %返回y为样值序列,Fs为采样率 null 3.图像信号和视频的获取 MATLAB可以从各种图像格式的文件中读取图像的点阵数据，也可以从avi格式的视频文件中获取数据，它们分别使用的命令是“imread”、“aviread”等。也可以将处理的图像数据保存为各种标准格式的图片或avi格式的视频文件，参见上节叙述。 6.2.2 确定信号的产生 1.正弦波和余弦波的产生 确定周期信号的产生一般需要考虑以下问题:首先，要设置取样频率，将模拟信号离散化，取样频率的设置要满足取样定理，为了信号波形显示光滑好看，一般取样频率设置为信号最高频率的5～20倍。然后，确定所要产生的信号的持续时间长度，再确定信号的其他参数，如幅度、相位等。最后，画出信号波形图，验证正确性。如果信号频率为音频(20Hz～20kHz)，可以利用“sound”函数从声卡输出信号波形。 null 利用函数sin和cos可以产生所需要的正弦波或余弦波。例如，要产生一个频率为150Hz，幅度为0.45，初始相位为35°的正弦波，信号持续时间为5s,作出波形图并从声卡输出。根据取样定理，采样率至少为信号最高频率的两倍，因此，本例中采样率必须大于300次/s。 为使得产生波形光滑，这里取采样率为2000次/s,信号持续时间为5s，则总的取样点数为点5×2000=10000。下面是相应的波形产生程序。 null 程序6-18 Fs=2000;%采样率 t=1/Fs:1/Fs:5;%信号持续时间范围 f=150;%正弦波频率 A=0.45;%幅度 Fai=35/180*pi;%相位 X=A*sin(2*pi*f.*t+Fai);%正弦波计算 plot(t(1:100),X(1:100));%作出波形图前100样点 xlabel(′time(sec)′); null ylabel(′sin2＼pift′); title(′150Hzsinwave′); disp(′按任意键开始播出5秒的150Hz正弦波...′);pause; sound(X′,Fs);%从声卡播放150Hz的正弦波 disp(′播放结束,下面将音频数据存盘为C:＼my50HzSIN.wav′); wavwrite(X,Fs,′C:＼my50HzSIN.wav′);%音频数据存盘 clear;%清除内存变量 ［R,Rs］=wavread(′C:＼mySIN.wav′);%读取刚刚存盘的音频数据 sound(R,Rs);%播放， %读者可以试一试用Windows媒体播放器来播放mySIN.wav文件 null 2.周期方波和周期三角波 利用MATLAB的信号处理工具箱中的“square”命令可以产生方波，而命令“sawtooth”可以产生三角波(锯齿波)。“square”命令的调用格式是:  x=square(2*pi*f*t+fai,duty); %t时间取样序列 %f方波的基频率 %fai方波的初相位 %duty占空比,为0～100之间的数 null %如占空比为20，表示在方波的一个周期内正电平占20% %x返回的幅度为1的矩形波样值序列   “sawtooth”命令的调用格式是:   x=sawtooth(2*pi*f*t+fai,width); %t时间取样序列 %f三角波的基频率 %fai三角波的初相位 %width宽度,为0～1之间的数 %如为0.3,表示在三角波的一个周期内上升沿占30% %x返回的幅度为1的锯齿波波样值序列 null 3.任意确定周期信号的产生 给定周期信号在一个周期内的波形样点值，就可以通过样值的周期重复来获得任意确定周期信号。设在一个周期内，信号波形的样值点序列为(列向量) (6-11) null 欲将之扩展为3个周期，可首先进行矩阵乘法运算: 然后利用MATLAB的“reshape”命令将结果矩阵按照列顺序变形为3n行1列的向量,即(6-12) (6-13) null 例如，给定的某心电图一个周期的波形数据为n(共37点采样点，采样间隔为0.02秒)。 下面程序作出的采样率为200和8000时的心电波形图如图6-21所示。 程序6-19 n=［ 78.0 -39.5 -70.7 -17.0 -17.0 -16.5 -16.5... -15.0 -15.0 -14.5 -10.0 -3.5 1.0 -3.5... -12.0 -16.0 -17.0 -18.0 -18.0 -18.0 -18.0... -18.0 -18.0 -18.0 -18.0 -16.5 -4.0 9.0... null 19.0 6.5 -11.0 -17.0 -17.0 -18.0 -23.0... -18.3 49.0］; %一个心跳周期的数据 Ts=0.02;%采样时间间隔 Fs=1/Ts;%采样率 m=10;%现在将该信号扩展为m个周期 nm=reshape(n′*ones(1,m),1,m*length(n)); nm=nm/max(nm);%归一化 t=Ts:Ts:length(nm)*Ts;%采样时间刻度 subplot(2,1,1);plot(t,nm,′.k′);axis(［02-22］);%波形图 %下面为了“听”心跳,需要采用插值方法将波形采样率提高到声卡所能够允许的值,如8000次/snull t8000=［Ts:1/8000:t(length(t))］;%新的采样时间刻度 n8000=interp1(t,nm,t8000,′spline′); subplot(2,1,2);plot(t8000,n8000,′.k′);axis(［02-22］);%波形图 sound(n8000,8000);%听一听心跳,注意,由于心跳频率很低,需将声音放大了仔细听! null图6-21 心电图插值前后的采样点比较null 4.脉冲信号的产生 MATLAB的信号处理工具箱还提供了各种脉冲信号的发生函数。当然，我们不利用这些函数，根据各种脉冲的数学表达式也能够很方便地计算出它们的波形。这些脉冲信号发生函数的用法详见帮助文档：   y=rectpuls(t,w);%产生矩形脉冲 y=tripuls(T,w,s)%产生三角形脉冲 y=diric(x,n)%产生狄拉克函数脉冲(冲激脉冲) yi=gauspuls(t,fc,bw,bwr);%产生高斯调制的正弦脉冲 y=pulstran(t,d,′func′);%产生脉冲串 y=sinc(x)%产生取样函数sinc脉冲 null 5.扫频信号的产生 利用MATLAB的信号处理工具箱中的“chirp”函数可以获得在设定频率范围内按照设定方式进行的扫频信号。其调用方式是： y=chirp(t,f0,t1,f1) y=chirp(t,f0,t1,f1,′method′) y=chirp(t,f0,t1,f1,′method′,phi) y=chirp(t,f0,t1,f1,′quadratic′,phi,′shape′)   详细用法参见帮助文档。 null 6.2.3 随机信号的产生 1.均匀分布的随机数序列的产生 命令“rand”可以产生在［0,1］区间上均匀分布的随机数序列。“rand”的一般调用格式是：   Y=rand(m,n);%产生m行n列的随机数矩阵   例如，产生一个在［-5,1］区间均匀分布的随机数矩阵，共5行5列。 null >>A=rand(5,5)*6-5 A= -2.1759-3.4137-4.7372-3.7770-1.9656 -0.0592-3.7619-4.9632-3.5947-4.6341 0.09990.5933-0.09730.3910-4.1883 0.7854-2.64000.4911-0.2735-4.8665 -0.1239-4.6258-0.5502-3.90210.2189 >>A=rand(1,100000)*6-5;%产生10万个随机数 >>max(A)%最大值 ans= 0.9999null >>min(A)%最小值,检验随机数值范围 ans= -4.9999 >>x=-7:0.1:3;hist(A,x);%作出A中随机数的发生频率分布 A中随机数的发生频率分布图如图6-22所示。 null图6-22 均匀分布的随机数的发生频率分布图null 又例如，已知一个二元离散信源模型为 (6-14) 现用MATLAB产生符合该信源统计特性的1000个输出符号，程序r=rand(1,1000);%产生1000个［0,1］区间均匀分布的随机数out=(r>0.8);%使得r中大于0.8的值出现时(概率为0.2)输出为1%否则输出为0 。 null 2.高斯分布的随机数序列的产生 利用命令“randn”可以产生零均值方差为1的标准正态分布的随机数序列。其用法是：   Y=randn(m,n);%产生m行n列的零均值方差为1的标准正态分布的随机数矩阵   例如，需要获得10000个均值为60、方差为3.7的正态分布随机数序列，可用下面的程序实现。 程序6-20 M=60;%设定均值 D=3.7;%设定方差null Y=M+sqrt(D)*randn(1,10000); M1=mean(Y)%检验均值是否符合设定要