MATLAB系統(tǒng)仿真

1、MATLAB編程與系統(tǒng)仿真重慶交通大學(xué)信息科學(xué)與工程學(xué)院通信工程系主講教師:譚晉Mobilemail:2009年10月ContentMatlab編程基礎(chǔ)p變量及其賦值變量及其賦值p運(yùn)算符與表達(dá)式運(yùn)算符與表達(dá)式p程序控制流程程序控制流程p自定義函數(shù)自定義函數(shù)p基本繪圖方法基本繪圖方法pM文件及其調(diào)試文件及其調(diào)試Matlab在信號(hào)處理中的應(yīng)用p信號(hào)及其表示、基本運(yùn)算信號(hào)及其表示、基本運(yùn)算p系統(tǒng)的響應(yīng)系統(tǒng)的響應(yīng)p傅里葉變換傅里葉變換p數(shù)字濾波器設(shè)計(jì)數(shù)字濾波器設(shè)計(jì)通信系統(tǒng)仿真p編譯碼編譯碼p模擬調(diào)制與解調(diào)模擬調(diào)制與解調(diào)p數(shù)字調(diào)制與解調(diào)數(shù)字調(diào)制與解調(diào)p性能仿真性能仿真Simul

2、ink應(yīng)用基礎(chǔ)課程考核方式課程考核方式考核方式p上機(jī)實(shí)驗(yàn)上機(jī)實(shí)驗(yàn)(實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)報(bào)告) 30分分p平時(shí)成績(jī)平時(shí)成績(jī)(考勤考勤+作業(yè)作業(yè))10分分p上機(jī)考試上機(jī)考試 60分分 上機(jī)考試兩道題上機(jī)考試兩道題: 25分分+35分分上機(jī)考試時(shí)間:上機(jī)課中最后一次上機(jī)時(shí)間Matlab編程基礎(chǔ)編程基礎(chǔ)變量及其賦值p標(biāo)識(shí)符號(hào)與數(shù)據(jù)格式標(biāo)識(shí)符號(hào)與數(shù)據(jù)格式標(biāo)識(shí)符號(hào)是標(biāo)志變量名、常量名、函數(shù)名和文件名的字符串的總稱。雙精度數(shù)據(jù)格式，64位二進(jìn)制。p矩陣及其元素的賦值矩陣及其元素的賦值賦值語句：變量=表達(dá)式（數(shù)）Matlab中變量都代表矩陣，其階數(shù)為MN賦值要求：u整個(gè)矩陣的值應(yīng)在中；u同一行中各元素之間以逗號(hào)“，”

3、或空格分開；u不同行的元素以“；”隔開。例P15-16Matlab編程基礎(chǔ)編程基礎(chǔ)矩陣元素的索引u訪問向量的元素只要使用相應(yīng)元素的索引即可賦值技巧u利用冒號(hào)“：”給全行的元素賦值例：a=eye(5); a(2,:)=ones(1,5) au利用行、列標(biāo)注構(gòu)成新的矩陣?yán)篵=a(2,3,4,1,2,3)零矩陣與空矩陣的區(qū)別？Matlab編程基礎(chǔ)編程基礎(chǔ)利用小矩陣組成大矩陣?yán)篴=ones(2,3) b=zeros(3,2) c=a;a;a d=b;b e=c,d常用生成矩陣的方法u利用生成矩陣的函數(shù)（P17-18）u利用“：”運(yùn)算符冒號(hào)運(yùn)算符創(chuàng)建了具有10個(gè)元素的向量。利用冒號(hào)運(yùn)算符創(chuàng)建向量的基

4、本語法如下：X = J:INC:KMatlab編程基礎(chǔ)編程基礎(chǔ)Matlab內(nèi)部的特殊變量和常數(shù)例：ans、eps、realmax、realmin、pi、Inf 、i、j、NaN復(fù)數(shù)的賦值方式u可將矩陣元素逐個(gè)賦予復(fù)數(shù)；u將矩陣的實(shí)部和虛部分別賦值。Matlab編程基礎(chǔ)編程基礎(chǔ)運(yùn)算符與數(shù)學(xué)表達(dá)式p算術(shù)運(yùn)算符（算術(shù)運(yùn)算符（P20表表2.2）矩陣加減：矩陣加減：“+”、“-”功能：兩矩陣對(duì)應(yīng)元素相加/減。A、B矩陣的階數(shù)必須相等或其中之一為標(biāo)量，標(biāo)量可與任意大小的矩陣相加/減。矩陣相乘：矩陣相乘：“*”矩陣元素相乘：矩陣元素相乘：“.*”功能：矩陣A、B對(duì)應(yīng)元素相乘。A、B行列必相同或其中之一為標(biāo)量

5、。Matlab編程基礎(chǔ)編程基礎(chǔ)矩陣右除矩陣右除B/A與左除與左除AC如果A為一非奇異矩陣，則右除：B/A=B*inv(A)如果A為一非奇異矩陣，則左除：AC=inv(A)*C矩陣元素右除矩陣元素右除A./B與左除與左除A.BA,B兩矩陣必同行列，或其中之一為標(biāo)量。矩陣冪矩陣冪“”及矩陣元素冪及矩陣元素冪“.”矩陣轉(zhuǎn)置矩陣轉(zhuǎn)置“”“”及矩陣非共軛轉(zhuǎn)置及矩陣非共軛轉(zhuǎn)置“.”Matlab編程基礎(chǔ)編程基礎(chǔ)p關(guān)系運(yùn)算符關(guān)系運(yùn)算符、=、=、=Matlab關(guān)系運(yùn)算符用來比較兩個(gè)同樣大小的數(shù)組，或用來比較一個(gè)數(shù)組和一個(gè)標(biāo)量；當(dāng)數(shù)組和標(biāo)量進(jìn)行比較時(shí)，結(jié)果與數(shù)組大小一樣；關(guān)系比較的結(jié)果只有兩種可能，即1或0。例：

6、比較下面兩段語句的結(jié)果x=(-3:3)/3 x=x+(x=0)*epssin(x)./x sin(x)./x例：統(tǒng)計(jì)數(shù)組A中值大于5.5元素的個(gè)數(shù)Matlab編程基礎(chǔ)編程基礎(chǔ)p邏輯運(yùn)算符邏輯運(yùn)算符&(與與) |(或或) (非非)真值表（P24）p其他邏輯函數(shù)（其他邏輯函數(shù)（P25表表2.6,2.7）p數(shù)學(xué)表達(dá)式的數(shù)學(xué)表達(dá)式的MATLAB描述（描述（P25）不要忽略數(shù)乘不能調(diào)用系統(tǒng)中不存在的函數(shù)Matlab編程基礎(chǔ)編程基礎(chǔ)流程控制pif 語句語句格式格式:格式1If 表達(dá)式 語句組AEnd例例1_1:隨機(jī)生成一個(gè)隨機(jī)生成一個(gè)(0,1)之間的數(shù)之間的數(shù),如果該數(shù)大于等如果該數(shù)大于等于于0.

7、5則輸出該數(shù)則輸出該數(shù). Matlab編程基礎(chǔ)編程基礎(chǔ)格式2If 表達(dá)式 語句組AElse 語句組BEnd例1_2: 隨機(jī)生成一個(gè)(0,1)之間的數(shù),如果該數(shù)大于等于0.5,則將該數(shù)變?yōu)?,否則將該數(shù)變成-1.例1_3: 隨機(jī)生成100個(gè)(0,1)之間的數(shù),將大于等于0.5的變?yōu)?,小于0.5的變?yōu)?.Matlab編程基礎(chǔ)編程基礎(chǔ)格式3If 表達(dá)式1 語句組AElseif 表達(dá)式2 語句組BElse 語句組CEnd例1_4: 隨機(jī)生成一個(gè)(0,1)之間的數(shù),如果該數(shù)所處區(qū)間為(0,0.3)則將該數(shù)變成0,如果該數(shù)所處區(qū)間為(0.3,0.6),則將該數(shù)變成1,否則將該數(shù)變成2. Matlab編程

8、基礎(chǔ)編程基礎(chǔ)pSwitch語句語句Switch語句稱為條件選擇語句語句稱為條件選擇語句,其中關(guān)鍵字包括其中關(guān)鍵字包括:switch, case, otherwise以及以及end.主要用于有選擇性的程序設(shè)計(jì)主要用于有選擇性的程序設(shè)計(jì),實(shí)現(xiàn)程實(shí)現(xiàn)程序的多分支選擇序的多分支選擇.格式格式:Switch 選擇表達(dá)式 case 情況表達(dá)式1 語句組1 case情況表達(dá)式2 語句組2 otherwise 語句組nend 例例1_5：隨機(jī)生成一個(gè)：隨機(jī)生成一個(gè)1到到12之間的之間的整數(shù)，按如下規(guī)則進(jìn)行轉(zhuǎn)換：整數(shù)，按如下規(guī)則進(jìn)行轉(zhuǎn)換：如果該數(shù)在如果該數(shù)在14之間，轉(zhuǎn)換成之間，轉(zhuǎn)換成-1；如果該數(shù)在如果該數(shù)在

9、58之間，轉(zhuǎn)換成之間，轉(zhuǎn)換成0；如果該數(shù)在如果該數(shù)在912之間，轉(zhuǎn)換成之間，轉(zhuǎn)換成1。Matlab編程基礎(chǔ)編程基礎(chǔ)pWhile語句語句While語句為條件循環(huán)語句，循環(huán)次數(shù)不確定，取語句為條件循環(huán)語句，循環(huán)次數(shù)不確定，取決于邏輯表達(dá)式。決于邏輯表達(dá)式。格式：格式：While 表達(dá)式 語句組AEndMatlab編程基礎(chǔ)編程基礎(chǔ)pFor語句語句是一種循環(huán)次數(shù)確定的循環(huán)語句。是一種循環(huán)次數(shù)確定的循環(huán)語句。格式：格式： for V=初值：增量：終值初值：增量：終值 語句組語句組A end功能：將語句組功能：將語句組A循環(huán)執(zhí)行循環(huán)執(zhí)行N次，其中：次，其中： N=1+（終值（終值-初值）初值）/增量增量M

10、atlab編程基礎(chǔ)編程基礎(chǔ)pFor循環(huán)注意事項(xiàng)循環(huán)注意事項(xiàng)不能通過在循環(huán)體內(nèi)重新賦值給循環(huán)變量來終止FOR循環(huán),只能使用Break語句;例:for n=1:10 x(n)=sin(n*pi/10); n=10; end for n=1:10 x(n)=sin(n*pi/10); if n=2 break;endMatlab編程基礎(chǔ)編程基礎(chǔ)循環(huán)被執(zhí)行前,應(yīng)預(yù)先分配數(shù)組;當(dāng)有一個(gè)等效的數(shù)組方法來解答給定的問題時(shí),應(yīng)避免使用循環(huán);循環(huán)可按需要進(jìn)行嵌套;例:(前面while循環(huán)的例子)Matlab編程基礎(chǔ)編程基礎(chǔ)例: 隨機(jī)生成01000之間的隨機(jī)數(shù)100個(gè),并按從小到大的順序輸出.例: 隨機(jī)生成511

11、00之間的隨機(jī)數(shù)50個(gè),請(qǐng)輸出序列中不存在的51100之間的數(shù).例: 請(qǐng)驗(yàn)證兩個(gè)正交的復(fù)噪聲信號(hào)之和的包絡(luò)服從Rayleigh分布Matlab編程基礎(chǔ)編程基礎(chǔ)數(shù)據(jù)的輸入/輸出及文件的讀寫(P31)psave, load命令命令基本數(shù)學(xué)函數(shù)(P42)基本繪圖方法p圖形窗口控制命令圖形窗口控制命令figure(N)功能:新創(chuàng)建一個(gè)圖形窗口,并自動(dòng)/N編號(hào)圖形窗口內(nèi)容保持命令hold on (off)功能:保持當(dāng)前圖形窗口的內(nèi)容,或取消hold on功能,恢復(fù)系統(tǒng)默認(rèn)狀態(tài).圖形窗口劃分命令格式: sumplot(m,n,p)功能: 將圖形窗口分成n*m個(gè)子圖形窗口,并選擇第p個(gè)子圖形窗口作為當(dāng)前圖形

12、窗口,供繪圖函數(shù)作圖使用.Matlab編程基礎(chǔ)編程基礎(chǔ)p二維圖形的繪制二維圖形的繪制二維圖形繪制命令plot格式一: plot(y)功能:輸入一個(gè)數(shù)據(jù)的情況. 如果y是一個(gè)數(shù)組,函數(shù)plot(y)給出笛卡兒坐標(biāo)下的二維圖.以y中元素的下標(biāo)為X坐標(biāo),以y中元素值為Y坐標(biāo).例:給出一個(gè)周期的sin(x)的圖形格式二: plot(x,y)功能:以x數(shù)組中元素值為橫坐標(biāo),以y元素值為縱坐標(biāo)作圖.例1_10:給出一個(gè)從(-pi,pi)的cos(x)圖形Matlab編程基礎(chǔ)編程基礎(chǔ)線型,標(biāo)記符及顏色的設(shè)置(P50)格式一:plot(y,字符串)格式二:plot(x,y,字符串)其中:字符串中是關(guān)于線形,顏

13、色及標(biāo)記符號(hào)的說明圖形的標(biāo)注、網(wǎng)格及圖例說明(1)添加圖形標(biāo)題命令title格式: title(string)功能:在當(dāng)前坐標(biāo)系的頂部加一個(gè)文本串string,作為該圖的標(biāo)題.(2)添加坐標(biāo)軸標(biāo)志函數(shù)xlabel,ylabel,zlabel格式:xlabel(string) ylabel(string) zlabel(string)Matlab編程基礎(chǔ)編程基礎(chǔ)(3)設(shè)置網(wǎng)格線命令grid格式: grid on(off)功能: 對(duì)當(dāng)前坐標(biāo)圖加上網(wǎng)格式或撤銷網(wǎng)格線.(4)圖形標(biāo)注函數(shù)legend格式:legend(string1,string2,string3,.)功能:在當(dāng)前圖中添加圖例.p坐標(biāo)

14、軸刻度設(shè)置坐標(biāo)軸刻度設(shè)置格式: axis(xmin,xmax,ymin,ymax,)功能:對(duì)x,y,等坐標(biāo)軸進(jìn)行標(biāo)定.Matlab編程基礎(chǔ)編程基礎(chǔ)p多條曲線的繪制多條曲線的繪制格式:plot(x1,y1,x2,y2,x3,y3,)功能:為每一對(duì)X-Y繪制一條圖線.可以在每個(gè)xi,yi后加線形,標(biāo)記和顏色控制.例1_11:在一張圖中畫出sin(t),cos(t),sin(t)+cos(t)的圖形Matlab編程基礎(chǔ)編程基礎(chǔ)M文件命令語句構(gòu)成的程序文件命令語句構(gòu)成的程序文件,可分為兩種可分為兩種:p主程序文件主程序文件由用戶為解決特定問題,而將命令語句序列放到m文件中而形成的文件主程序文件結(jié)構(gòu)u有

15、關(guān)程序的功能,使用方法等內(nèi)容的注釋部分以”%”開始直至行末的部分為注釋部分u程序主體由若干條命令語句組成,實(shí)現(xiàn)程序設(shè)計(jì)功能,通常用clear, close all開始u主程的運(yùn)行方式例1_12編寫程序畫出sin(6pix),cos(2pix)以及兩函數(shù)之各的曲線.Matlab編程基礎(chǔ)編程基礎(chǔ)p函數(shù)文件結(jié)構(gòu)函數(shù)文件結(jié)構(gòu)函數(shù)定義行function 變量列表=函數(shù)名(參數(shù))H1行注釋行,主要說明該函數(shù)的功能,可以在help命令中顯示.函數(shù)幫助文本注釋行,用來比較詳細(xì)地說明函數(shù)的相關(guān)信息,可在help命令中顯示函數(shù)體用來完成函數(shù)功能的語句序列,并以end語句為結(jié)束標(biāo)志.例例1_13生成四個(gè)長(zhǎng)度為64的

16、序列,分別為s1,s2,s3,s4,對(duì)每個(gè)序列乘上1,-1,j,-j中的任意一個(gè)因子后求和.求所有可能的和序列中峰均功率比最低的序列及相應(yīng)的因子.系統(tǒng)仿真基礎(chǔ)系統(tǒng)仿真基礎(chǔ)通信系統(tǒng)仿真發(fā)射機(jī)接收機(jī)信道干擾發(fā)射機(jī)發(fā)射機(jī)抽樣量化編碼信源編碼信道編碼調(diào)制信源線性分組碼循環(huán)碼BCHRS卷積碼TurboASKFSKPSKMSKGMSKQAM信源編碼抽樣量化編碼信源編碼信源function signal=informationSource(t)function sampleSignal=sampling()function quantizeData=quantizing(sampleSingal)funct

17、ion bitstream=sourceCoding(quantizeData)信源信源function signal=informationSource(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 q

18、uantizeData=quantizing(sampleSingal) dv=2*a/(M-1); quantizeData=sign(sampleSignal).* fix(abs(sampleSignal+dv/2)/dv)*dv; end020406080100120-2.5-2-1.5-1-0.500.511.522.5量 化 電 平 數(shù) 32020406080100120-2.5-2-1.5-1-0.500.511.522.5量 化 電 平 數(shù) 256編碼編碼function bitstream=sourceCoding(quantizeData,dv,M) N=length(qu

19、antizeData); bitstream=; N1=fix(log2(M); for i=1:N NData=fix(abs(quantizeData(i)/dv); bindata=convert2bin(NData,N1); if (quantizeData(i)=0) bindata(1)=0; else bindata(1)=1; end bitstream=bitstream,bindata; endend找到量化電平編號(hào)將該編碼轉(zhuǎn)換成二進(jìn)制加上符號(hào)位處理下一個(gè)量化電平編碼編碼將一個(gè)整數(shù)轉(zhuǎn)換成二進(jìn)制數(shù)據(jù)function bindata=convert2bin(NData,N) b

20、indata=zeros(1,N); i=N; while (NData=0) bindata(i)=rem(NData,2); NData=fix(NData/2); i=i-1; end if (NData(1)=1) NData=0,ones(1,N-1); endend信道編碼信道編碼信道編碼線性分組碼循環(huán)碼BCHRS卷積碼Turbo為方便,使用奇偶校驗(yàn)碼將原比特流,每7位一組,添加校驗(yàn)位1位.采用偶校驗(yàn).該碼不具糾錯(cuò)功能,但具檢一位錯(cuò)的功能.信道編碼信道編碼function datastream=channelcodingEV(bitstream) N=fix(length(bits

21、tream)/7); datastream=; for i=1:N temp=bitstream(i-1)*7+1:i*7); datastream=datastream,temp,rem(sum(temp),2); endend調(diào)制調(diào)制ASKFSKPSKMSKGMSKQAM本例采用2ASK調(diào)制方式,其它調(diào)制方式可以類似進(jìn)行仿真.Bit速率為1000bps;載波頻率為10000Hz;由每一位bit占用時(shí)間:1/1000=1ms;載波的一個(gè)周期所占用時(shí)間:1/10000=0.1ms因此,一個(gè)bit位有10個(gè)完整的載波周期.載波的一個(gè)周期以16點(diǎn)抽樣.即 t=2*pi/16:2*pi/16:2*p

22、i; carrier=cos(t);調(diào)制時(shí):1)如果某位為1,則復(fù)制10次carrier;2)如果某位為0,則復(fù)制10*16個(gè)0.調(diào)制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)=

23、1) transmittedSignal=transmittedSignal,carrier; else transmittedSignal=transmittedSignal,zeros(1,length(carrier); end endend發(fā)射機(jī)發(fā)射機(jī)Transmitterfunction transmittedSignal=transmitter() sampleSignal=sampling() a=2; M=256; quantizeData,dv=quantizing(sampleSignal,a,M); bitstream=sourceCoding(quantizeData,

24、dv,M); datastream=channelcodingEV(bitstream) bitRate=1000; carrierFrequency=10000; N1=16; transmittedSignal=BAskModu (datastream,bitRate,carrierFrequency,N1);end信源編碼信道編碼調(diào)制信道仿真信道仿真本例使用加性高斯白噪聲p直接疊加平均功率為直接疊加平均功率為sig的噪聲信號(hào)的噪聲信號(hào) 如如: signal noise=sqrt(sig2)*1/sqrt(2)*randn(1,length(signal); signal=signal+n

25、oise;p利用awgn函數(shù) signal=sqrt(2)*signal; signal=awgn(signal,20,0)信道仿真信道仿真function signal=gussianChannel(transmittedSignal,sig) N=length(transmittedSignal); noise=sqrt(sig)*1/sqrt(2)*randn(1,N); signal=transmittedSignal+noise;end接收機(jī)設(shè)計(jì)接收機(jī)設(shè)計(jì)譯碼反量化解調(diào)信道譯碼信源譯碼再現(xiàn)線性分組碼譯碼循環(huán)碼譯碼BCH譯碼RS譯碼卷積碼譯碼Turbo譯碼ASK解調(diào)FSK解調(diào)PSK解調(diào)

26、MSK解調(diào)GMSK解調(diào)QAM解調(diào) bitstream=demoASK(receivedSignal,bitRate,carrierFrequency,N);datastream=deEVChannelCoding(bitstream);restoreOriginalSignal=D2A(datastream,dv);ASK解調(diào)解調(diào)在仿真時(shí)假設(shè)能精確同步!ASK解調(diào)解調(diào)BPF整流LPF抽樣判決數(shù)據(jù)速率=1000bps;載波頻率=10000Hz;抽樣頻率=載波頻率*16;信號(hào)0點(diǎn)帶寬(單邊) =2/ =2 /(1/1000) =1000Hz信號(hào)的通帶= 9000Hz,11000Hz根據(jù)信號(hào)速率可得

27、:基帶信號(hào)的帶寬=0,1000Hz假設(shè)所有同步有能精確同步!濾波器設(shè)計(jì)濾波器設(shè)計(jì)幾個(gè)基本概念p通帶、通帶截止頻率、通帶容限通帶、通帶截止頻率、通帶容限p阻帶、阻帶容限、阻帶截止頻率阻帶、阻帶容限、阻帶截止頻率p過渡帶過渡帶p通帶最大衰減通帶最大衰減p阻帶最小衰減阻帶最小衰減濾波器設(shè)計(jì)方法濾波器設(shè)計(jì)方法FDAToolIIR(MATLAB直接設(shè)計(jì))p利用利用IIR濾波器階次估計(jì)函數(shù)求最小階數(shù)和截止濾波器階次估計(jì)函數(shù)求最小階數(shù)和截止頻率（頻率（P162，P163）FIR(MATLAB直接設(shè)計(jì))pfir1pfir2ASK解調(diào)解調(diào)function bitstream=demoASK(receivedSignal,bitRate,carrierFrequency,N) load mypass den=1; signal1=filter(bandpassNum,den,receivedSignal); signal1=