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

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

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

3、ctpuls非周期方波sawtooth 鋸齒波或三角波sinc抽樣函數(shù)square方波tripuls非周期的三角波vco壓控振蕩信號的基本運算運算類型數(shù)學描述MALAB實現(xiàn)信號加x(n)=x1(n)+x2(n) x=x1+x2信號乘x(n)=x1(n)x2(n) x=x1. *x2倍率運算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)具有時間基的信號相加如果序列長度不等，可以通過補零的方式使x1(n)、x2(

4、n)成為具有相同維數(shù)的序列向量y1(n)、y2(n)，然后對y1(n)、y2(n)相加或相乘。大部分信號的表示都是由其時間序號和信號值兩部分組成的，實現(xiàn)信號的乘加運算的關鍵在于使兩者的時間基對齊。序列平移序列的平移可以看作是將序列的時間序號向量平移，而對應原時間序號的序列樣值不變。要將序列左移n0個單位時，則將時間序號向量都減小n0個單位；若要右移n0個單位時，則將時間序號向量都增大n0個單位。對具有對應時間向量的信號而言，序列的平移其實不需要改變信號序列本身，只需要改變時間向量值即可。序列翻轉(zhuǎn)序列的翻轉(zhuǎn)可用MATLAB中的fliplr函數(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)序列);離散傅里葉變換運用傅里葉分析，可以從頻域的觀點分析和處理信號，達到時域分析不能獲得的效果。離散傅里葉變換DFT是Matlab數(shù)字信號處理的主要工具之一，由Matlab提供的一組內(nèi)置函數(shù)實現(xiàn)。這些函數(shù)包括fft、ifft、fftshift和ifftshift。fft和ifft函數(shù)以著名的fftw為核心編寫，執(zhí)行效率較高。信號處

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

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

8、1:50) %時域波形title(Signal Corrupted with Zero-Mean Random Noise)xlabel(time (ms)NFFT = 2nextpow2(L); %為了加快計算速度，選擇FFT點數(shù)為2的冪次。Y = fft(y,NFFT)/L; %進行FFT變換f = Fs/2*linspace(0,1,NFFT/2+1); %頻域橫坐標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實現(xiàn)移頻subplot(1,3,3)plot(f1,yy) %繪制移頻后的雙邊頻譜title(Double-Sided Amplitude Spectrum of y(t)xlabel(Frequency (Hz)ylabel(|Y(f)|)信號的能量數(shù)字復信號能量可以通過對于序列中所有樣點的能量進行累加得到可以選用三個不同的Matlab語句實現(xiàn)：sum(

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

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

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

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

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

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

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

17、應。如果已知沖激響應，則可以利用MATLAB提供的conv函數(shù)執(zhí)行標準的一維卷積運算，conv函數(shù)的調(diào)用格式為：w = conv(u,v)，如果向量u的長度為m，向量v的長度為n ，則w為一個長度為m+n-1的向量，其第 k個元素為：帶有時間序號的卷積函數(shù)conv假定u(k)、v(k)都是從k=0開始，這就限制了它的應用范圍。因此，要對從任意時間序號值的序列進行卷積運算，且正確表示出計算結(jié)果，還須構(gòu)造序列對應序號向量。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中線性時不變系統(tǒng)分析實例若離散線性時不變系統(tǒng)的單位沖激響應為：h(n)=(0.8)nu(n)-u(n-20)輸入的離散信號為：X(n)= (n)+2 (n-20)+4 (n-40)利用卷積法求解輸出y(n)線性時不變系統(tǒng)的頻域分析線性時不變系統(tǒng)的頻域分析與時域分析之間具有對偶關系，是一種經(jīng)常遇到的系統(tǒng)分析方法。輸入信號通過線性時不變系統(tǒng)，所得到的輸出信號為根據(jù)卷積定理，時域卷積對應頻域相乘Matlab線性系統(tǒng)頻域分析實例輸入信號x(t)1

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

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