第6章 信号处理

第6章信号处理6.1MATLAB的数据可视化表达6.2MATLAB仿真中的信号处理6.3信号特征参数的计算6.4信号的频谱分析6.5线性系统6.6数字滤波器设计与实现6.1MATLAB的数据可视化表达 6.1.1二维图形的绘制 1.常用的二维图形绘图函数 基本的二维绘图函数有： ·plot：绘制二维曲线； ·title：给图形加标题； ·grid：显示网格线； ·xlabel：给x轴加标记； ·ylabel：给y轴加标记； ·text：在坐标图中加入文字注释。 ［例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轴标注，可以有汉字 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所示。图6-1基本的二维绘图函数用法 2.图形的线型和颜色控制 在命令“plot”的高级用法中，可以设置作图的线型，标记类型，线和标记的颜色、粗细等特征。用命令“docLineSpec”和“docplot”可以查询详细的帮助文档。作图时常用的线型、标记以及颜色的定义参见表6-1。表6-1作图时常用的线型、标记以及颜色的定义 利用线型属性和标记属性可以随心所欲地设计作图图式。常用的线型属性和标记属性参见表6-2。表6-2常用的线型属性和标记属性 程序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所示(黑白印刷时颜色表现不出来)。图6-2二维绘图线型和标记的特征控制(1) ［例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所示。图6-3二维绘图线型和标记的特征控制(2) 3.图形的标注 图形的标注可以用“text”函数。其用法是： text(x,y,′字符串′); text(...′PropertyName′,PropertyValue...);

其中，字符串中若有特殊符号，如希腊字母、箭头等，需要采用LaTeX格式表示。′PropertyName′为字符属性名称，PropertyValue为相应的属性取值。用“doctext-props” 可以得到详细的用法手册。对于常用的一些，举例如下： ·改变字符大小，属性为′FontSize′，取值为：10，12，16等。 ·改变字符字体，属性为′FontName′，取值为：′Courier′，′宋体′，′黑体′等。 ·改变字体背景颜色，属性为′BackgroundColor′，取值为：［R,G,B］和′r′，′b′，′k′，′w′(红，蓝，黑，白)等,参见手册中的“docColorSpec”。 命令“xlabel”，“ylabel”，“title”等也可用类似方法修改字体属性，详见帮助文档。

修改上例为： 程序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,... %字号为18号 ′BackgroundColor′,′w′);%背景为白,字符部分将盖住网格线,使字符更清晰 text(2.3,sin(2.3),′＼leftarrowsin2＼pi2.3′,′BackgroundColor′,［0.80.80.8］); %′BackgroundColor′,［0.70.70.7］使得背景为灰色 结果如图6-4所示。图6-4修改作图字体属性和线型属性 4.坐标轴的控制方法 MATLAB中对作图坐标轴的设置十分灵活，功能十分强大。用命令“docaxes”可以获得完整的帮助文档。与坐标轴设置相关的命令有axis，axex，get，set，gca等等。 下面对常用的坐标轴设置作简要介绍。 1)作图坐标范围设置(axis命令)

axis(［xminxmaxyminymax］)%用于设置x,y坐标作图范围 axisoff%用于不显示坐标 axison%用于显示坐标(默认) 2)网格的控制(grid命令等)

gridon%用于显示网格 gridoff%用于不显示网格 set(gca,′XGrid′,′on′)%用于只显示X方向网格 set(gca,′YGrid′,′on′)%用于只显示Y方向网格 set(gca,′GridLineStyle′,′:′)%用于设置网格的线型,具体如下: %′-′表示实线;′--′表示虚线;′:′表示虚点线(默认);′-.′表示点画线 3)坐标轴线型的控制 set(gca,′LineWidth′,2)%控制坐标轴线宽度为2(默认为1) 4)坐标形式的控制

set(gca,′box′,′on′)%用于显示封闭形式的坐标(默认) set(gca,′box′,′off′)%用于显示开放形式的坐标 5)坐标刻度方向控制

set(gca,′TickDir′,′in′)%坐标刻度朝内(默认) set(gca,′TickDir′,′out′)%坐标刻度朝外 6)坐标颜色控制 set(gca,′Color′,′y′)%坐标面背景颜色设置,本例为:黄 set(gca,′XColor′,′k′)%设置横坐标轴,刻度,字符的颜色 set(gca,′YColor′,′r′)%设置纵坐标轴,刻度,字符的颜色

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′)%纵坐标反方向(由右到左为增) 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处标度

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′}) %同上,对纵坐标作标记 ［例4］二元信号的误码概率曲线的计算公式为

现在要求用MATLAB作出这两条曲线，曲线使用宽度为2的粗实线，颜色为黑；手工确定作图坐标范围并手工作出坐标刻度；进行标注，注意标注中需要写入公式等特殊字符，需用LaTeX格式;网格线需要设定为细实线；横坐标表示每比特SNR，用分贝表示,纵坐标是对数刻度的。 作图程序如下： 程序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像素

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，…处标度

%下面是在图中写字 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)′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(10,2e-2,′＼rho-r=0′,... ′FontSize′,12,... ′BackgroundColor′,′w′); text(10,9e-3,′Orthogonal′,... ′FontSize′,12,... ′BackgroundColor′,′w′); 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所示。图6-5二元信号的误码概率曲线作图实例 ［例5］学习连续信号及其采样后的离散信号的表示方法。 以取样函数f(x)=sinx/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］);%至此作图坐标有何变化? 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); %对离散信号作图形式进行修正,修正标记符,标记大小 set(h(2),′LineWidth′,2); %修正到标记的竖线线型 set(h(3),′Color′,′r′); %修正基线的颜色属性 %用get(h(1)),get(h(2)),get(h(3))可以获得可修改的属性和取值 结果如图6-6所示。图6-6连续信号与离散信号在同一图中作出(注意boxoff的坐标形式) 此外，我们还可以进一步对坐标轴标度进行手工设定。例如，将横坐标标度的字符进行任意设置，接上例，如果继续执行以下两句指令:

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所示。图6-7连续信号与离散信号在同一图中作出(对坐标标度进行了修改)

［例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所示。图6-8正弦波的plot，stairs，bar作图表达的比较 ［例7］特殊数学函数和参数方程所表达的函数曲线的绘制问题。请在同一图上绘制以下3个方程的函数曲线：

x2+3y2=5(6-1) 和

x2+y2=1(6-2) 以及

x=sin3tcost

y=sin3tsintt∈(0,π)

(6-3) 分析：显然，前两个方程是隐函数表示，不易得到函数的显式表达。第一个方程表示一个椭圆，而第二个方程是一个单位圆，最后一个方程则是由参数方程所表达的曲线。将前两个方程化为参数方程，即采用极坐标表示，则可以方便地进行作图。为了利用极坐标，设

x=rcost

y=rsint 代入第一个方程，可得

r2(1+2sin2t)=4(6-5)(6-4)r∈［0,+∞),t∈［0,2π) 即 所以第一个方程对应的参数方程为

(6-6)(6-7)

类似地推导可得，第二个方程(即单位圆)对应的参数方程为

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);% plot(x,y,′-k′,′LineWidth′,3);%以3点宽黑实线对第一个方程曲线作图 holdon;%保持图不被擦除 plot(cos(t),sin(t),′-k′,′LineWidth′,1);%对单位圆作图 axis(［-33-22］);%设定作图的坐标范围 axisequal%使得纵、横坐标比例相同,保证图不失真 %请比较一下该句执行前后的图形变化 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所示。图6-9隐函数和参数方程曲线的作图 利用极坐标作图命令“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所示。图6-10极坐标作图的例子

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π) 的程序如下： 程序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所示。图6-11三维曲线作图指令plot3的例子 采用命令“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所示。图6-12三维曲线作图指令stem3的例子 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)：绘制三维表面图，带光照效果。 (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) 其中，绘制范围为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所示。图6-13三维曲面作图指令mesh的例子 如果使用指令“meshc、meshz”代替“mesh”指令，即

figure(1);meshc(x,y,z);figure(2);meshz(x,y,z);

则分别作出带等高线图的以及给出零基准平面的三维网格图，如图6-14(a)、(b)所示。 图6-14带等高线图的以及给出零基准平面的三维网格图 (a）三维曲面作图指令meshc的例子； (b）三维曲面作图指令meshz的例子 如果使用指令“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)所示。图6-15三维表面图、带等高线的以及具有光照效果的三维表面图(a）指令surf的例子；(b）指令surfc的例子；(c）指令surfl的例子 另外，使用“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所示。图6-16平面等高线图、三维等高线图以及三维瀑布图(a）指令contour的例子；(b）指令contour3的例子；(c）指令waterfall的例子 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′) 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′)图6-17用subplot指令创建的子图 6.1.4创建图形窗口 如果在一个程序中需要创建多个图形窗口来作多幅图，那么就需要创建不同的图形窗口。创建图形窗口的指令是“figure(n)”。例如，下面程序创建了三个图形窗口，其中第三个图形窗口中又分割为上、下两个子图，如图6-18所示。图6-18用figure指令创建多个作图窗口 程序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中下图 6.1.5图像表现 MATLAB可以通过“imread”命令从图像文件中读取图像数据矩阵，并通过“imshow”命令显示出来，还可以由“imwrite”将矩阵数据写为图像格式，为图像信号处理提供方便。在Toolbox＼images＼imdemos子目录下有多幅图片和图像处理的例子供测试用，例如图片cameraman.tif等。如下程序将演示如何读取灰度图，如何进行对比度和亮度设置、图像二值化、灰度直方图绘制以及格式转换、存盘等等，运行结果如图6-19所示。图6-19图像处理的输出结果 程序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类型，以便显示 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 6.1.6动画 利用MATLAB中的动态图形指令可以方便地得到图形或曲线的动态变化效果。 1.影片动画 使用指令“moviein、getframe”可以制作动画，然后用指令“movie”播出动画。首先，使用“moviein(n)”指令预先分配一个n列矩阵，每列将存储一帧动画数据；然后,使用“getframe”指令将当前图形窗口的像素以列的形式存储到分配的矩阵中；最后,通过“movie”指令按照一定的速度播放这些存储的帧，于是就形成了动画。 例如： 程序6-16 M=moviein(15);%预先分配一个能够存储15帧的矩阵 Z=peaks;surf(Z);%第0帧作图，并通过以下两句来固定作图坐标系 axistight;%或使用axismanual set(gca,′nextplot′,′replacechildren′);%保持作图的帧使用相同坐标系 fork=1:15%循环，形成15帧图形数据并记录 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.avi

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之间的数。

例如： 程序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所示。(动画的速度与计算机配置有关。）图6-20用comet3指令动画绘制三维曲线的瞬间6.2MATLAB仿真中的信号处理 6.2.1数据采集 1.MATLAB中信号的表示 连续信号通常指的是时间连续、幅度连续的信号，也称为模拟信号，在数学上表示为一个时间连续函数f(t)。数字信号是指时间离散而且幅度也离散的信号，在数学上可以表示为有限精度的一个离散时间序列{x1,x2,…,xn,…}。 对模拟信号进行时间域取样，就获得了时间域离散的取样信号。若取样速率满足奈奎斯特取样定理，即取样速率大于等于模拟信号的最高频率的2倍，那么模拟信号可以由其取样值序列构成的时间离散信号无失真地表达。 对时间离散信号再进行量化，就获得了幅度离散的样值序列，也就是数字信号。量化将引入量化误差，形成量化噪声成分。当量化足够精细时，量化噪声可以忽略。这时我们就认为，所获得的数字信号就无失真地代表了原来的模拟信号。 在数字计算机中所能够存储的是数字序列，也即模拟信号必须通过取样和量化后，变成相应的数字信号，才能被计算机存储和处理。最后，可以将处理结果序列以取样速率送入理想低通滤波器，恢复为所需的模拟信号。 对于音频信号来说，实现模拟音频信号与数字音频信号之间的转化的模块就是声卡，MATLAB可以方便地对声卡进行诸如采样率等输入/输出参数的配置，从而实现录音和放音，而且其声卡操作函数是跨操作平台的，也是硬件无关的。

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类型 %详细用法参见“docwavrecord”帮助 %常用用法是: y=wavrecord(n，Fs);%回车后立即开始录音，直到录音结束才返回提示符 %该语句以双精度记录n点采样速率为Fs的单声道数字音频。y的值域在-1到1之间特别需要强调的是，录音输出序列y已经是一个n×1的数字信号序列。采用均匀量化规则，对于“double”(默认)，“single”，“int16”数据类型，每个样值的量化精度为16bit，这对工程应用来说是足够的。一般情况下，我们可以认为这样的量化过程的量化噪声几乎可以忽略。 当我们需要研究8bitA律PCM的语音质量时，就可以将上述的录音结果视为量化之前的取样样值序列。 使用命令“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”格式是通行的。 利用命令“wavwrite(y，Fs，′filename′);”就可以将向量y存储为取样率为Fs的wav音频文件。wav文件可以采用Windows中的多媒体播放器进行播放。如果要保存为“*.au”格式，则采用命令“auwrite(y，Fs，′aufile′);”即可。注意，数字音频信号的向量y是n行1列的(对于单声道而言)或n行2列的(对于双声道立体声而言)。另外，y的取值范围要在［-1,1］区间。

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为采样率

3.图像信号和视频的获取 MATLAB可以从各种图像格式的文件中读取图像的点阵数据，也可以从avi格式的视频文件中获取数据，它们分别使用的命令是“imread”、“aviread”等。也可以将处理的图像数据保存为各种标准格式的图片或avi格式的视频文件，参见上节叙述。 6.2.2确定信号的产生 1.正弦波和余弦波的产生 确定周期信号的产生一般需要考虑以下问题:首先，要设置取样频率，将模拟信号离散化，取样频率的设置要满足取样定理，为了信号波形显示光滑好看，一般取样频率设置为信号最高频率的5～20倍。然后，确定所要产生的信号的持续时间长度，再确定信号的其他参数，如幅度、相位等。最后，画出信号波形图，验证正确性。如果信号频率为音频(20Hz～20kHz)，可以利用“sound”函数从声卡输出信号波形。 利用函数sin和cos可以产生所需要的正弦波或余弦波。例如，要产生一个频率为150Hz，幅度为0.45，初始相位为35°的正弦波，信号持续时间为5s,作出波形图并从声卡输出。根据取样定理，采样率至少为信号最高频率的两倍，因此，本例中采样率必须大于300次/s。 为使得产生波形光滑，这里取采样率为2000次/s,信号持续时间为5s，则总的取样点数为点5×2000=10000。下面是相应的波形产生程序。 程序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)′);

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文件 2.周期方波和周期三角波 利用MATLAB的信号处理工具箱中的“square”命令可以产生方波，而命令“sawtooth”可以产生三角波(锯齿波)。“square”命令的调用格式是:

x=square(2*pi*f*t+fai,duty); %t时间取样序列 %f方波的基频率 %fai方波的初相位 %duty占空比,为0～100之间的数

%如占空比为20，表示在方波的一个周期内正电平占20% %x返回的幅度为1的矩形波样值序列

“sawtooth”命令的调用格式是:

x=sawtooth(2*pi*f*t+fai,width); %t时间取样序列 %f三角波的基频率 %fai三角波的初相位 %width宽度,为0～1之间的数 %如为0.3,表示在三角波的一个周期内上升沿占30% %x返回的幅度为1的锯齿波波样值序列 3.任意确定周期信号的产生 给定周期信号在一个周期内的波形样点值，就可以通过样值的周期重复来获得任意确定周期信号。设在一个周期内，信号波形的样值点序列为(列向量) (6-11) 欲将之扩展为3个周期，可首先进行矩阵乘法运算: 然后利用MATLAB的“reshape”命令将结果矩阵按照列顺序变形为3n行1列的向量,即(6-12)(6-13) 例如，给定的某心电图一个周期的波形数据为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.51.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.09.0...

19.06.5-11.0-17.0-17.0-18.0-23.0... -18.349.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 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);%听一听心跳,注意,由于心跳频率很低,需将声音放大了仔细听!图6-21心电图插值前后的采样点比较

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脉冲

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′)

详细用法参见帮助文档。 6.2.3随机信号的产生 1.均匀分布的随机数序列的产生 命令“rand”可以产生在［0,1］区间上均匀分布的随机数序列。“rand”的一般调用格式是：

Y=rand(m,n);%产生m行n列的随机数矩阵

例如，产生一个在［-5,1］区间均匀分布的随机数矩阵，共5行5列。

>>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 >>min(A)%最小值,检验随机数值范围 ans= -4.9999 >>x=-7:0.1:3;hist(A,x);%作出A中随机数的发生频率分布 A中随机数的发生频率分布图如图6-22所示。图6-22均匀分布的随机数的发生频率分布图 又例如，已知一个二元离散信源模型为

(6-14)

现用MATLAB产生符合该信源统计特性的1000个输出符号，程序r=rand(1,1000);%产生1000个［0,1］区间均匀分布的随机数out=(r>0.8);%使得r中大于0.8的值出现时(概率为0.2)输出为1%否则输出为0。 2.高斯分布的随机数序列的产生 利用命令“randn”可以产生零均值方差为1的标准正态分布的随机数序列。其用法是：

Y=randn(m,n);%产生m行n列的零均值方差为1的标准正态分布的随机数矩阵

例如，需要获得10000个均值为60、方差为3.7的正态分布随机数序列，可用下面的程序实现。 程序6-20 M=60;%设定均值 D=3.7;%设定方差 Y=M+sqrt(D)*randn(1,10000); M1=mean(Y)%检验均值是否符合设定要求 D1=var(Y)%检验方差是否符合设定要求 hist(Y,［45:0.1:75］);%检验随机数分布情况 M1= 60.0128 D1= 3.7269 Y中随机数的发生频率分布图如图6-23所示。图6-23高斯分布的随机数的发生频率分布图 3.其他分布随机数序列的产生 在MATLAB的统计工具箱中，提供了几乎所有概率分布形式的随机数矩阵的产生方法。其调用方法是：

y=random(′name′,A1,A2,A3,m,n); 其中，′name′为分布名称，含义参见表6-3。A1，A2，A3分别是相应分布的参数，m，n为生成随机数矩阵的行数和列数。更详细用法请参考帮助文档。表6-3random函数中参数′name′的含义 另外，在统计工具箱中各种分布的随机数序列产生也可使用各自对应的命令。命令有如下一些： betarnd，binornd，chi2rnd，exprnd，frnd，gamrnd，geornd，hygernd，lognrnd，nbinrnd，ncfrnd，nctrnd，ncx2rnd，normrnd，poissrnd，raylrnd，trnd，unidrnd，unifrnd，weibrnd。读者如果用到时可查询相关帮助文档。6.3信号特征参数的计算 6.3.1离散信号的能量和功率计算方法 时间连续信号f(t)在时间t∈［t1,tn］内的能量定义为该信号在1Ω电阻上消耗的能量，即 (6-15) 在时间f∈［t1,tn］内的平均功率定义为(6-16) 设函数f(t)在t∈［t1,tn］内的离散点{t1,t2,…,ti,…,tn}上的值为f(ti)，且这些离散点间隔距离足够小，则以上求积分可以近似为 (6-17) 以及(6-18) 其中，Δt=ti-ti-1为相邻离散点的间隔。 事实上，从以上推导中已经得出了时间离散信号的功率计算方法，即对于离散时间信号f(k),k=1,2,…,n，其采样间隔为Δt，则其能量计算公式为 而平均功率的计算公式为(6-19)(6-20) 6.3.2离散信号的直流分量和交流分量 离散时间信号f(k)(k=1,2,…,n)的直流分量就是信号样值的平均值，即 而f(k)(k=1,2,…,n)的交流分量为(6-21)(6-22) 从而有交流分量的功率计算公式： 序列xi的方差定义为 其中，是xi的平均值，即(6-23)(6-24)(6-25) 对比交流功率计算公式，显然两者是一致的，于是我们得到一个重要概念，离散时间序列的交流功率就是该序列的方差。MATLAB中使用命令“mean”求序列的均值，用命令“var”来求序列的方差。 现在举例说明信号直流分量、交流分量以及功率的计算方法。设模拟信号为f(t)=1+sin2π100t+2cos2π150t,其中，时间t∈［0,5］s，采样率为1000次/s。显然，信号f(t)的直流分量为1，直流功率为PDC=12=1;交流分量有2个，频率分别是100Hz和150Hz，功率分别是PAC1=12/2=0.5和PAC2=22/2=2,交流总功率为PAC=0.5+2=2.5。对f(t)取样后得到序列f(n)，下面从序列f(n)来计算上述信号参数。程序如下：

程序6-21 Fs=1000;Ts=1/Fs;%采样率和采样间隔 t=0:Ts:5;%信号持续时间

f=1+sin(2*pi*100*t)+2*cos(2*pi*150*t);%离散信号序列 f-DC=sum(f)./length(f)%计算信号的平均值——直流分量 f-AC=f-f-DC;%计算信号的交流分量 Average-Power-of-f=sum(f.^2)./(length(f)-1)%信号总功率 Power-of-f-DC=f-DC.^2%信号直流分量的功率 Power-of-f-AC=sum(f-AC.^2)./(length(f-AC)-1) %信号交流分量的功率 %下面是计算结果,对比以上分析的理论结果 f-DC= 1.0004 Average-Power-of-f= 3.5018 Power-of-f-DC= 1.0008 Power-of-f-AC= 2.5008 %事实上,用MATLAB提供的mean,var函数可以更加方便地计算直流分量和交流功率(方差) f-DC=mean(f) Power-of-f-AC=var(f) %计算结果: f-DC= 1.0004 Power-of-f-AC= 2.5008 6.3.3离散信号的统计特征参数的计算 离散信号的统计特征参数除了均值、方差之外，常用的还有标准差、相关系数、协方差矩阵、互相关系数、互协方差函数、中值、最小值、最大值等等。命令“s=std(X)”用来求离散序列X的标准差。标准差即方差开平方根的结果。命令“C=cov(x，y)”计算序列x和序列y的协方差矩阵。命令“R=corrcoef(x，y)”计算序列x和序列y的互相关系数，而命令“c=xcorr(x，y)”用于估计互协方差函数。计算序列的中值、最小值、最大值分别采用命令“median”、“min”和“max”，这些命令的详细情况和所使用的计算方法请参考在线帮助。6.4信号的频谱分析 6.4.1频谱分析的原理 1.模拟信号的频谱分析——连续时间傅立叶变换(CTFT)对于时间连续信号x(t)，其频域分析可以采用连续时间傅立叶变换(CTFT)来进行： (6-26) 其中正变换为

反变换为 连续时间傅立叶变换特别适合于对时间连续信号的理论分析(如信号与系统课程中的内容)，但是，由于其变换两边的函数x(t)和X(ω)都是连续函数，不适合于计算机处理，因此在进行数值计算的时候，必须将其离散化，然后利用离散傅立叶变换(DFT)来近似实现。(6-27)(6-28) 2.离散信号的频谱分析——离散傅立叶变换(DFT)

离散傅立叶变换的定义是，对于N点的离散序列x(n)，其离散傅立叶变换(DFT)也是N点离散序列，设为X(k)，即 (6-29) 其中正变换为(6-30)k=0,1,…,N-1 反变换为 3.利用离散傅立叶变换(DFT)来近似计算连续时间傅立叶变换(CTFT) 设x(t)是一个连续时间信号，为了对其进行离散化处理，首先根据取样定理确定采样时间间隔T，然后再确定截取的时间段L，由这两个参数来确定离散时间序列x(n)的点数N。例如，根据实际经验，若信号x(t)的极大部分能量集中在频率段［0,fm］，且x(t)的非零值区间集中在时间段［t0,t1］，则选取采样时间间隔T≤，选取截取的时间段L≥|t2-t1|。于是得到离散序列的点数N=［L/T］+1。n=0,1,…,N-1(6-31) 用t=nT,n=0,1,…,N-1以及N=［L/T］+1对连续信号x(t)进行离散化，对x(t)的连续时间傅立叶变换中的积分做近似求和计算，有 但是，上式的结果X(ω)仍然是连续函数，计算机不能直接加以处理。为了实现数值计算，还需要对X(ω)作离散化处理，为此，在频率ω=kω0= ,k=0,1,…,N-1处对频谱X(ω)进行离散化，利用上式，有(6-32)

(6-33)

对比DFT计算公式，显然有 该式表明，利用DFT(FFT)计算连续时间傅立叶变换的频谱时，除了计算时域样点的离散傅立叶变换的频谱X(k)，还要将X(k)乘以取样时间间隔T，才能得到结果。(6-34) 4.快速傅立叶变换 根据离散傅立叶变换(DFT)的公式，有 这可以用MATLAB编程实现，但是直接利用DFT定义来计算是低效率的。20世纪60年代末期人们找到了一种DFT的高效算法，这就是我们常说的快速傅立叶变换算法(FFT)，关于FFT的算法原理请参考数字信号处理相关书籍。利用FFT算法可以大大提高DFT的计算效率。MATLAB提供了函数fft可以直接利用。,k=0,1，…,N-1(6-35)

另外，对于二维信号的离散傅立叶变换的快速变换可用指令“fft2”完成，更一般地，n维FFT指令是“fftn”。指令“fft”的常用语法是：

Y=fft(X)%若X为1行N列或N行1列序列,则计算N点DFT

Y=fft(X,n)%计算n点DFT %若X序列的点数小于n,则自动在X序列后补零 %若X序列的点数大于n,则自动对X序列截断离散傅立叶变换DFT(FFT)的反变换函数是“ifft”。其调用语法是:

y=ifft(X)%若X为1行N列或N行1列序列,则计算N点逆DFT

y=ifft(X,n)%计算n点逆DFT %若X序列的点数小于n,则自动在X序列后补零 %若X序列的点数大于n,则自动对X序列截断 5.确定信号的频谱计算 我们利用FFT来对确定信号的频谱进行近似计算，并与理论计算作对比。下面举例说明。已知时间连续信号x(t)=e-tu(t),-∞<t<+∞，求其傅立叶变换，并利用FFT进行数值计算。根据连续信号傅立叶变换的公式计算得信号的频谱为 于是就可以作出频谱曲线,如图6-24所示。 图6-24理论计算的时域波形x(t)=e-tu(t)及其频谱程序如下：

程序6-22

t=-1:0.01:10; x-t=exp(-t).*(t>0); subplot(3,1,1);plot(t,x-t);%时域波形 axis(［-110-0.11.1］);xlabel(′timesec′);(6-36) w=-40:0.1:40; X-w=1./(1+j*w); subplot(3,1,2);plot(w,abs(X-w));%频域幅度谱 axis(［-404001.1］);xlabel(′freqrad′); subplot(3,1,3);plot(w,angle(X-w));%频域相位谱 axis(［-4040-pi/2pi/2］);xlabel(′freqrad′); 我们再用FFT进行数值计算。从图6-24中可见，该信号的时域为无限长的，其频谱也是无限宽的，因此离散化的时候必须对其进行时域截断和频域截断。从图中分析，时间大于5s后波形近似为0，频率大于40rad/s的部分值较小，可以截断。更加严格地，应该选取大于95%或99%的总能量时域和频域部分。 于是我们得到对x(t)的截断参数：时域:t∈［0,5］，频域:ω∈［0,40］(单边)。相应的采样间隔T为(6-37) 采样点数为 (6-38) 程序如下： 程序6-23

w-m=40;%截断频率 T=pi/w-m;%采样间隔 L=5; t=0:T:L;%时域截断 x-t=exp(-t).*(t>0);%信号序列 N=length(x-t);%序列长度(点数) %------------------------ X-k=fft(x-t);%FFT %------------------------ w0=2*pi/(N*T);%离散频率间隔 kw=2*pi/(N*T).*［0:N-1］;%离散频率样点 X-kw=T.*X-k;%乘以T得到连续傅立叶变换频谱的样值 %------------------------- plot(kw,abs(X-kw),′.′,′MarkerSize′,14);%作出数值计算的幅度谱点 axis(［-4090-0.11.1］);xlabel(′freqrad/s′); holdon;%保持前面作图曲线不被擦除 w=-40:0.1:40; X-w=1./(1+j*w);%理论计算频谱表达式 plot(w,abs(X-w));%作图对比 结果如图6-25所示，曲线为理论计算结果，小点为DFT计算结果。可见利用DFT(FFT)很好地估计了连续时间信号的频谱。值得注意的是，频率在0~40rad/s和40~80rad/s的DFT计算的幅度谱以40rad/s(截断频率处)为轴对称。而理论频率谱的计算则是以0频率为对称的,这是由于DFT固有的时域频域离散化的结果。另外，利用函数“fftshift”可以将fft的计算输出零频搬移到输出中心，从而与理论结果的范围匹配。图6-25频谱的理论计算值与利用DFT(FFT)计算值之间的对比 6.信号的时域频域能量和功率关系——帕斯瓦尔定理 对于时间连续信号及其傅立叶变换x(t)X(ω)，有能量关系 称为连续时间傅立叶变换的帕斯瓦尔能量定理。 对于N点的离散序列及其DFT:x(n)X(k)，有能量关系(6-39)(6-40)

下面用实例加以验证，设有一连续信号x(t)=1+sin2π50t,求其功率。显然，该信号的理论计算功率为(6-41)(6-42) 由于该信号中最高频率成分为50Hz，故进行离散化时采样频率大于100Hz即可，这里我们取采样率为200Hz，信号时间截取为0~10s。程序如下： 程序6-24 fs=200;Ts=1/fs; t=0:Ts:10; x-t=1+sin(2*pi*50*t); N=length(x-t); P-av=sum(x-t.*x-t)./N%由时域计算信号功率 P-av= 1.4998 X-k=fft(x-t,N); P-av=sum((abs(X-k)).^2)./(N.^2)%由频域计算信号功率 P-av= 1.4998

6.4.2频谱分析的若干问题 通过前面分析我们知道，对于连续非周期信号的频谱的数值计算必须首先对信号进行时域取样，得到时间离散化的信号，时域取样必须满足或近似满足取样定理，根据时域频域的对应关系，时域取样将导致得到的抽样信号频谱周期化。然后，为了使得周期化后的抽样信号频谱便于计算机处理，还必须将这个频谱进行离散化，方法是对该频谱进行抽样。根据时域和频域的对应关系，频域离散化将对应于时域信号的周期化。因此，对连续非周期信号的频谱进行数值计算时，要确定如何截取信号的时间段，如何选择时域采样率，以及如何在某时间段上对信号进行截取。

截取信号的时间段的长度决定了时域周期化的周期，相应地对应于频域抽样的间隔，即频率分辨率。时域采样率决定了频域周期化的周期，即频谱数值计算的范围。而在某时间段上对信号进行截取的方法，也即不同窗函数的应用，决定了信号频谱估计的精度和实用范围。 设连续非周期信号为f(t)，要分析该信号在频率范围［0,fm］的频谱,而要求分析的频率分辨率为ΔfHz。 1.根据该信号的频率范围确定采样率 要分析该信号在频率范围［0,fm］的频谱，采样率fs必须满足采样定理：fs≥2fm。相应地，采样时间间隔T（也称为时间分辨率）满足：。 2.根据频率分辨率要求确定分析信号f(t)的截取时间段长度要使得分析的频率分辨率达到Δf，即每隔频率Δf计算一个频率点，那么对信号的截取时间长度L必须满足：L≥1/Δf