《信息與通信系統(tǒng)仿真》課件6 Matlab信號(hào)處理

1、別志松 Matlab信號(hào)處理Matlab信號(hào)處理現(xiàn)代信息與通信系統(tǒng)仿真中，無(wú)論發(fā)射機(jī)、接收機(jī)的建模還是一些常用的分析方法本質(zhì)上都是在進(jìn)行數(shù)字信號(hào)處理。信號(hào)處理是信息與通信系統(tǒng)建模與仿真的基本組成部分。Matlab中提供了信號(hào)處理工具箱，該工具箱中的函數(shù)與Matlab本身所攜帶的一些內(nèi)置函數(shù)和統(tǒng)計(jì)工具箱中的有關(guān)函數(shù)一起實(shí)現(xiàn)Matlab強(qiáng)大的信號(hào)處理功能。提綱信號(hào)處理基礎(chǔ)線性時(shí)不變系統(tǒng)的表示與分析濾波器的設(shè)計(jì)、分析與實(shí)現(xiàn)隨機(jī)信號(hào)分析和經(jīng)典譜估計(jì)信號(hào)的表示向量表示單通道信號(hào)在時(shí)間軸上的采樣；多通道信號(hào)：矩陣。在這種表示方法中，實(shí)際上是把列號(hào)(或行號(hào))隱含作為對(duì)應(yīng)的采樣時(shí)間，若需要指定信號(hào)的真實(shí)下標(biāo)，

2、則需構(gòu)造另一個(gè)向量作為時(shí)間向量與之對(duì)應(yīng)。若要產(chǎn)生具有特定采樣率的信號(hào)，需要定義時(shí)間軸向量。常見(jiàn)典型信號(hào)的生成t = (0:0.001:1);imp= 1; zeros(99,1); %單位沖擊unit_step = ones(100,1); % 單位階躍(0點(diǎn)跳變.)ramp_sig= t; % 單位斜變quad_sig=t.2; % 二次序列sq_wave = square(4*pi*t); % 周期為0.5的方波函數(shù)名稱波形描述chirp調(diào)頻余弦信號(hào)diric Dirichlet 或周期 sinc 函數(shù)gauspuls高斯調(diào)制的余弦脈沖gmonopuls高斯單脈沖pulstran脈沖串re

3、ctpuls非周期方波sawtooth 鋸齒波或三角波sinc抽樣函數(shù)square方波tripuls非周期的三角波vco壓控振蕩信號(hào)的基本運(yùn)算運(yùn)算類型數(shù)學(xué)描述MALAB實(shí)現(xiàn)信號(hào)加x(n)=x1(n)+x2(n) x=x1+x2信號(hào)乘x(n)=x1(n)x2(n) x=x1. *x2倍率運(yùn)算y(n)kx(n) y=k*x位移y(n)=x(n-k)y=zeros(1,k) x翻轉(zhuǎn)y(n)=x(-n) y=fliplr(x)采樣和y(n)=x(n)sum (x(n1:n2) 采樣積y(n)=x(n)prod(x(n1:n2)具有時(shí)間基的信號(hào)相加如果序列長(zhǎng)度不等，可以通過(guò)補(bǔ)零的方式使x1(n)、x2(

4、n)成為具有相同維數(shù)的序列向量y1(n)、y2(n)，然后對(duì)y1(n)、y2(n)相加或相乘。大部分信號(hào)的表示都是由其時(shí)間序號(hào)和信號(hào)值兩部分組成的，實(shí)現(xiàn)信號(hào)的乘加運(yùn)算的關(guān)鍵在于使兩者的時(shí)間基對(duì)齊。序列平移序列的平移可以看作是將序列的時(shí)間序號(hào)向量平移，而對(duì)應(yīng)原時(shí)間序號(hào)的序列樣值不變。要將序列左移n0個(gè)單位時(shí)，則將時(shí)間序號(hào)向量都減小n0個(gè)單位；若要右移n0個(gè)單位時(shí)，則將時(shí)間序號(hào)向量都增大n0個(gè)單位。對(duì)具有對(duì)應(yīng)時(shí)間向量的信號(hào)而言，序列的平移其實(shí)不需要改變信號(hào)序列本身，只需要改變時(shí)間向量值即可。序列翻轉(zhuǎn)序列的翻轉(zhuǎn)可用MATLAB中的fliplr函數(shù)實(shí)現(xiàn)。n=0:10;x=n.*n-n;y=fliplr

5、(x);n1=-fliplr(n);subplot(1,2,1),stem(n,x);xlabel(n);ylabel(x(n);title(原序列);subplot(1,2,2),stem(n1,y);xlabel(n);ylabel(x(-n);title(翻轉(zhuǎn)序列);離散傅里葉變換運(yùn)用傅里葉分析，可以從頻域的觀點(diǎn)分析和處理信號(hào)，達(dá)到時(shí)域分析不能獲得的效果。離散傅里葉變換DFT是Matlab數(shù)字信號(hào)處理的主要工具之一，由Matlab提供的一組內(nèi)置函數(shù)實(shí)現(xiàn)。這些函數(shù)包括fft、ifft、fftshift和ifftshift。fft和ifft函數(shù)以著名的fftw為核心編寫(xiě)，執(zhí)行效率較高。信號(hào)處

6、理工具箱中的許多函數(shù)，如譜分析、濾波器實(shí)現(xiàn)等功能，都用到了FFT。fft與ifft函數(shù)分別可以計(jì)算DFT及IDFTDFT和IDFT計(jì)算為符合MATLAB的表示方式，上式中的序列下標(biāo)從1開(kāi)始。fft函數(shù)和ifft函數(shù)用法如果x為向量，則函數(shù)fft(x) 計(jì)算向量的DFT;若x為矩陣，則計(jì)算每列的DFT。FFT函數(shù)可以指定FFT點(diǎn)數(shù)N，這種情況下如果向量長(zhǎng)度大于N則截?cái)啵∮贜則補(bǔ)零。ifft與fft函數(shù)的使用方法完全相同，需要注意的是，fft函數(shù)執(zhí)行完畢信號(hào)的平均功率放大N倍，ifft之后信號(hào)平均功率則減小為原來(lái)的N分之一。fftshift在執(zhí)行FFT后，零頻的位置為向量中第一個(gè)數(shù)據(jù)，Matla

7、b提供了fftshift函數(shù)，其作用是將fft之后數(shù)據(jù)的零頻移到數(shù)據(jù)的中間，實(shí)際上就是將向量的左右兩個(gè)半邊交換位置。函數(shù)ifftshift具有類似功能，當(dāng)向量長(zhǎng)度為偶數(shù)時(shí)兩者功能完全相同， 長(zhǎng)度為奇數(shù)時(shí)兩者稍有區(qū)別。fft示例Fs = 1000; % 采樣頻率T = 1/Fs; % 采樣時(shí)間L = 1000; % 信號(hào)長(zhǎng)度t = (0:L-1)*T; % 時(shí)間向量x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); % 信號(hào)：兩個(gè)正弦波疊加y = x + 2*randn(size(t); %加上白噪聲subplot(1,3,1)plot(Fs*t(1:50),y(

8、1:50) %時(shí)域波形title(Signal Corrupted with Zero-Mean Random Noise)xlabel(time (ms)NFFT = 2nextpow2(L); %為了加快計(jì)算速度，選擇FFT點(diǎn)數(shù)為2的冪次。Y = fft(y,NFFT)/L; %進(jìn)行FFT變換f = Fs/2*linspace(0,1,NFFT/2+1); %頻域橫坐標(biāo)subplot(1,3,2)% Plot single-sided amplitude spectrum.plot(f,2*abs(Y(1:NFFT/2+1) %利用其中一半繪制單邊頻譜title(Single-Sided

9、Amplitude Spectrum of y(t)xlabel(Frequency (Hz)ylabel(|Y(f)|)f1= Fs*linspace(-1,1,NFFT)/2;yy=abs(fftshift(Y); %利用fftshift實(shí)現(xiàn)移頻subplot(1,3,3)plot(f1,yy) %繪制移頻后的雙邊頻譜title(Double-Sided Amplitude Spectrum of y(t)xlabel(Frequency (Hz)ylabel(|Y(f)|)信號(hào)的能量數(shù)字復(fù)信號(hào)能量可以通過(guò)對(duì)于序列中所有樣點(diǎn)的能量進(jìn)行累加得到可以選用三個(gè)不同的Matlab語(yǔ)句實(shí)現(xiàn)：sum(

10、x.*conj(x); sum(abs(x).2);x*x;樣點(diǎn)功率數(shù)字復(fù)信號(hào)的信號(hào)功率則總可以與樣點(diǎn)平均功率聯(lián)系起來(lái)，平均樣點(diǎn)功率為sum(x.*conj(x) /length(x); sum(abs(x).2) /length(x);x*x/length(x);多速率信號(hào)處理信號(hào)處理中經(jīng)常遇到過(guò)采樣、降采樣、抽取和插值等處理，這些處理統(tǒng)稱為多速率信號(hào)處理。MATLAB內(nèi)置函數(shù)和信號(hào)處理工具箱提供了一組完整的函數(shù)，能夠完成各種多速率信號(hào)處理任務(wù)。函數(shù)名稱功能描述decimate 抽取，以低通濾波的方式減低采樣率downsample降采樣，以整數(shù)倍降低采樣率interp低通濾波形式的插值，整數(shù)

11、倍內(nèi)插resample任意倍數(shù)的信號(hào)重采樣upfirdn升采樣、FIR濾波及降采樣級(jí)聯(lián)操作upsample升采樣，增加采樣率整數(shù)倍抽取給定有限長(zhǎng)N的離散時(shí)間序列信號(hào)x(n)，n=0,1,N-1，若要將其采樣速率減小M倍 （假定N/M為整數(shù)），最簡(jiǎn)單的方法是將x(n)中每M個(gè)點(diǎn)抽取一個(gè)， 組成新的離散序列 ，然而，根據(jù)采樣定理，當(dāng)新的采樣率低于原信號(hào)頻寬的2倍時(shí)，會(huì)發(fā)生混疊，所以在抽取之前先對(duì)x(n)做低通濾波。MATLAB的decimate函數(shù)實(shí)現(xiàn)信號(hào)的整數(shù)倍抽取，常用的調(diào)用方式為y = decimate(x,r)，在抽取之前默認(rèn)采用8階切比雪夫 I型低通濾波器壓縮頻帶，截止頻率為 0.8*(

12、Fs/2)/r，還可以指定所采用的低通濾波器的類型及階數(shù)，調(diào)用形式可以為y = decimate(x,r,n)或y = decimate(x,r,fir)。插值與之相對(duì)應(yīng)的是interp函數(shù)。MATLAB的函數(shù)interp(x,r)進(jìn)行信號(hào)的整數(shù)倍內(nèi)插，并在插值后進(jìn)行低通濾波。y = interp(x,r) 將 x 的采樣率提高 r倍，插值后的向量y比原始輸入x長(zhǎng) r倍。還可以用y = interp(x,r,l,alpha)指定濾波器長(zhǎng)度l及截止頻率因子alpha。默認(rèn)情況下l =4， alpha = 0.5。resampleMatlab還提供了非整數(shù)倍的插值和抽取函數(shù)，對(duì)于給定的有限長(zhǎng)N的離

13、散時(shí)間序列信號(hào)x(n)，n=0,1,N-1，若希望將原來(lái)的采樣率Fs變?yōu)樵瓉?lái)的L/M倍， 則可以先對(duì)原信號(hào)做L倍插值，然后做M倍抽取，MATLAB的函數(shù)resample實(shí)現(xiàn)這一功能。resample(data,P,Q) 重采樣data，先插值 P 倍，然后再按Q倍抽取，默認(rèn)采用10階濾波器實(shí)現(xiàn)，也可通過(guò)調(diào)用方式 resample(data,P,Q,order) 指定了插值之后和抽取之前所用濾波器的階數(shù)。線性時(shí)不變系統(tǒng)的表示與分析線性時(shí)不變系統(tǒng)的基本模型線性系統(tǒng)是指系統(tǒng)滿足均勻性和疊加性，若輸入x1(n)與 x2(n)的 輸出分別為y1(n)和 y2(n)，線性系統(tǒng)對(duì)輸入ax1(n)+bx2(n

14、)的輸出為 ay1(n)+by2(n)。時(shí)不變是指系統(tǒng)響應(yīng)與激勵(lì)加于系統(tǒng)的時(shí)刻無(wú)關(guān)。若輸入x1(n)的 輸出為y1(n)，則輸入x(n-m)時(shí)輸出y(n-m)。MATLAB的信號(hào)處理工具箱提供了多種線性系統(tǒng)的表示模型，包括傳遞函數(shù)、極點(diǎn)增益法、狀態(tài)空間法、帶余數(shù)的部分分式展開(kāi)法等，以方便用戶針對(duì)具體應(yīng)用選擇最適合的表示方法，并提供了這些模型之間的轉(zhuǎn)換函數(shù)。傳遞函數(shù)表示法其中常量b(i) 和 a(i)是濾波器系數(shù)，n和m 的最大值是濾波器的階數(shù)。在Matlab中，可以用兩個(gè)向量存儲(chǔ)系數(shù)（行向量），一個(gè)向量存儲(chǔ)分子，一個(gè)存儲(chǔ)分母，這兩個(gè)向量實(shí)質(zhì)上對(duì)應(yīng)的是分子和分母兩個(gè)多項(xiàng)式。傳遞函數(shù)表示法舉例例如

15、傳遞函數(shù)為H(z)=(z2-0.5z+2)/(z2+0.4z+1)的線性系統(tǒng)用兩個(gè)向量表示： b=1 -0.5 2 a=1 0.4 1調(diào)用matlab的函數(shù)，可計(jì)算該系統(tǒng)的頻率響應(yīng)： h,f=freqz(b,a)plot(f,abs(h); %繪出幅頻特性；零極點(diǎn)增益表示法仍然是傳遞函數(shù)不過(guò)形式不同狀態(tài)空間表示法數(shù)字濾波器或數(shù)字系統(tǒng)通?？梢员硎緸橐浑A差分方程組的形式，組成一組狀態(tài)方程和輸出方程多個(gè)模型之間的轉(zhuǎn)換函數(shù)matlab工具箱提供了多個(gè)線性系統(tǒng)模型之間相互轉(zhuǎn)換的函數(shù)。tf2ss 由傳遞函數(shù)模型轉(zhuǎn)換成狀態(tài)空間模型，Ss2tf由狀態(tài)空間模型轉(zhuǎn)換成傳遞函數(shù)模型；tf2zp由傳遞函數(shù)模型轉(zhuǎn)換成零

16、極點(diǎn)增益模型，zp2ft由零極點(diǎn)增益模型轉(zhuǎn)換成傳遞函數(shù)模型；ss2zp 由狀態(tài)空間模型轉(zhuǎn)換成零極點(diǎn)增益模型，zp2ss由零極點(diǎn)增益模型轉(zhuǎn)換成狀態(tài)空間模型等。卷積任意序列都可以表示為單位抽樣序列的移位加權(quán)和形式如果已知線性時(shí)不變系統(tǒng)的沖激響應(yīng)h(n)，則對(duì)任意輸入x(n)的輸出y(n)均可以用卷積實(shí)現(xiàn)：y(n)=x(n)*h(n)Matlab中卷積的實(shí)現(xiàn)在Matlab中，如果已知線性系統(tǒng)的任何表示方式，都可以得到其沖激響應(yīng)，例如，對(duì)單極點(diǎn)濾波器b = 1, a = 1 -0.9，可以用h = filter(b,a,imp)得到其沖激響應(yīng)，由零極點(diǎn)模型、狀態(tài)空間模型等也能夠通過(guò)模型變換得到其沖激響

17、應(yīng)。如果已知沖激響應(yīng)，則可以利用MATLAB提供的conv函數(shù)執(zhí)行標(biāo)準(zhǔn)的一維卷積運(yùn)算，conv函數(shù)的調(diào)用格式為：w = conv(u,v)，如果向量u的長(zhǎng)度為m，向量v的長(zhǎng)度為n ，則w為一個(gè)長(zhǎng)度為m+n-1的向量，其第 k個(gè)元素為：帶有時(shí)間序號(hào)的卷積函數(shù)conv假定u(k)、v(k)都是從k=0開(kāi)始，這就限制了它的應(yīng)用范圍。因此，要對(duì)從任意時(shí)間序號(hào)值的序列進(jìn)行卷積運(yùn)算，且正確表示出計(jì)算結(jié)果，還須構(gòu)造序列對(duì)應(yīng)序號(hào)向量。functiony,ny=convsubscrip(x,nx,h,nh)% modified convolution function% convsubscrip(x,nx,h,

18、nh)ny1=nx(1)+nh(1);ny2=nx(length(x)+nh(length(h);ny=ny1:ny2;y=conv(x,h); MATLAB中線性時(shí)不變系統(tǒng)分析實(shí)例若離散線性時(shí)不變系統(tǒng)的單位沖激響應(yīng)為：h(n)=(0.8)nu(n)-u(n-20)輸入的離散信號(hào)為：X(n)= (n)+2 (n-20)+4 (n-40)利用卷積法求解輸出y(n)線性時(shí)不變系統(tǒng)的頻域分析線性時(shí)不變系統(tǒng)的頻域分析與時(shí)域分析之間具有對(duì)偶關(guān)系，是一種經(jīng)常遇到的系統(tǒng)分析方法。輸入信號(hào)通過(guò)線性時(shí)不變系統(tǒng)，所得到的輸出信號(hào)為根據(jù)卷積定理，時(shí)域卷積對(duì)應(yīng)頻域相乘Matlab線性系統(tǒng)頻域分析實(shí)例輸入信號(hào)x(t)1

19、）確定該信號(hào)的離散傅立葉變換2）假如該信號(hào)通過(guò)一個(gè)帶寬為1.5Hz的理想低通濾波器，確定濾波器的輸出并繪出曲線。2）假如該信號(hào)通過(guò)沖激響應(yīng)為h(t)的濾波器，求輸出信號(hào)，其中濾波器的設(shè)計(jì)流程濾波器設(shè)計(jì)圖形界面工具fdatool隨機(jī)信號(hào)分析與經(jīng)典譜估計(jì)我們通常假設(shè)所處理的隨機(jī)過(guò)程為廣義平穩(wěn)隨機(jī)過(guò)程，許多統(tǒng)計(jì)量都采用時(shí)間平均來(lái)代替統(tǒng)計(jì)平均。對(duì)平穩(wěn)隨機(jī)變量序列X(n)，可以用以下特征量描述其概率分布特性。1) 均值：即數(shù)學(xué)期望值，可以用MATLAB函數(shù) mean(x)進(jìn)行計(jì)算。2) 均方值：是隨機(jī)變量x(n)平方的數(shù)學(xué)期望，實(shí)際上等效于信號(hào)的樣點(diǎn)功率，可以用mean(x.*x)得到。3) 方差：計(jì)算

20、函數(shù)為 var(x)4) 標(biāo)準(zhǔn)差：可以采用std(x) 或sqrt(var(x)得到xcorr相關(guān)函數(shù)xcorr函數(shù)用于估計(jì)隨機(jī)序列信號(hào)的互相關(guān)序列，其自相關(guān)序列作為一種特殊情況進(jìn)行處理。隨機(jī)過(guò)程的互相關(guān)本來(lái)應(yīng)當(dāng)計(jì)算：但是由于序列長(zhǎng)度為有限值，因此只能計(jì)算估計(jì)值。不同調(diào)用方式所采用的估計(jì)公式有所不同。xcorr調(diào)用方式c=xcorr(x,y,maxlags,option)option指定互相關(guān)的歸一化選項(xiàng)，它可以是：biased：計(jì)算互相關(guān)函數(shù)的有偏互相關(guān)估計(jì)；unbiased：計(jì)算互相關(guān)函數(shù)的無(wú)偏互相關(guān)估計(jì)；coeff：歸一化方式；none:默認(rèn)狀態(tài)，函數(shù)執(zhí)行非歸一化計(jì)算相關(guān)。協(xié)方差MATLAB信號(hào)處理工具箱提供了xcov用于估計(jì)自協(xié)方差序列與互協(xié)方差序列，其用法與xcorr函數(shù)基本類似。譜分析對(duì)于一個(gè)隨機(jī)信號(hào)，它本身的傅里葉