Matlab-信号与系统实验.doc

实验一：利用Matlab来实现函数的绘制，主要思想就是我们给定一个定义域，然后再通过定义一个对应法则，使之生成一个对应的值域，两个集合存储于向量中，然后再一一对应，如果将两个向量看作两个集合，那么就可以将其视作两个集合满射。Exp：There are -5， 2， 3， 13， 10 我们定义对应法则f(x) = 2 * x 1 We get -11，3， 5，25， 19 这就是值域了然后画图过程就是建立平面坐标系，x轴用来标识定义域的值，y轴用来标识值域的值，绘制出来的图形就可以算是对f(x) 的仿真了。该实验要用到的函数：1. sin( x )这里对x变量要理解，x可以支持符号变量，即我们平时说的x，真正的未知数，然而由于这一实验我们运用的是如何通过Matlab中的向量（集合）来仿真信号，所以会用到一有限的定义域（向量），我们要做的就是将这个有限的集合当作变量传递给sin( x )，其内部是这样计算的。K = a, b, c, d, , z / K 是我们所给定的有限的定义域。sin( K ) 的返回值为 T = sin(a), sin(b), sin(c), sin(d), , sin(z) 很显然，K与T的长度也会是相同的。 在整个Matlab 编程时一定要有向量（向量也算是一维的矩阵，Matlab可是叫矩阵实验室呀）的思想。在仿真的过程中一定要考虑到这个函数是不是支持向量以点列独立的身份代入运算，如果是的话，那么是最好进行仿真的，给定一个有限的定义域（例如 K = -10 : 0.01 : 10， 产生向量 K，其值从-10 到 10，每个相邻的数差值为0.01，K的长度为 2001 ）括号内的K的长度达到了2001，在很小的区间内打印2001个点，就一定能够模拟出原函数的连续图形了。我们来看下面这个函数 function y = jump( x ) %阶跃函数 if( x 0 ) y = 1; else y = 0.5; endend咋看起来是没有什么问题的，学过c语言的我们知道，这个阶跃函数应该是对的。当问题就是，这个函数不支持输入的向量或者矩阵以点列的形式来处理，为什么，我们来一条条语句进行分析。首先是if( x 0 的结果是这样来判定的，首先 x 0 会返回一个boolean值的与x等长的向量，这个向量中的每一个元素对应与x中每一个元素进行一次 “ 0” 的判断，如果返回向量中有一个大于零的元素，那么x 0为假，所以K 代入后就是既不大于零又不小于零，所以结果就是返回的就是一个 1 * 1 的矩阵，y = 0.5；所以结果就只有一个点，这和我们的期望相去甚远啊。那么怎么改呢，有一种直接解决麻烦的改法： function y = jump( x ) %阶跃函数 len = length( x ); for( i = 1 : len ) if( x(i) 0 ) %直接对每一个元素进行判断 y(i) = 1; %赋值同样是对应赋值，不要动不动就y= 1了，这时赋值矩阵了 elseif( x(i) 0; %大于操作符是直接对每一个元素进行判定，返回等长的向量end这样写的话，就是利用了 操作符的性质，这在前面已经进行了说明。这个写法注意函数的返回值是boolean类型，也即逻辑值，也就是说，如果要使得x = 0时，y = 0.5的话，还要将先其转化为double类型。 function y = jump( x ) %阶跃函数 y = x 0; y = double( y ); fix = ( x = 0 ) .* 0.5; %乘号*默认是进行矩阵相乘的，直接用的话肯定不行（不满足条件） y = y + fix; %所以这里改成了数乘符，该运算符能够对向量中的每个元素进行操作end综上，最后这样写是最好不过的了 function y = jump( x ) %阶跃函数 y = ( x 0 ) + ( x = 0 ) .* 0.5; % y = ( x = 0 ) - ( x= 0 ) .* 0.5; end一时兴起，找来了源码，有谁能看懂吗？function Y = heaviside(X)%HEAVISIDE Step function.% HEAVISIDE(X) is 0 for X 0, and .5 for X = 0.% HEAVISIDE(X) is not a function in the strict sense.% See also DIRAC. % Copyright 1993-2008 The MathWorks, Inc.% $Revision: 1.1.6.1 $ $Date: 2009/03/09 20:41:30 $ Y = zeros(size(X);Y(X 0) = 1;eng = symengine;if strcmp(eng.kind,maple) Y(X = 0) = nan;else Y(X = 0) = .5;end.一个sin(x)就写了这么长，这叫我情何以堪.2.plot( )这个函数的用法比较多，这里用到了plot(x)，即传送一个变量，系统将自动生成横坐标（1 - length( x )），画出二维图形。plot( x, y )，以x为横坐标，y为纵坐标，画出二维坐标，此时x，y的长度必须相等。3.stem( )用来输出形状是火柴棍的图形，其他与plot差不多。t = -5 : 5;x = t;subplot( 1, 2, 1 );stem( t );subplot( 1, 2, 2 );stem( x, t );4.axis( a, b, c, d )该函数中的ax是axes（轴）的意思，后面是个谓词，说明这个函数是个表属性的函数。用于输出图形的区域，参数为一个向量，a, b 分别为横坐标左/右界，c, d分别为纵坐标的下/上界。切忌不要写成 axis( a, b, c, d )，这相当于给axis函数传递了四个矩阵，对吧。5.title( . )不多做解释了，这个函数就是在图形的上端输出你所想要打印的字符，后面还可以接一些参数，用来补充说明字符的颜色，大小等。6.xlabel( . ) / ylabel( . ) 用来在坐标轴旁边标注上注释，参数类型与title类似。7.grid用来画网格，常用的是 grid on / grid off操作。实验二：这次实验才算是初步接触到了信号，上一次的实验可以说是对函数的熟悉，并没有考虑到信号的这个层次。实验分两部分，第一部分即对已存在的信号 f(t) 通过编程得到信号 f( at + b ) 以及周期延拓。在没有接触到符号变量之前，我们还是用设出一个足够精度的定义域来仿真原函数的图形及其平移及尺度变换后的图形。t = -2 : 0.01 : 2; subplot( 2, 2, 1 );plot( t, exp( t ) );title( f(t) = et );axis( -2 2 0 60 );grid on; subplot( 2, 2, 2 );plot( t, exp( t + 2 ) );title( f(t) = e(t+2) );axis( -2 2 0 60 );grid on; subplot( 2, 2, 3 );plot( t, exp( -t ) );title( f(t) = e-t );axis( -2 2 0 60 );grid on; subplot( 2, 2, 4 );plot( t, exp( -t + 2 ) );title( f(t) = e(-t+2) );axis( -2 2 0 60 );grid on;周期延拓就是将已知区间的图形复制到新的定义域内去。比如说我要以 et （0t2）为周期，设定一个函数。解决这个问题的话，我们可以先设出一个周期的定义域，并生成这个周期内的对应定义域的值域集合。仿真这个函数就变成了生成更广的定义域区间的值，这部分的值是在已经得到的一个周期内的值域的重复。那么接下来有两条思路，一条是直接扩大区间进行赋值操作，另外一种则是定义一个函数，利用该取余操作直接生成一个定义域区间够大的一组集合。法一：t = 0 : 0.1 : 2; y = exp( t ); for( i = 0 : 3 ) for( j = 1 : length( y ) ) t( i * 21 + j ) = i * 2 + t( j ); z( i * 21 + j ) = y( j ); end end plot( t, z );法二：function Y = f( X ) Y = exp( mod( X, 2 ) );end注意函数应该保留在与函数同名的 “.m”文件中，这里应保留文件名为f.m t = 0 : 0.1 : 8; y = f( t ); plot( t, y );指导书上的方法也算是第二种方法。本次实验第二部分即卷积的实现，关于卷积，是一个函数在另一个函数上对一些点的连续积累效应。指导书中所说的以多项式的形式来模拟卷积的过程是很贴切的，但是我们最后求的是面积，所以要乘以一个增量。现将 f(t) = et * ( jump(t) - jump(t-4) )与 g(t) = jump(t) jump(t-4) 进行卷积。步骤如下：t = 0 : 0.1 : 4; ft = exp( t ) .* ( jump(t) - jump(t-4) ); gt = ( jump(t) - jump(t-4) ); y = 0.1 * conv( ft, gt );%求面积所以要乘以底长 plot( 0 : 0.1 : 8, y );PS：谈下我对卷积的理解吧， ，首先我们应该认为这是一个关于时间 t 的函数，而 只是我们用来表示积分的积分变量，关于 我们知道的是完全的，类似于我们求傅里叶变换时，就是知道了时域的全部信息后去得到频域的频谱函数。现在我们只考虑f( )在无限区间内绝对可积，因此积分上下限可以用一定的区间来近似代替，此时我们来看g( t - ) 这个函数，同样要求无限区间内绝对可积连同前面的f( )以及积分变量来看，我们更愿意将其看作关于 的函数，事实上一旦给定一个 t，那么积分式中 t 的作用就等同于常数。于是我要这样变形一下，看作是关于 的函数后，那么这个函数就是先将函数g( )关于y轴反转，再将函数向右平移 t 个单位。这样当我们从取时间为 -inf 到 inf g( t- )也就是将反转后的图像，从最左边平移到最右边的一个过程。在这个过程中的每一个时刻，我们叫一声暂停，于是我们将看见两个函数在平面的位置，我们严格遵守所给定的表达式，于是乎我们将它俩乘积不为零的区域进行积分。考虑到卷积是个关于时间 t 的函数，在我们叫一声暂停的时候，其实这就是当 t 取一定值时的卷积结果。右图为两函数f(t) = 0.2*t.*(heaviside(t)-heaviside(t-6);g(t) = exp(-0.2 * t)*(heaviside(t)-heaviside(t-5);的卷积示意图( t = 2 )。定义 z(t) = conv( f(t), g(t) )，则该图表示t=2时，两函数的位置图。该图与函数的关系为：该图是代入t的值的函数的几何解释。现在我们来分析Matlab中的conv()函数是如何来计算卷积的。苦于找不到好的模型对应解释离散型的数据卷积（计算机内最终都是以离散数据来处理的），用该处引入的多项式分析的话就非常和谐了。我们假设在两函数上( f( t ) 和 g( t )，暂不考虑g( t ) 的变化 )分别等距离dx取上100个点，计算出对应的值，保存并且按 x 轴方向一次组成一个一维向量。f( t ) 为 Fn * A(100 n + 1)g( t ) 为 Gn * A( n-1 )我们将取到的 Yn 看作是任意实数A的(n-1)次方。让我们想想最后得到 z(t) = conv( f(t), g(t) )中的 t 会有多少个取值。按照积分的思想，那么就有200个取值。自己可以画图得到这个答案，也可以这样理解（g( t )在运动过程中，每次有不同的乘积值时，我们就暂停这个过程，那么这个区间就是两者第一次相遇到最后一次相遇），我们开看当我们暂停时，乘积的值在这个模型下会有哪些特征。g( t )中在运动中首先与f( t )发生有效乘积的是g( t ) 的低次项（经过了反转），而f( t )相迎的同样是低次项。当第一次发生相遇时，相乘后的多项式的指数是 0，现在来分析第二次相遇，g( t )的0次项（用eg(0)表示）与f( t )的ef( 1 )相乘，eg(1 ) 与 ef( 0 )相乘，相加之后的指数是1。总结出来的公式就是，当第n次相遇： 依次类推，相同指数的项就代表了它们发生在同一时刻，啊哈，神奇的方法啊。得到Zn后，我们的工作并没有完成，卷积不是计算相乘之后的新的函数的面积吗？按理说应该是在两个系数相乘之后就马上乘以最小间距，但是我们把这一过程放到最后来做是一样的，即 conv( f( t ), g( t ) ) = Zn(t) * dx; 这里设n是关于t的函数。至此，卷积过程结束。实验三：该次实验终于来到了我们期待已久的频域，信号与系统的神奇之初就是将我们从时间域中解放出来了，来到了频域，我们会发现很多东西和之前的直觉感受不同，我们不再关心分析一个信号的各种特性，哪里得到最值，哪里斜率又开始增大.在频域里面我们只关心一个信号由哪些频率的基准函数组成。这就像我们不必了解一个多边形的具体形状，知道他们的边的方程或者密集的点集也就可以了。正交函数让我们也找到了这样的一个坐标系，三角函数，复指数函数。比如我们已知一个信号的周期为T，则这个信号的基频率就是 2*pi/T，假设我只取在这个频率左右M （Hz）的信号，则我们的任务就是计算这些合法分量的坐标值就行了，而这其中，对于周期信号和非周期信号的处理方式又稍有不同，周期信号的频谱为离散型（本身正交函数就是周期函数，所以只需要离散的频率的分解即可，可以理解为在基波上进行各种修饰，而且进行修饰的频率分量的周期又是要必须整除掉基周期的，换句话说，基周期必须是最大的周期，基频率必须是最小的频率）；而对与非周期函数，其频谱为连续型（基频率是无穷大的，其他分析类似），而且此时也只有分析相对大小才有意义，即频谱密度函数。理解并记住书本上的关于傅里叶级数的公式，运用Matlab进行仿真。 本次实验最值得称道的就是在不定义符号变量的情况下，利用矩阵的特性来求和。对于连续的信号，采用的策略是使得dt 尽可能的小，这样就能模拟出连续的情况来了。上面两个公式的重要性就不言而喻了。进行一个例题分析，求 f( t ) = e( -2|t| ) 的傅里叶变换。dt = 0.001;dw = 0.1;t = -100 : dt :100;y = exp( -2 * abs( t ) );N = 60;for k = -N : N F( k + N + 1 ) = y * exp( -j * k * dw * t ) * dt;endplot( -N*dw:dw:N*dw, real( F ) );1.由于没有指定时间范围，所以该函数是一个时间无穷的非周期函数，时间应取得尽可能的长。2.按理说dw应该为无穷小量，这里用0.1代替，太小运算速度很慢。3.我们认为取定频率范围为w = w | -6 = w = 6 121个点来作出图像。4.由于F向量只能从1开始编号，所以要给k一定的偏移量。5.最后输出图像，由于函数是偶函数，所以傅里叶系数中的虚部均为零。上面写的是非周期的傅里叶变换，周期函数的傅里叶级数求法与上面差不多，不同的是其dw用w0代替，而w0 相对于dw 来说一般会比较大，计算时还要除以周期T，当然输出时用stem()函数，这样能够凸显出其拥有离散的频谱图。在上面的基础上我们来进行傅里叶逆变换，知道了上面的编码流程后，这个逆变换是十分容易实现的，我们也可以根据最后出来的图形来验证我们的变换是否成功。还是贴出完整的代码吧，以便于大家能够运行这个程序，呵呵。clear all dt = 0