窄帶隨機(jī)信號分析與仿真

1、隨機(jī)信號實(shí)驗(yàn)實(shí)驗(yàn)報告題 目 窄帶隨機(jī)信號的仿真與分析實(shí)驗(yàn)十三 窄帶隨機(jī)信號的仿真與分析實(shí)驗(yàn)?zāi)康睦糜嬎銠C(jī)仿真窄帶隨機(jī)信號，考察其數(shù)字特征，以此加深對窄帶隨機(jī)信號的特點(diǎn)及分析方法。熟悉常用的信號處理仿真軟件平臺：matlab、c/c+。 實(shí)驗(yàn)原理如果一個隨機(jī)過程的功率譜密度，只分布在高頻載波附近的一個窄頻率范圍內(nèi)，在此范圍之外全為零，則稱之為窄帶過程。窄帶過程是在信息傳輸系統(tǒng)，特別是接收機(jī)中經(jīng)常遇到的隨機(jī)信號。當(dāng)窄帶系統(tǒng)(接收機(jī))的輸入噪聲(如熱噪聲)的功率譜分布在足夠?qū)挼念l帶(相對于接收機(jī)帶寬)上時，系統(tǒng)的輸出即為窄帶過程。窄帶信號的確切定義如下：一個實(shí)平穩(wěn)隨機(jī)x(t)，如果它的功率譜密度sx

2、()具有下述性質(zhì) 而且?guī)挐M足，則稱此過程為窄帶平穩(wěn)隨機(jī)過程。窄帶平穩(wěn)隨機(jī)過程的功率譜密度函數(shù)如圖24所示。 圖24 窄帶隨機(jī)過程的功率譜密度從示波器觀看窄帶隨機(jī)過程的一個樣本函數(shù)，可看到如圖25所示的波形，從這個波形可以看出，窄帶隨機(jī)過程可以表示成具有角頻率以及慢變幅度與相位的正弦振蕩，這就是說可以把它寫成 z(t)=b(t)cost + (t) 式中，b(t)是隨機(jī)過程的慢變幅度，(t)是過程的慢變相位。稱之為準(zhǔn)正弦振蕩，這就是窄帶過程的數(shù)學(xué)模型。 圖25 窄帶過程的某個樣本函數(shù) 實(shí)驗(yàn)任務(wù)與要求 用matlab或c/c+語言編寫和仿真程序。 產(chǎn)生滿足下列條件的窄帶隨機(jī)信號： 其中a(t)包

3、絡(luò)頻率為1khz,幅值為1v。載波頻率為：4khz，幅值為1v，是一個固定相位，n(t)為高斯白噪聲，采樣頻率設(shè)為16khz。實(shí)際上，這是一個帶有載波的雙邊帶調(diào)制信號，如下圖：計算窄帶隨機(jī)信號的均值、均方值、方差、概率密度、頻譜及功率譜密度、相關(guān)函數(shù)，用圖示法來表示。 窄帶系統(tǒng)檢測框圖如圖26所示： 圖26 窄帶系統(tǒng)檢測 x(t)：窄帶隨機(jī)信號、：載波 低通濾波器設(shè)計：低通濾波器技術(shù)要求：通帶截止頻率1khz阻帶截止頻率2khz。過渡帶:1khz阻帶衰減：35db通帶衰減：1db采樣頻率：44.1khz 計算a點(diǎn)、b點(diǎn)、ac(t)、as(t)、y(t)的均值、均方值、方差、頻譜及功率譜密度、相

4、關(guān)函數(shù)，用圖示法來表示。下圖為ac(t)的時域、頻域波形： 注意觀察、在系統(tǒng)中起到什么作用，調(diào)制信號的頻譜結(jié)構(gòu)有沒有變化。 y(t)還原了什么信號。 按要求寫實(shí)驗(yàn)報告。實(shí)驗(yàn)步驟：1. 根據(jù)實(shí)驗(yàn)要求得到x(t)的表達(dá)式：x=(1+a).*cos(2*pi*4000*t+2)+noisy/10;其中：noisy為高斯白噪聲，由wgn函數(shù)生成a=cos(2*pi*1000*t)均值：ex=mean(x);方差：dx=var(x);可得：x(t)的均值為3.1440e-004x(t)的方差為0.7667x(t)的時域波形，均值還有方差如圖01所示：圖01-x(t)的波形（藍(lán)色），方差（上方的直線），均

5、值（下方的直線）2. 求x(t)的概率密度方法是將最大最小區(qū)間分成14等份,然后分別計算各個區(qū)間的個數(shù)，如圖02中柱形條所示，利用曲線擬合，得到合適的概率密度函數(shù)。為了得到光滑的曲線，我們利用了多項(xiàng)式擬合，經(jīng)過測試，9次擬合曲線比較符合我們的要求，獲得的曲線如圖02中曲線所示：圖02-x(t)概率密度分布3. 對x(t)進(jìn)行頻譜分析在matlab中，利用fft函數(shù)可以很方便得求得x(t)的頻譜，然后用abs和angle函數(shù)求得幅值和相位，畫出圖像如圖03所示：圖03-x(t)的頻譜分析圖4. 求x(t)的自相關(guān)函數(shù)我們的方法是用xcorr函數(shù)求出自相關(guān)序列，得到自相關(guān)函數(shù)的圖像如圖04所示：圖

6、04-x(t)的自相關(guān)函數(shù)圖像5. 求x(t)功率譜密度有了自相關(guān)函數(shù)，對它進(jìn)行fft變換，可得x(t)的功率譜密度，如圖05所示：圖05-x(t)的功率譜密度6. 建立濾波器下面我們建立一個巴特沃思濾波器對產(chǎn)生的x(t)進(jìn)行檢測。濾波器的幅度譜和相位譜如圖06，圖07所示：圖06-巴特沃思濾波器幅度譜圖07-巴特沃思濾波器相位譜7. 求ac的統(tǒng)計特性ac為x(t)*2*cos(2*4000*t)通過步驟6中建立的低通濾波器的信號，其各項(xiàng)統(tǒng)計特性如下：ac輸出信號均值eh=-0.0024ac的均方值8.8873e-004ac的方差8.8383e-004各圖像如下列圖所示：圖08-ac(t)的時

7、域圖圖09-ac(t)的頻譜圖圖10-ac(t)的自相關(guān)函數(shù)圖11-ac(t)的功率譜密度8. 求as的統(tǒng)計特性as為x(t)*2*sin(2*4000*t)通過步驟6中建立的低通濾波器的信號，其各項(xiàng)統(tǒng)計特性如下：as輸出信號均值eh= -0.0018as的均方值0.0039as的方差0.0039各圖像如下列圖所示：圖12-as(t)的時域圖圖13-as(t)的頻譜圖圖14-as(t)的自相關(guān)函數(shù)圖15-as(t)的功率譜密度9. 求輸出y(t)的統(tǒng)計特性最后，我們來求y(t)的統(tǒng)計特性,y(t)的表達(dá)式如下：y(t)=ac(t)*cos(2*4000*t)+as(t)*sin(2*4000*

8、t)所用的方法和前述一樣，其統(tǒng)計特性如下y(t)的均值：0.1325y(t)的均方值：0.0720y(t)的方差：0.0544各圖像如下列圖所示：圖16-y(t)的時域圖圖17-y(t)的頻譜圖圖18-y(t)的自相關(guān)函數(shù)圖19-y(t)的功率譜密度附件1：心得體會終于在小組的努力下完成了這個我們大學(xué)真正獨(dú)立完成的隨機(jī)信號實(shí)驗(yàn)，辛勤的汗水換回那一張一張圖片是那么的讓人賞心悅目。雖然這不是我們的第一次合作，但合作的重要性又一次讓我們明白其價值的不菲。對于我們這些只知道念書讀書的大學(xué)生來說，這次實(shí)驗(yàn)的意義遠(yuǎn)不止完成這個實(shí)驗(yàn)，更重要是讓我們明白了知識與實(shí)踐的橋梁不是我們一直認(rèn)為的那么簡單，只有在那一

9、次次的調(diào)試我們自己的程序，才會明白一切的一切都是需要我們?nèi)ブ匦聦W(xué)習(xí)。兩年的時光足以忘記很多東西，大一學(xué)習(xí)的matlab再次回來自己的手邊是那么陌生，自以為學(xué)會編程的我們面對那一個一個的錯誤才會明白知識的無止境。實(shí)驗(yàn)中存在的問題：1 對于matlab太陌生，沒有系統(tǒng)的學(xué)習(xí)過那些真正用于信號分析的那些工具箱。2 對于隨機(jī)信號的知識書本上的那些知識，沒有學(xué)寬，學(xué)出去。3 對某些小的知識點(diǎn)了解不夠問題解決方法：1. 合理的分工，大大提高了工作效率。2. 合理的利用可以使用的一切資源，對于遇到的問題集體討論并上網(wǎng)尋求最佳的答案。3. 對于無法解決的問題，進(jìn)一步咨詢了老師，爭取做到答案的完美。在三人的共同

10、努力下，我們完成了這次實(shí)驗(yàn)，希望老師驗(yàn)收并加以指正。附件2：程序% 說明：此程序是隨機(jī)實(shí)驗(yàn)十三 窄帶隨機(jī)信號的仿真與分析的程序 % 班級：010852班% 作者：鄧曉光 劉春燕 譚正中 %- 生成x(t)并求得x(t)的統(tǒng)計特性 -%clcfs=16000; %設(shè)定采樣頻率n=1300; %取的樣本點(diǎn)數(shù)n=0:n-1; %取樣點(diǎn)序列l(wèi)ength(n)t=n/fs; %獲得以1/16000為時間間隔采樣序列noisy=wgn(1,n,0); %產(chǎn)生高斯白噪聲a=cos(2*pi*1000*t); %獲取a(t)的采樣點(diǎn)x=(1+a).*cos(2*pi*4000*t+2)+noisy/10; %

11、獲取x(t)的采樣點(diǎn)figure(1); %畫圖subplot(2,1,1);plot(n,x);%以t為橫坐標(biāo)畫出x(t)的時域圖型axis(0 140 -3 3);%axis tight; xlabel(采樣點(diǎn));ylabel(x(t)/v);title(圖1-窄帶隨機(jī)信號波形);grid on;%求x(t)均值disp(x(t)的均值為);ex=mean(x);hold on;plot(n,ex,r.);disp(ex);%求x(t)方差disp(x(t)的方差為)dx=var(x);hold on;plot(n,dx,g.)disp(dx);each=linspace(min(x),m

12、ax(x),14);%將最大最小區(qū)間分成19等份,然后分別計算各個區(qū)間的個數(shù)nr=hist(x,each);%計算各個區(qū)間的個數(shù)nr=nr/length(x);%計算各個區(qū)間的個數(shù)歸一化subplot(2,1,2);%多項(xiàng)式擬合bar(each,nr);hold on;plot(each,nr,g)hold on;p=polyfit(each,nr,9);eachi=-2:0.1:2;nri=polyval(p,eachi);plot(eachi,nri,r)axis tight; title(x(t)概率密度分布);xlabel(x(t);ylabel(p(x);grid on;%對x(t)

13、進(jìn)行頻譜分析fx=fft(x,n); %對x(t)進(jìn)行fft變換,在016000區(qū)間內(nèi)得到2n-1個頻率值magn=abs(fx); %求x(t)幅值xangle=angle(fx);%求x(t)相位labelang=(0:length(fx)-1)*16000/length(x);%在016000區(qū)間內(nèi)求橫坐標(biāo)刻度figure(2); %subplot(2,2,3);plot(labelang,magn*10,labelang,xangle);%在016000區(qū)間內(nèi)做頻譜和相位圖%plot(labelang,magn);%在02區(qū)間內(nèi)做頻譜圖axis(0 16000 -0.5 600) ;

14、xlabel(頻率/hz);ylabel(幅值 );title(x(t)頻譜圖);grid on;%求x(t)的自相關(guān)函數(shù)c,lags=xcorr(x,coeff);%求出自相關(guān)序列figure(3); subplot(2,1,1);plot(lags/fs,c); %在時域內(nèi)畫自相關(guān)函數(shù)%axis(-0.003 0.003 -0.5 1);axis tight; xlabel();ylabel( rx();title(x(t)的自相關(guān)函數(shù));grid on;%求x(t)的功率譜密度long=length(c);sx=fft(c,long); %求出自相關(guān)序列 index=0:round(n/

15、2-1);labelx=(0:long-1)*2*pi/long;plot_magn=10*log10(abs(sx);subplot(2,1,2);plot(labelx,plot_magn); %畫功率譜密度axis tight; xlabel(w);ylabel( s(w);title(x(t)的功率譜密度);grid on;%窄帶系統(tǒng)檢測z1=x.*cos(2*pi*4000*t);z2=x.*sin(2*pi*4000*t);ac=z1.*x;as=z2.*x;y=ac.*cos(2*pi*4000*t)+as.*sin(2*pi*4000*t);f_p=1000;f_s=1600;

16、r_p=1;r_s=35;%設(shè)定濾波器參數(shù):通、阻帶截止頻率，通、阻帶衰減ws=2*f_s/fs;wp=2*f_p/fs; %頻率歸一化n,wn=buttord(wp,ws,r_p,r_s);%采用巴特沃思濾波器b,a=butter(n,wn);%求得濾波器傳輸函數(shù)的多項(xiàng)式系數(shù)figure(4);h,w=freqz(b,a);%求得濾波器傳輸函數(shù)的幅頻特性subplot(2,1,1);plot(w*fs/(2*pi),abs(h);%在 02 區(qū)間內(nèi)作幅度譜title(低通濾波器幅度譜)grid on;subplot(2,1,2);plot(w*fs/(2*pi),angle(h);%在 02

17、 區(qū)間內(nèi)作相位譜title(低通濾波器相位譜);grid on;%- 生成ac(t)并求得ac(t)的統(tǒng)計特性 -%信號通過濾波器mc=filter(b,a,ac);%對信號進(jìn)行濾波%求輸出信號均值disp(ac輸出信號均值);eh=mean(mc)%求輸出信號的均方值disp(ac的均方值是);e2h=mc*mc/n%求輸出信號方差disp(ac輸出信號方差)dh=var(mc)%作輸出信號時域波形figure(6)subplot(2,1,1);n=0:n-1;plot(n,mc);%作信號時域波形axis(0 300 -0.1 0.1);% axis tight; xlabel(采樣點(diǎn) )

18、;ylabel(幅值) title(ac時域譜);grid on;%進(jìn)行 fft 變換并做頻譜圖yc=fft(mc,length(mc);%對濾波后信號進(jìn)行 fft 變換longc=length(yc);%求傅里葉變換后的序列長度labelx=(0:longc-1)*16000/longc;magn1=abs(yc);% 求濾波后信號幅值subplot(2,1,2);plot(labelx,magn1);%做頻譜圖axis tight; xlabel(頻率(hz); ylabel(幅值 ); title(ac頻譜圖);grid on;%求輸出信號的自相關(guān)函數(shù)c1,lags1=xcorr(mc,

19、coeff);%求出自相關(guān)序列figure(7);subplot(2,1,1);plot(lags1/fs,c1);%在時域內(nèi)畫自相關(guān)函數(shù)xlabel(); ylabel(rx(); axis tight;title(ac的自相關(guān)函數(shù));grid on;%求輸出信號雙邊功率譜sac=fft(c1,length(c1);%對自相關(guān)函數(shù)進(jìn)行傅里葉變換magnc=abs(sac);%求雙邊功率譜幅值subplot(2,1,2);long=length(sac);%求傅里葉變換后的序列長度labelc=(0:long-1)*16000/long;plot(labelc,10*log10(magnc);

20、%畫自相關(guān)函數(shù)譜xlabel(頻率(hz);ylabel(功率譜(dbw);axis tight;title(ac雙邊功率譜);grid on;%- 生成as(t)并求得as(t)的統(tǒng)計特性 -%信號通過濾波器ms=filter(b,a,as);%對信號進(jìn)行濾波%求輸出信號均值disp(as輸出信號均值);eh=mean(ms)%求輸出信號的均方值disp(as的均方值是);e2h=ms*ms/n%求輸出信號方差disp(as輸出信號方差)dh=var(ms)%作輸出信號時域波形figure(8)subplot(2,1,1);n=0:n-1;plot(n,ms);%作信號時域波形axis(0

21、300 -0.1 0.1);% axis tight; xlabel(采樣點(diǎn) );ylabel(幅值) title(as時域譜);grid on;%進(jìn)行 fft 變換并做頻譜圖ys=fft(ms,length(ms);%對濾波后信號進(jìn)行 fft 變換longs=length(ys);%求傅里葉變換后的序列長度labelx=(0:longs-1)*16000/longs;magn2=abs(ys);% 求濾波后信號幅值subplot(2,1,2);plot(labelx,magn2);%做頻譜圖axis tight; xlabel(頻率(hz); ylabel(幅值 ); title(as頻譜圖

22、);grid on;%求輸出信號的自相關(guān)函數(shù)c2,lags2=xcorr(ms,coeff);%求出自相關(guān)序列figure(9);subplot(2,1,1);plot(lags2/fs,c2);%在時域內(nèi)畫自相關(guān)函數(shù)xlabel(); ylabel(rx(); axis tight;title(as的自相關(guān)函數(shù));grid on;%求輸出信號雙邊功率譜sas=fft(c2,length(c2);%對自相關(guān)函數(shù)進(jìn)行傅里葉變換magnc=abs(sac);%求雙邊功率譜幅值subplot(2,1,2);long=length(sas);%求傅里葉變換后的序列長度labels=(0:long-1)*16000/long;plot(labelc,10*log10(magnc);%畫自相關(guān)函數(shù)譜xlabel(頻率(hz);ylabel(功率譜(dbw);axis tight;title(as雙邊功率譜);grid on;%- 生成y(t)并求得y(t)的統(tǒng)計特性 -%求輸出信號均值disp(y輸出信號均值);eh=mean(y)%求輸出信號的均方值disp(y的均方值是);e2h=y*y/n%求輸出信號方差disp(y輸出信號方差)