MATLAB系统仿真

1、MATLAB编程与系统仿真重庆交通大学信息科学与工程学院通信工程系主讲教师:谭晋Mobilemail:2009年10月Content Matlab编程基础p变量及其赋值变量及其赋值p运算符与表达式运算符与表达式p程序控制流程程序控制流程p自定义函数自定义函数p基本绘图方法基本绘图方法pM文件及其调试文件及其调试 Matlab在信号处理中的应用p信号及其表示、基本运算信号及其表示、基本运算p系统的响应系统的响应p傅里叶变换傅里叶变换p数字滤波器设计数字滤波器设计 通信系统仿真p编译码编译码p模拟调制与解调模拟调制与解调p数字调制与解调数字调制与解调p性能仿真性能仿真 S

2、imulink应用基础课程考核方式课程考核方式 考核方式p上机实验上机实验(实验报告实验报告) 30分分p平时成绩平时成绩(考勤考勤+作业作业)10分分p上机考试上机考试 60分分 上机考试两道题上机考试两道题: 25分分+35分分 上机考试时间:上机课中最后一次上机时间Matlab编程基础编程基础 变量及其赋值p标识符号与数据格式标识符号与数据格式标识符号是标志变量名、常量名、函数名和文件名的字符串的总称。双精度数据格式，64位二进制。p矩阵及其元素的赋值矩阵及其元素的赋值赋值语句：变量=表达式（数）Matlab中变量都代表矩阵，其阶数为MN赋值要求：u整个矩阵的值应在中；u同一行中各元素之

3、间以逗号“，”或空格分开；u不同行的元素以“；”隔开。例P15-16Matlab编程基础编程基础矩阵元素的索引u访问向量的元素只要使用相应元素的索引即可赋值技巧u利用冒号“：”给全行的元素赋值例：a=eye(5); a(2,:)=ones(1,5) au利用行、列标注构成新的矩阵例：b=a(2,3,4,1,2,3)零矩阵与空矩阵的区别？Matlab编程基础编程基础利用小矩阵组成大矩阵例：a=ones(2,3) b=zeros(3,2) c=a;a;a d=b;b e=c,d常用生成矩阵的方法u利用生成矩阵的函数（P17-18）u利用“：”运算符冒号运算符创建了具有10个元素的向量。利用冒号运算

4、符创建向量的基本语法如下：X = J:INC:KMatlab编程基础编程基础Matlab内部的特殊变量和常数例：ans、eps、realmax、realmin、pi、Inf 、i、j、NaN复数的赋值方式u可将矩阵元素逐个赋予复数；u将矩阵的实部和虚部分别赋值。Matlab编程基础编程基础 运算符与数学表达式p算术运算符（算术运算符（P20表表2.2）矩阵加减：矩阵加减：“+”、“-”功能：两矩阵对应元素相加/减。A、B矩阵的阶数必须相等或其中之一为标量，标量可与任意大小的矩阵相加/减。矩阵相乘：矩阵相乘：“*”矩阵元素相乘：矩阵元素相乘：“.*”功能：矩阵A、B对应元素相乘。A、B行列必相同

5、或其中之一为标量。Matlab编程基础编程基础矩阵右除矩阵右除B/A与左除与左除AC如果A为一非奇异矩阵，则右除：B/A=B*inv(A)如果A为一非奇异矩阵，则左除：AC=inv(A)*C矩阵元素右除矩阵元素右除A./B与左除与左除A.BA,B两矩阵必同行列，或其中之一为标量。矩阵幂矩阵幂“”及矩阵元素幂及矩阵元素幂“.”矩阵转置矩阵转置“”“”及矩阵非共轭转置及矩阵非共轭转置“.”Matlab编程基础编程基础p关系运算符关系运算符、=、=、=Matlab关系运算符用来比较两个同样大小的数组，或用来比较一个数组和一个标量；当数组和标量进行比较时，结果与数组大小一样；关系比较的结果只有两种可能

6、，即1或0。例：比较下面两段语句的结果x=(-3:3)/3 x=x+(x=0)*epssin(x)./x sin(x)./x例：统计数组A中值大于5.5元素的个数Matlab编程基础编程基础p逻辑运算符逻辑运算符&(与与) |(或或) (非非)真值表（P24）p其他逻辑函数（其他逻辑函数（P25表表2.6,2.7）p数学表达式的数学表达式的MATLAB描述（描述（P25）不要忽略数乘不能调用系统中不存在的函数Matlab编程基础编程基础 流程控制pif 语句语句格式格式:格式1If 表达式 语句组AEnd例例1_1:随机生成一个随机生成一个(0,1)之间的数之间的数,如果该数大于等如果

7、该数大于等于于0.5则输出该数则输出该数. Matlab编程基础编程基础格式2If 表达式 语句组AElse 语句组BEnd例1_2: 随机生成一个(0,1)之间的数,如果该数大于等于0.5,则将该数变为1,否则将该数变成-1.例1_3: 随机生成100个(0,1)之间的数,将大于等于0.5的变为1,小于0.5的变为0.Matlab编程基础编程基础格式3If 表达式1 语句组AElseif 表达式2 语句组BElse 语句组CEnd例1_4: 随机生成一个(0,1)之间的数,如果该数所处区间为(0,0.3)则将该数变成0,如果该数所处区间为(0.3,0.6),则将该数变成1,否则将该数变成2.

8、 Matlab编程基础编程基础pSwitch语句语句Switch语句称为条件选择语句语句称为条件选择语句,其中关键字包括其中关键字包括:switch, case, otherwise以及以及end.主要用于有选择性的程序设计主要用于有选择性的程序设计,实现程实现程序的多分支选择序的多分支选择.格式格式:Switch 选择表达式 case 情况表达式1 语句组1 case情况表达式2 语句组2 otherwise 语句组nend 例例1_5：随机生成一个：随机生成一个1到到12之间的之间的整数，按如下规则进行转换：整数，按如下规则进行转换：如果该数在如果该数在14之间，转换成之间，转换成-1；如

9、果该数在如果该数在58之间，转换成之间，转换成0；如果该数在如果该数在912之间，转换成之间，转换成1。Matlab编程基础编程基础pWhile语句语句While语句为条件循环语句，循环次数不确定，取语句为条件循环语句，循环次数不确定，取决于逻辑表达式。决于逻辑表达式。格式：格式：While 表达式 语句组AEndMatlab编程基础编程基础pFor语句语句是一种循环次数确定的循环语句。是一种循环次数确定的循环语句。格式：格式： for V=初值：增量：终值初值：增量：终值 语句组语句组A end功能：将语句组功能：将语句组A循环执行循环执行N次，其中：次，其中： N=1+（终值（终值-初值）

10、初值）/增量增量Matlab编程基础编程基础pFor循环注意事项循环注意事项不能通过在循环体内重新赋值给循环变量来终止FOR循环,只能使用Break语句;例:Matlab编程基础编程基础循环被执行前,应预先分配数组;当有一个等效的数组方法来解答给定的问题时,应避免使用循环;循环可按需要进行嵌套;例:(前面while循环的例子)Matlab编程基础编程基础 例: 随机生成01000之间的随机数100个,并按从小到大的顺序输出. 例: 随机生成51100之间的随机数50个,请输出序列中不存在的51100之间的数. 例: 请验证两个正交的复噪声信号之和的包络服从Rayleigh分布Matlab编程基

11、础编程基础 数据的输入/输出及文件的读写(P31)psave, load命令命令 基本数学函数(P42) 基本绘图方法p图形窗口控制命令图形窗口控制命令figure(N)功能:新创建一个图形窗口,并自动/N编号图形窗口内容保持命令hold on (off)功能:保持当前图形窗口的内容,或取消hold on功能,恢复系统默认状态.图形窗口划分命令格式: sumplot(m,n,p)功能: 将图形窗口分成n*m个子图形窗口,并选择第p个子图形窗口作为当前图形窗口,供绘图函数作图使用.Matlab编程基础编程基础p二维图形的绘制二维图形的绘制二维图形绘制命令plot格式一: plot(y)功能:输入

12、一个数据的情况. 如果y是一个数组,函数plot(y)给出笛卡儿坐标下的二维图.以y中元素的下标为X坐标,以y中元素值为Y坐标.例:给出一个周期的sin(x)的图形格式二: plot(x,y)功能:以x数组中元素值为横坐标,以y元素值为纵坐标作图.例1_10:给出一个从(-pi,pi)的cos(x)图形Matlab编程基础编程基础线型,标记符及颜色的设置(P50)格式一:plot(y,字符串)格式二:plot(x,y,字符串)其中:字符串中是关于线形,颜色及标记符号的说明图形的标注、网格及图例说明(1)添加图形标题命令title格式: title(string)功能:在当前坐标系的顶部加一个文

13、本串string,作为该图的标题.(2)添加坐标轴标志函数xlabel,ylabel,zlabel格式:xlabel(string) ylabel(string) zlabel(string)Matlab编程基础编程基础(3)设置网格线命令grid格式: grid on(off)功能: 对当前坐标图加上网格式或撤销网格线.(4)图形标注函数legend格式:legend(string1,string2,string3,.)功能:在当前图中添加图例.p坐标轴刻度设置坐标轴刻度设置格式: axis(xmin,xmax,ymin,ymax,)功能:对x,y,等坐标轴进行标定.Matlab编程基础编程

14、基础p多条曲线的绘制多条曲线的绘制格式:plot(x1,y1,x2,y2,x3,y3,)功能:为每一对X-Y绘制一条图线.可以在每个xi,yi后加线形,标记和颜色控制.例1_11:在一张图中画出sin(t),cos(t),sin(t)+cos(t)的图形Matlab编程基础编程基础 M文件命令语句构成的程序文件命令语句构成的程序文件,可分为两种可分为两种:p主程序文件主程序文件由用户为解决特定问题,而将命令语句序列放到m文件中而形成的文件主程序文件结构u有关程序的功能,使用方法等内容的注释部分以”%”开始直至行末的部分为注释部分u程序主体由若干条命令语句组成,实现程序设计功能,通常用clear

15、, close all开始u主程的运行方式例1_12编写程序画出sin(6pix),cos(2pix)以及两函数之各的曲线.Matlab编程基础编程基础p函数文件结构函数文件结构函数定义行function 变量列表=函数名(参数)H1行注释行,主要说明该函数的功能,可以在help命令中显示.函数帮助文本注释行,用来比较详细地说明函数的相关信息,可在help命令中显示函数体用来完成函数功能的语句序列,并以end语句为结束标志.例例1_13生成四个长度为64的序列,分别为s1,s2,s3,s4,对每个序列乘上1,-1,j,-j中的任意一个因子后求和.求所有可能的和序列中峰均功率比最低的序列及相应的

16、因子.系统仿真基础系统仿真基础 通信系统仿真发射机接收机信道干扰发射机发射机抽样量化编码信源编码信道编码调制信源线性分组码循环码BCHRS卷积码TurboASKFSKPSKMSKGMSKQAM 信源编码抽样量化编码信源编码信源function signal=informationSource(t)function sampleSignal=sampling()function quantizeData=quantizing(sampleSingal)function bitstream=sourceCoding(quantizeData)信源信源 function signal=informa

17、tionSource(t) signal=sin(2*pi*t)+cos(4*pi*t); end 00.10.20.30.40.50.60.70.80.91-2-1.5-1-0.500.511.5抽样抽样 function sampleSignal=sampling() t=0:0.01:1 sampleSignal=informationSource(t) end020406080100120-2-1.5-1-0.500.511.5量化量化 function quantizeData=quantizing(sampleSingal) dv=2*a/(M-1); quantizeData=s

18、ign(sampleSignal).* fix(abs(sampleSignal+dv/2)/dv)*dv; end020406080100120-2.5-2-1.5-1-0.500.511.522.5量 化 电 平 数 32020406080100120-2.5-2-1.5-1-0.500.511.522.5量 化 电 平 数 256编码编码function bitstream=sourceCoding(quantizeData,dv,M) N=length(quantizeData); bitstream=; N1=fix(log2(M); for i=1:N NData=fix(abs(

19、quantizeData(i)/dv); bindata=convert2bin(NData,N1); if (quantizeData(i)=0) bindata(1)=0; else bindata(1)=1; end bitstream=bitstream,bindata; endend找到量化电平编号将该编码转换成二进制加上符号位处理下一个量化电平编码编码 将一个整数转换成二进制数据function bindata=convert2bin(NData,N) bindata=zeros(1,N); i=N; while (NData=0) bindata(i)=rem(NData,2);

20、 NData=fix(NData/2); i=i-1; end if (NData(1)=1) NData=0,ones(1,N-1); endend信道编码信道编码信道编码线性分组码循环码BCHRS卷积码Turbo为方便,使用奇偶校验码将原比特流,每7位一组,添加校验位1位.采用偶校验.该码不具纠错功能,但具检一位错的功能.信道编码信道编码function datastream=channelcodingEV(bitstream) N=fix(length(bitstream)/7); datastream=; for i=1:N temp=bitstream(i-1)*7+1:i*7);

21、datastream=datastream,temp,rem(sum(temp),2); endend调制调制ASKFSKPSKMSKGMSKQAM本例采用2ASK调制方式,其它调制方式可以类似进行仿真.Bit速率为1000bps;载波频率为10000Hz;由每一位bit占用时间:1/1000=1ms;载波的一个周期所占用时间:1/10000=0.1ms因此,一个bit位有10个完整的载波周期.载波的一个周期以16点抽样.即 t=2*pi/16:2*pi/16:2*pi; carrier=cos(t);调制时:1)如果某位为1,则复制10次carrier;2)如果某位为0,则复制10*16个0

22、.调制function transmittedSignal=BAskModu(signal,bitRate,carrierFrequency,N1) M=carrierFrequency/bitRate; t=(2*pi/N1):(2*pi/N1):2*pi; carrier=; for j=1:M carrier=carrier,cos(t); end transmittedSignal=; for i=1:length(signal) if (signal(i)=1) transmittedSignal=transmittedSignal,carrier; else transmitted

23、Signal=transmittedSignal,zeros(1,length(carrier); end endend发射机发射机Transmitterfunction transmittedSignal=transmitter() sampleSignal=sampling() a=2; M=256; quantizeData,dv=quantizing(sampleSignal,a,M); bitstream=sourceCoding(quantizeData,dv,M); datastream=channelcodingEV(bitstream) bitRate=1000; carri

24、erFrequency=10000; N1=16; transmittedSignal=BAskModu (datastream,bitRate,carrierFrequency,N1);end信源编码信道编码调制信道仿真信道仿真 本例使用加性高斯白噪声p直接叠加平均功率为直接叠加平均功率为sig的噪声信号的噪声信号 如如: signal noise=sqrt(sig2)*1/sqrt(2)*randn(1,length(signal); signal=signal+noise;p利用awgn函数 signal=sqrt(2)*signal; signal=awgn(signal,20,0)信

25、道仿真信道仿真function signal=gussianChannel(transmittedSignal,sig) N=length(transmittedSignal); noise=sqrt(sig)*1/sqrt(2)*randn(1,N); signal=transmittedSignal+noise;end接收机设计接收机设计译码反量化解调信道译码信源译码再现线性分组码译码循环码译码BCH译码RS译码卷积码译码Turbo译码ASK解调FSK解调PSK解调MSK解调GMSK解调QAM解调 bitstream=demoASK(receivedSignal,bitRate,carri

26、erFrequency,N);datastream=deEVChannelCoding(bitstream);restoreOriginalSignal=D2A(datastream,dv);ASK解调解调在仿真时假设能精确同步!ASK解调解调BPF整流LPF抽样判决数据速率=1000bps;载波频率=10000Hz;抽样频率=载波频率*16;信号0点带宽(单边) =2/ =2 /(1/1000) =1000Hz信号的通带= 9000Hz,11000Hz根据信号速率可得:基带信号的带宽=0,1000Hz假设所有同步有能精确同步!滤波器设计滤波器设计 几个基本概念p通带、通带截止频率、通带容限通

27、带、通带截止频率、通带容限p阻带、阻带容限、阻带截止频率阻带、阻带容限、阻带截止频率p过渡带过渡带p通带最大衰减通带最大衰减p阻带最小衰减阻带最小衰减滤波器设计方法滤波器设计方法 FDATool IIR(MATLAB直接设计)p利用利用IIR滤波器阶次估计函数求最小阶数和截止滤波器阶次估计函数求最小阶数和截止频率（频率（P162，P163） FIR(MATLAB直接设计)pfir1pfir2ASK解调解调function bitstream=demoASK(receivedSignal,bitRate,carrierFrequency,N) load mypass den=1; signal1=filter(bandpassNum,den,receivedSignal); signal1=abs(si