lab01Matlab平台的使用记得带耳机2012年

lab01Matlab平台的使用记得带耳机2012年 Matlab平台的使用 1.总论 Matlab是信号处理中功能强大且广泛使用的工具。本实验的主要目标是使你们熟悉使用Matlab。我们假定你们有过其他某种语言的编程经验，应该已经了解变量，循环和函数等基本概念。这里介绍一下Matlab的基本编程结构。 本实验的三个目标: 1.学习基本的Matlab命令和语法，包括“帮助”系统 2.学习编写你自己的Matlab的脚本文件，并像命令一样运行它们 3.学习一点Matlab的高级编程技术，例如向量化 1.1要提交的内容 提交的实验报告要求回答所有下面被方框框住的问，还有方框部分要求提交的Matlab代码。实验报告要求包括: 1.姓名，学号，日期和实验报告的题名。 2.一两句介绍本实验的内容 3.所有方框部分的回答和代码，答案要求简洁清楚 4.几句结论，重点谈谈你遇到的问题以及怎么解决的，或者你用了什么较好的方法，或者你有什么收获和感受 实验报告文档名为“班号-学号-姓名”。只交电子版。 一周内，提交给各班学习委员，学习委员再打包提交给我，wukn@hust.edu.cn 2. 初识Matlab 1. 启动Matlab，运行Matlab的帮助。在命令窗口输入“doc”，或在菜单中选“help”。 2. 在命令窗口中，试试: helpwin helpwin plot helpwin colon %<--- a VERY IMPORTANT notation helpwin ops helpwin zeros helpwin ones lookfor filter %<--- keyword search 3.运行Matlab的demo:在命令窗口输入“demo” 4.把Matlab当计算器用，试试: pi*pi - 10 sin(pi/4) ansˆ2 %<--- "ans" holds the last result 5.在Matlab中给变量赋值，试试: xx = sin( pi/5 ); cos( pi/5 ) %<--- assigned to what? yy = sqrt( 1 - xx*xx ) 平方根 ans 注意观察语句末尾有分号和无分号有什么差别。 有分号是赋值语句 6.Matlab中的复数运算的基本操作，试试: zz = 3 + 4i, ww = -3 + 4j real(zz), imag(zz) abs([zz,ww]) %<-- Vector constructor 取模运算 conj(zz+ww) 复共轭 angle(zz) 幅角 exp( j*pi ) 指数运算 exp(j*[ pi/4, 0, -pi/4 ]) 计算三个结果 3. 再识Matlab 3.1Matlab数组索引 (a).冒号的含义 试试如下命令的结果: jkl = 0 : 6 冒号 :就是colon向量 jkl = 2 : 4 : 17 从2开始，2+4，2+4+4~~~ jkl = 99 : -1 : 88 ttt = 2 : (1/9) : 4 tpi = pi * [ 0:0.1:2 ]; %注意这里有分号 (b).从向量提取数字和,或插入数字，思考如下xx的定义: xx = [ zeros(1,3), linspace(0,1,5), ones(1,4) ] 为向量赋值 xx(4:6) 取xx中第4到第6个向量 size(xx) 向量的维度 length(xx) 向量的长度 xx(2:2:length(xx)) 从2开始间隔2取向量 xx(2:2:end) 从2开始间隔2取向量 解释以上代码最后4行的结果 (c).观察如下赋值的结果: yy = xx; yy(4:6) = pi*(1:3) 现在写一条语句,xx用(b)中定义的方法，把xx的偶数索引的元素(即xx(2),xx(4)等)的 ,值替换为常数,。使用向量替换，不要用循环。xx(2:2:end)=pi^pi 3.2Matlab脚本文件 (a)看看关于向量的实验。把向量看作一个数字集合，试试: xk = cos( pi*(0:11)/4 ) %<---comment: compute cosines 解释一下余弦函数的值是怎样存储在向量xk中的。xk(1)是什么意思,有定义xk(0)吗,xk(1)是xk的首向量，没有定义xk(0) 注:语句后没有分号，所以语句会把xk的计算结果依次显示在命令窗口，在百分号,后面的文字是注释，可以忽略。 (b)Matlab中可以写循环，但是这不是最有效率的解决办法。尽量避免循环，并使用冒号会提高效率。下面的代码是用一个循环来计算余弦函数的值的: yy = [ ]; %<--- initialize the yy vector to be empty for k=-5:5 yy(k+6) = cos( k*pi/3 ) end yy 解释为什么需要写yy(k+6)。如果用yy(k)会发生什么,k+6使yy的赋值从1开始，如果用yy(k)则yy不能成功赋值，因为向量必须从一个正整数开始 不使用循环，而用3.1的方法，利用冒号，重新实现这个功能。 yy = [ ]; %<--- initialize the yy vector to be empty k=-5:5 yy(k+6) = cos( k*pi/3 ) yy (c)Matlab中实数和复数的绘制都很容易。基本绘图命令绘制向量yy和向量xx的关系，试试: xx = [-3 -1 0 1 3]; yy = xx.*xx - 3*xx; plot( xx, yy ) zz = xx + yy*sqrt(-1) plot( zz ) %<---- complex values: plot imag vs. real 在命令窗口输入“helpwin arith”，了解当xx是一个向量时，操作xx.*xx是如何工作的。当xx是一个矩阵时，操作xx.*xx又是如何工作的。 当xx是一个向量时，xx.*xx，是每个向量分别相乘，当xx是一个矩阵时，是矩阵乘法 当你没把握某个命令时，可随时使用帮助系统。 (d)使用Matlab的编译器创建脚本文件，命名为mylab1.m，文件内容如下: clear all; close all; tt = -1 : 0.01 : 1; xx = cos(5*pi*tt); zz = 1.4*exp(j*pi/2)*exp(j*5*pi*tt); plot(tt, xx, 'b-', tt, real(zz), 'r--') %<--- plot a sinusoid grid on title('TEST PLOT of a SINUSOID') xlabel('TIME (sec)') 其中红色虚线('r--')为real(zz)。 zz = 1.4*exp(j*pi/2)*exp(j*5*pi*tt) = 1.4*(0 + j)*(cos(5*pi*tt)+jsin(5*pi*tt)) = - 1.4*sin(5*pi*tt) + j*cos(5*pi*tt) real(zz)=-1.4*sin(5*pi*tt) =1.4cos(5*pi*tt+pi/2) 则幅度 A=1.4，初相位 φ=-π/2。 解释为什么real(zz)的图形是一个正弦曲线。它的相位和幅度是多少,由所绘图形的时移量计算相位。 (e)在命令窗口输入“mylab1”即可运行上述m文件，即 mylab1 %<---will run the commands in the file type mylab1 %<---will type out the contents of % mylab1.m to the screen 3.3 Matlab与声音相关的函数命令 这一节练习是关于声音信号的，所以需要带耳机来机房听效果。 (a)运行Matlab的声音demo:在命令窗口输入“xpsound”。你可以听到多种声音，并可以查看对应的时间序列，功率谱密度和声谱图(声谱图详见7.3)。 1(b)现在在Matlab中生成一个音调(例如正弦)，用soundsc()命令听一听 在3.2(d)的第2行和第3行创建了一个向量xx为2.5Hz的正弦序列。而你们要生成的正弦信号的频率为2000Hz，时长为0.9秒，采样率fs为11025个样本,秒。采样率说明了采样点之间的时间间隔，所以时间向量应定义为: tt = 0:(1/fs):dur; 其中fs为采样率，dur为时长(以秒为单位)。可以在帮助系统中查找sound()和soundsc()了解它们的用法和含义。那么，你们的tt向量的长度(即样本数)是多少, 4.用Matlab处理正弦信号 在实验报告中要包括本节的图形和实验结果描述。写一个Matlab的m文件来完成步骤(a)到步骤(d)。你们的实验报告也要包括这个m文件的代码。 (a)生成一个时间向量tt，需要覆盖下面的(b)中频率为4000Hz的正弦信号的两个周期。 d)的方式。如果我们使用T表示正弦的周期，定义向量tt的开始时tt的定义方式使用3.2( 间为,T，结束时间为，T，那两个周期就会包括tt=0。最后，确保你们正弦波的每个周期有至少25个样本。也就是说，当你们使用冒号操作符定义时间变量时，要使增量足够小，才能保证每个周期可以产生至少25个样本。 (b)生成以下两个4000Hz的正弦序列: x1(t) = A1 cos(2π(4000)(t ? tm1)) x2(t) = A2 cos(2π(4000)(t ? tm2)) 其中A1是你的年龄，A2=1.2A1，tm1 = (37.2/M)T ，tm2 = ?(41.3/D)T ，这里M和D分别是你生日的月和日，T是周期 绘制以上两个信号在-T?t?T的图形。为了输出下面(d)的要求，使用subplot(3,1,1)和subplot(3,1,2)来使图在一个图形窗口内，详见helpwin subplot 一般用法是: subplot(2,1,2); plot(xx); (c)创建第3个正弦信号为 x3(t) = x1(t) + x2(t). 在Matlab中，这表示把两个正弦向量中的值对应相加。绘制 x3(t) 在-T?t?T的图形，使用subplot(3,1,3)绘图。 (d)对每个图都添加一个图名，图名都要包含你的姓名，用title函数 详见 helpwin title, helpwin subplot 1 soundsc(xx,fs)函数需用两个参数，,xx是播放文件的数据的向量，fs是播放样本的采样率。另外，soundsc(xx,fs)会自动把xx的值尺度变化成,1到1之间，然后调用sound(xx,fs)播放信号。 clear all; close all; tt=-1/4000: 0.000000001 : 1/4000; x1=20*cos(2*pi*(4000)*(tt-37.2*4000/6)) ; x2=24*cos(2*pi*(4000)*(tt+41.3*4000/25)) ; x3=x1+x2; subplot(3,1,1); plot(tt,x1); title('luye x1'); subplot(3,1,2); plot(tt,x2); title('luye x2'); subplot(3,1,3); plot(tt,x3); title('luye x3'); grid on; xlabel('TIME (sec)') ; 4.1 复数幅度 用一条语句生成以上正弦信号x1(t) 的值，可以使用复数幅度表示方法: j,tx1(t) = Re{Xe } 其中X和w为某个常数。 5.复指数简介 本节目标是让大家熟悉复数，并了解怎么使用它们表示正弦信号(例如 x(t) = Acos(ωt + φ) )为 j,j,t复指数形式(例如z(t) = Aee)，方法如下所示: j,j,tx(t) = Acos(ωt + φ) = Re{ Aee } 5.2 Matlab中的复数 Matlab可用于计算复值式子，并可以向量显示结果，以下是一些Matlab的复数操作符: conj—复共轭 abs—幅度 angle—以弧度为单位的角度或相位 real –实部 imag—虚部 ,1i,j –已预定义为 x = 3 + 4i jθexp(j*theta) – 复指数函数 e 每个函数都以向量或矩阵作为输入参量。 注意:Matlab中没有 mag()和phase() 函数。 5.3 向量化 Matlab最重要的部分就是矩阵,向量语法。大多数情况下，循环都可以用向量操作代替，因为形如exp() 和cos() 等函数都被定义为向量输入，例如 cos(vv) = [cos(vv(1)), cos(vv(2)), cos(vv(3)), ... cos(vv(N))] 其中vv是N个元素的行向量。向量化可用于简化代码。如果有如下代码来绘制某个信号: M = 200; for k=1:M x(k) = k; y(k) = cos( 0.001*pi*x(k)*x(k) ); end plot( x, y, 'ro-' ) 你可以把循环替换掉，而只用如下3行代码即可: M = 200; y = cos( 0.001*pi*(1:M).*(1:M) ); plot( 1:M, y, 'ro-' ) 用这种向量化的方法编写2,3行代码完成以下Matlab代码，不使用循环。(注:当xx是向量时，xx*xx和xx.*xx是不同的) %--- make a plot of a weird signal N = 200; for k=1:N xk(k) = k/50; rk(k) = sqrt( xk(k)*xk(k) + 2.25 ); sig(k) = exp(j*2*pi*rk(k)); end plot( xk, real(sig), 'mo-' ) N = 200; rk = sqrt( ((1:N)/50).*((1:N)/50) + 2.25 ); sig = exp(j*2*pi*rk); plot( ((1:N)/50), real(sig), 'mo-' ) 5.4函数 函数是一种特别的M-文件，可以接受输入(矩阵或向量)，也可以得到输出。关键词function必须是文件的第一个词，并且定义了函数，M-文件的第一行定义函数如何传递输入和输出参量。文件扩展名必须使用小写的“m”，例如my_func.m。 如下函数有一些小错误，请找找错在哪里:(至少有3处错误) matlab mfile [xx,tt] = badcos(ff,dur) %BADCOS Function to generate a cosine wave % usage: % xx = badcos(ff,dur) % ff = desired frequency in Hz % dur = duration of the waveform in seconds % tt = 0:1/(100*ff):dur; %-- gives 100 samples per period badcos = cos(2*pi*freeq*tt); 错误有:第一个词必须是“function”。而且“freeq”在使用前没有定义。最后，第一行指明要“xx”作输出，因此“xx”需要在函数体中出现在至少一个赋值行的左边。 正确的函数写法应该是: function [xx,tt] = goodcos(ff,dur) tt = 0:1/(100*ff):dur; %-- gives 100 samples per period xx = cos(2*pi*ff*tt); 可在命令窗口输入函数的文件名(注意不是函数名)来调用它(也可以在另一个m文件调用它)，例如，如果以上goodcos函数的文件名为testcos.m，且ff=1，dur=10，则调用时应输入testcos(1,10)。你可以尝试一下，如果你输入goodcos(1,10),matlab会有什么反应。Undefined command/function 'goodcos'. 6.复指数 6.1 生成正弦信号的M-文件 写一个可以生成单一正弦信号x(t) = Acos(ωt+φ)的函数，使用4个输入参量:幅度A，频率,，相位,和时长dur。函数应当返回两个输出参量:正弦信号的值x和对应的时间t。确保函数生成的正弦信号在每个周期有20个值，函数名为one_cos()。提示:可借鉴上面的goodcos()函数。 绘制你们的one_cos()函数，参数选为:A=95，,=200,弧度,秒，,=,/5弧度，时长为0.025秒。推导所绘图形的周期和相位是否正确。如果周期以毫秒为单位是多少, function [x,t] = one_cos(A,w,B,dur) %幅度A，频率w，初相位B和时长dur t = 0:1/(20*(w/(2*pi))):dur; %-- gives 20 samples per period x = A*cos(w*t+B); plot(t,x) title('one_cos()函数') xlabel('时间t(s)') 若周期为毫秒，修改代码: function [x,t] = one_cos_ms( A, w, B, dur ) %幅度A，频率w，初相位B，时长dur t = 0 : 1/(20*(w/(2*pi))) : (dur); %-- gives 20 samples per period x = A*cos(w * t+B); plot(t*1000,x) title('one_cos_ms()函数' ) xlabel('时间(ms)') 7线性调频脉冲chirp 线性调频脉冲chirp信号(或线性扫频信号)是一个正弦信号，其频率从一个初值线性变化到终值。Chirp信号的表达式可写为: 2x(t)=cos(ψ(t))， 其中 ψ(t) = 2π μ t+ 2π f t + φ 0 对ψ(t) 求导后得到随时间线性变化的瞬时频率: f(t) = 2μt + f i 0 f (t) 的斜率等于 2μ，截距为f0。如果信号于t=0开始，f0也是初始频率。这种时变角度造i 成的频率变化称为频率调制。这类信号是频率调制FM信号的一个例子。更一般的，我们通常考虑它们为更大一类的角度调制信号。最后，因为频率的线性变化可以生成一个类似汽笛声或鸟叫声(chirp)，所以线性FM信号也可以称为chirp信号。 7.1chirp信号的Matlab合成方法 以下的Matlab代码可以合成一个chirp信号 fsamp = 11025; dt = 1/fsamp; dur = 1.8; tt = 0 : dt : dur; psi = 2*pi*(100 + 3000*tt + *tt.*tt); xx = real( 7.7*exp(j*psi) ); soundsc( xx, fsamp ); (a)确定合成信号的总时长(秒)，确定tt向量的长度(样本数) (b)在Matlab中，只能合成离散时间信号，所以对于chirp信号，我们处理为: 2x(t)=Acos( 2π μ t+ 2π f t+ φ) nn0n 其中，t=nTs表示离散时间值。在以上的Matlab代码中，tn的值是多少,A, ,, f和,的值是n0多少, (c)确定以上Matlab代码合成的频率的范围(Hz为单位)。手绘瞬时频率随时间的变化情况。听到的最小频率和最大频率是多少, (d)听听信号的频率是上升还是下降(使用soundsc())。注意soundsc()需要知道信号创建时的采样率。详见帮助。 7.2chirp的函数 使用以下代码来写一个可以合成chirp信号的函数，根据注释的提示编写。 function [xx,tt] = mychirp( f1, f2, dur, fsamp ) %MYCHIRP generate a linear-FM chirp signal % % usage: xx = mychirp( f1, f2, dur, fsamp ) % % f1 = starting frequency % f2 = ending frequency % dur = total time duration % fsamp = sampling frequency (OPTIONAL: default is 11025) % % xx = (vector of) samples of the chirp signal % tt = vector of time instants for t=0 to t=dur % if( nargin < 4 ) %-- Allow optional input argument fsamp = 11025; end 请生成一个chirp，其频率起始于2500Hz，终止于500Hz。时长应当为1.5秒。听听你的chirp。 7.3 声谱图 考虑信号的频谱非常有用。信号频谱是信号所含频率的表示方式。固定单一频率的正弦信号的频谱包含两个成分，一个在2,f0，一个在-2,f0。对于更复杂的信号，频谱会非常有趣，FM信号的频谱可以认为是时变的。一种表示时变频谱信号的方式就是声谱图。 声谱图可通过估计信号短时的频率成分得到。把频谱在每个部分的幅度画成关于频率和时间的二维的密度或颜色，即可得到一幅声谱图。 关于声谱图需要知道: 1.Matlab的函数specgram可计算声谱图。详见帮助，了解语法和所需参量。 2.声谱图是数值计算出来的，仅提供一个信号的时变频率成分的估计。 23.一般调用形式为specgram(xx,1024,fs)。第二个参量是窗长度，可以选用不同的长度观察 3不同尺度的声谱图。如果窗更长，例如1024或2048，声谱图可以查看单独的谱线。 为了看看声谱图是什么样，可以在命令窗口运行如下代码: fs=8000; xx = cos(3000*pi*(0:1/fs:0.5)); specgram(xx,1024,fs); 注意声谱图在正弦的频率处有一条水平的线。 7.4chirp的声谱图 使用mychirp函数合成chirp信号。实验参数为: 1.总时长为3秒，D/A转换率为fs=11025Hz 2.瞬时频率起始于5000Hz，终止于300Hz 听一听信号。Chirp信号的声音频率是否线性变化，频率变高还是变低, 7.5 有趣的chirp 再合成一个chirp信号，使用如下参数: 1.总时长为3秒，采样率为fs=11025Hz 2.频率起始于3000Hz，终止于-2000Hz(负频率) 听一听信号。频率是怎么变化的, 显示这个chirp信号的声谱图。 使用频谱理论(正频率成分和负频率成分)解释你听到的声音和看到的声谱图。 fsamp=11025; 2 如果第2项为“空矩阵”，则默认为使用 256 3 通常窗长度选为2的幂次，因为用了FFT计算。 dt=1/fsamp; dur=3; tt=0:dt:dur; psi=2*pi*(3000 + (4000/3)*tt-1000*tt.*tt); xx=real(exp( j* psi)); soundsc(xx,fsamp); specgram(xx,1024,fsamp); 负频率是由于傅立叶变换中出现 j 导致的。