第一篇：隨機信號處理案例——雙耳時間差的計算（小編推薦）

隨機信號案例——相關法計算雙耳時間差ITD

1.人耳對聲源的定位

在自然聽音中，人的聽覺系統對聲源的定位取決于多個因素——雙耳接收到的號差異用來決定聲源的水平位置，由外耳對高頻信號的反射所引起的耳郭效應決定聲源的垂直位置，而人耳的某些心理聲學特性對于聲源的定位也起到很大的作用。2.雙耳效應

在自然聽音環境中，雙耳信號之間的差異對于聲源的定位是非常重要的。該因素可以在直達聲場的聽音環境中得到最好解釋，如圖2-1所示。

圖2-1 聲源S與鏡像聲源S′引入最大程度相似的雙耳因素 聲源位于水平面上，水平方位角為θ，與人頭中心的距離為r，到達左右耳的距離分別為SL和SR。由于SL>SR，聲音首先到達右耳，從而在到達雙耳的時間先后上形成時間差。這種時間差被定義為雙耳時間差（interaural time difference，ITD），它與聲源的水平方位角θ有關。當θ = 0°時，= 0；當θ = ±90°時，達到最大值，對一般人頭來說，為0.6～0.7ms 的數量級。在低中頻（f <1.5kHz）情況下，雙耳時間差是定位的主要因素。3.頭相關傳輸函數簡介

頭相關傳輸函數(head-related transfer function, HRTF)描述了自由場聲波從聲源到雙耳的傳輸過程，它反映了頭部、耳廓和軀干等構成的生理系統對聲波散射(綜合濾波)的結果。HRTF 是聲源方向、距離、頻率的連續函數, 它是聲源到雙耳的頻域傳輸函數。自由場的情況下，HRTF 定義為HL?HL?r,?,?,f?，HR?HR?r,?,?,f?，其中r為聲源到頭中心的距離，f為聲波的頻率；方位角 0°≤θ< 360°和仰角?90°≤?≤ 90°表示聲源的方向, 其中φ = 0°和90°分別表示水平面和正上方, 而(θ= 0°，φ= 0°)和(θ= 90°，φ= 0°)分別表示水平面上正前和正右方向。HRTF的時域表示是頭相關脈沖響應hl?r,?,?,t?,和hr?r,?,?,t?簡記為 HRIR，它們與 HL, HR 互為 Fourier 變換。

4.ITD的相關法定義

ITD的定義四[2]（相關法）雙耳脈沖響應HRIR的歸一化互相關函數定義為：

+???LR（?）=????hL(t??)hR(t)dt??1/ ??22t）dt]??[?hL(t)dt][?h（R??????

（3-2-7）

（?）（?）（?）按定義，0≤|?LR|≤1。由式（3-2-7）可計算出函數?LR在|?LR|在|?|≤1ms范圍內的最大值，與此相應的?=?max即為相關法定義的雙耳時間差ITDcorre，即ITDmax??,????max

因而相關法是利用左、右耳HRIR的相似性求出ITD。實際中通常得

hnhn到的是經過離散時間采樣的HRIR，即L??和R??。因而（3-2-7）對連續時間t的積分將變成對離散時間n的求和。例如在44.1KHz的采樣率下，時間分辨率約為23?s。為了提高時間分辨率，在進行

hnhn（3-2-7）計算之前，可先對L??和R??進行過采樣處理。例如10倍過采樣可將時間分辨率變為2.3?s。下面圖a[1]是有26名女性受試者的平均ITD。

圖 a

不同緯度面φ的ITD與方位角θ的關系 5.MATLAB仿真實驗

本實驗中采用的數據庫中采樣率為44.1KHz，時間分辨率為Ts=23?s的512點的離散序列——HRIR序列。ITD的單位是?s。參數具體是：-45°≤φ≤90°，0°≤θ≤360°。而HRIR序列是按不同φ確定的不同緯度面上，θ以人頭正前方為0°開始的，每5°變化一個方向取得hL??,??和hR??,??離散值。θ=0:5:355，這樣對于給定俯仰角φ的緯度面上就有72個方向的hL??,?i?和hR??,?i?離散值。為了方便記錄，將不同俯仰角?i下的雙耳時間差記為：ITD?i。（1）俯仰角φ=0°，方位角θ=0°；

程序如下： ITD0=[ ];Ts=23;load D:Signalshrtfselev0L0e000a.dat;hl0=L0e000a;load D:Signalshrtfselev0R0e000a.dat;hr0=R0e000a;c0=normxcorr2(hl0,hr0);[max_c0,imax]=max(abs(c0(:)));[yspeak,xspeak]=ind2sub(size(c0),imax(1));n0=[yspeak-size(hl0,1),xspeak-size(hl0,2)];t0=n0(1)*Ts;ITD0=[ITD0 t0];運行結果為： ITD0 =23；

（2）仰角φ=0°，方位角θ=5°； 運行結果為：ITD0 =[23 69]；（3）仰角φ=0°，方位角θ=10°;運行結果為：ITD0 =[23 69 92]；

這樣得到俯仰角φ=0°即水平面上的雙耳時間差 ITD01×72= [ 23 69 92 161

437 736 345 483 920 299 506 828 276

552 828 230

207 598 851 184

230 644 851 138

276 667 621

345 690 414

368713391

46 0-46-92-138-184-230-276-322-368-414-437-483-506-828-575-598-736-713-713-667-667-690-598-552-506-460-414-391-345-322-253-207-184-115-92-46 0 ]（4）同理可以變成計算出ITD15，ITD30，ITD45，ITD60，ITD75，ITD_15(φ=-15°)，ITD_30(φ=-30°)（5）得到的是7組離散的序列ITD?i，對其進行插值和平滑處理，基本可觀察出我們所需要的大體情況。程序如下： x=0:5:355;xi=0:0.01:355;yi0=interp1(x,ITD0,xi,'spline');yi15=interp1(x,ITD15,xi,'spline');yi30=interp1(x,ITD30,xi,'spline');yi45=interp1(x,ITD45,xi,'spline');yi60=interp1(x,ITD60,xi,'spline');yi75=interp1(x,ITD75,xi,'spline');yi_15=interp1(x,ITD_15,xi,'spline');yi_30=interp1(x,ITD_30,xi,'spline');yy0=smooth(yi0,0.1);yy15=smooth(yi15,0.1);yy30=smooth(yi30,0.1);yy45=smooth(yi45,0.1);yy60=smooth(yi60,0.1);yy75=smooth(yi75,0.1);yy_15=smooth(yi_15,0.1);yy_30=smooth(yi_30,0.1);figure(1);xi=0:0.01:355;plot(xi,yy0,'r',xi,yy15,'m',xi,yy30,'g',xi,yy45,'c',xi,yy60,'b',xi,yy75,'y',xi,yy_15,'k',xi,yy_30);grid on;axis([0 360-950 950]);xlabel('θ/(°)');ylabel('ITD/μs');仿真圖如下：

6.小結

通過計算左右耳接收到的信號之間的互相關，便可得到雙耳時間差ITD，通過觀察發現實驗中仿真圖與前文中所給的平均值還是有一定出入的。這也是必然會出現的結果，首先實驗所采用的是HRIR的時間分辨率低，還有考慮耳郭對高頻信號部分的影響。這在實驗中是沒有給予考慮的。

參考文獻

[1] Head-related transfer function database and its analyses

Xie BoSunt，ZHONG XiaoLi，RAO Dan&LIANG ZhiQiang Acoustics of China of 0641，China 2007；

[2] 頭相關傳輸函數與虛擬聽覺，謝菠蓀著，國防工業出版社2008； [3] MATLAB程序設計教程，李海濤，鄧櫻著，高等教育出版社2007；

第二篇：隨機信號處理教學文本

隨機信號處理教學大綱

課程名稱：隨機信號處理

學 時：45學時 開課學期：第六學期

適用專業：電子信息工程、電子科學與技術 課程類別：選修 課程性質：專業基礎課

先修課程：數字信號處理、概率論與數理統計、數字電路、計算機原理

教 材：《隨機信號處理》 張玲華，鄭寶玉著

清華大學出版社2003年9月第一版(一)本課程的地位、性質和任務

隨機信號是客觀世界中普遍存在的一類信號，對其特性的深入理解以及掌握相應的分析與處理方法，對電子信息工程專業的學生是非常重要的。本課程是電子信息工程、信息對抗技術專業的本科生掌握現代電子技術必備的一門學科基礎課。學習本課程的目的在于掌握信號統計分析與處理的理論和方法，通過學習，具備一定的隨機信號分析和處理的能力，為以后專業課學習打下基礎。(二)課程教學的基本要求：

通過該課程的學習，要求學生理解隨機信號的基本概念，掌握隨機信號的基本理論和分析處理方法，為學習“統計信號處理”或“信號檢測與估值”等后續課程以及將來的發展奠定堅實的基礎。

(三)課程主要內容及學時分配：

第1章 緒論(2學時)要求了解數字信號處理的基本概念，學科概貌，DSP的基本組成、特點等。主要包括下面幾部分內容：

1.1 數字信號處理的基本概念

1.2 數字信號處理的學科概貌（研究內容）1.3 數字信號處理系統的基本組成 1.4 數字信號處理的特點 1.5 本課程的特點

第1章 數字信號處理基礎（10學時）

要求掌握離散時間信號系統相關概念、數字濾波器的結構等內容。主要包括下面幾部分內容：

1.1 離散時間信號系統 1.2 數字濾波器的結構

2、《隨機過程理論及應用》，陸大鑫等，高等教育出版社，1987。

3、《Probability RandomVariable Radom process》帕布里斯（美）

4、《統計信號處理》 沈鳳麟，葉中付，錢玉美著 中國科技大學出版2001年3月(五)教學方法的原則性建議： 重點難點

1、隨機信號基本理論和概念的建立

2、基本隨機信號處理方法的掌握

3、現代譜估計理論和自適應信號處理技術

方法提示

授課、小結、習作討論、輔導與答疑相結合。

第三篇：隨機信號大作業

隨機信號大作業 第一章上機題：

設有隨機初相信號X(t)=5cos（t+），其中相位是在區間（0,2）上均勻分布的隨機變量。（1）試用Matlab編程產生其三個樣本函數。（2）產生t=0時的10000個樣本，并畫出直方圖 估計P(x)畫出圖形。

解：

（1）由Matlab產生的三個樣本函數如下圖所示：

程序源代碼：

clc clear m=unifrnd(0,2*pi,1,10);for k=1:3 t=1:0.1:10;X=5*cos(t+m(k));plot(t,X);hold on end xlabel(＇t＇);ylabel(＇X(t)＇);grid on;axis tight;（2）產生t=0時的10000個樣本，并畫出直方圖估計P(x)的概率密度并畫出圖形。

源程序代碼：

clear;clc;=2*pi*rand(10000,1);x=5*cos();figure(2),hist(x,20);hold on;第二章上機題：

利用Matlab程序設計一正弦型信號加高斯白噪聲的復合信號。

（1）分析復合信號的功率譜密度，幅度分布的特性；

（2）分析復合信號通過RC積分電路后的功率譜密度和相應的幅度分布特性；

（3）分析復合信號通過理想低通系統后的功率譜密度和相應的幅度分布特性。

解：

設正弦信號的頻率為10HZ，抽樣頻率為100HZ x=sin(2*pi*fc*t)正弦曲線圖：

程序塊代碼：

clear all;fs=100;fc=10;n=201;t=0:1/fs:2;x=sin(2*pi*fc*t);y=awgn(x,10);m=50;i=-0.49:1/fs:0.49;for j=1:m R(j)=sum(y(1:n-j-1).*y(j:199),2)/(n-j);Ry(49+j)=R(j);Ry(51-j)=R(j);end subplot(5,2,1);plot(t,x,＇r＇);title(＇正弦信號曲線＇);ylabel(＇x＇);xlabel(＇t/20pi＇);grid;（1）正弦信號加上高斯白噪聲產生復合信號y：

y=awgn(x,10)對復合信號進行傅里葉變換得到傅里葉變換：

Y(jw)=fft(y)復合信號的功率譜密度函數為：

G(w)=Y(jw).*conj(Y(jw)/length(Y(jw)))復合信號的曲線圖，頻譜圖和功率譜圖：

程序塊代碼：

plot(t,y,＇r＇);title(＇復合信號曲線＇);ylabel(＇y＇);xlabel(＇t/20pi＇);grid;程序塊代碼：

FY=fft(y);FY1=fftshift(FY);f=(0:200)*fs/n-fs/2;plot(f,abs(FY1),＇r＇);title(＇復合信號頻譜圖＇);ylabel(＇F(jw)＇);xlabel(＇w＇);grid;程序塊代碼：

P=FY1.*conj(FY1)/length(FY1);plot(f,P,＇r＇);title(＇復合信號功率譜密度圖＇);ylabel(＇G(w)＇);xlabel(＇w＇);grid;（2）正弦曲線的復合信號通過RC積分電路后得到信號為：

通過卷積計算可以得到y2 即：y2= conv2(y,b*pi^-b*t)y2的幅度分布特性可以通過傅里葉變換得到Y2(jw)=fft(y2)y2的功率譜密度G2(w)=Y2(jw).*conj(Y2(jw)/length(Y2(jw)))復合信號通過RC積分電路后的曲線頻譜圖和功率譜圖：

程序塊代碼：

b=10;y2=conv2(y,b*pi^-b*t);Fy2=fftshift(fft(y2));f=(0:400)*fs/n-fs/2;plot(f,abs(Fy2),＇r＇);title(＇復合信號通過RC系統后頻譜圖＇);ylabel(＇Fy2(jw)＇);xlabel(＇w＇);grid;程序代碼：

P2=Fy2.*conj(Fy2)/length(Fy2);plot(f,P2,＇r＇);title(＇復合信號通過ＲＣ系統后功率密度圖＇);ylabel(＇Gy2(w)＇);xlabel(＇w＇);grid;（3）復合信號 y通過理想濾波器電路后得到信號y3 通過卷積計算可以得到y3 即：y3=conv2(y,sin(10*t)/(pi*t))y3的幅度分布特性可以通過傅里葉變換得到 Y3(jw)=fft(y3)，y3的功率譜密度 G3(w)=Y3(jw).*conj(Y3(jw)/length(Y3(jw)))復合信號通過理想濾波器后的頻譜圖和功率密度圖：

程序塊代碼：

y3=conv2(y,sin(10*t)/(pi*t));Fy3=fftshift(fft(y3));f3=(0:200)*fs/n-fs/2;plot(f3,abs(Fy3),＇r＇);title(＇復合信號通過理想濾波器頻譜圖＇);ylabel(＇Fy3(jw)＇);xlabel(＇w＇);grid;程序塊代碼：

P3=Fy3.*conj(Fy3)/length(Fy3);plot(f3,P3,＇r＇);title(＇理想信號通過理想濾波器功率密度圖＇);ylabel(＇Gy3(w)＇);xlabel(＇w＇);grid;

第四篇：《隨機信號分析》實驗報告

《隨機信號分析》實驗報告

學號：

姓名：

2009年12月21日

實驗一：平穩隨機過程的數字特征

1、實驗目的“正文、小四宋體1.5倍行距”

2、實驗任務

3、實驗流程

4、實驗結果

5、實驗代碼

“代碼、五號宋體1倍行距”

1、實驗目的“正文、小四宋體1.5倍行距”

2、實驗任務

3、實驗流程

4、實驗結果

5、實驗代碼

“代碼、五號宋體1倍行距”

1、實驗目的“正文、小四宋體1.5倍行距”

2、實驗任務

3、實驗流程

4、實驗結果

5、實驗代碼

“代碼、五號宋體1倍行距”

1、實驗目的“正文、小四宋體1.5倍行距”

2、實驗任務

3、實驗流程

4、實驗結果

5、實驗代碼

“代碼、五號宋體1倍行距”

第五篇：隨機信號分析實驗報告

Ｈ ａ ａr r ｂ ｂi in n

Ｉ Ｉn ns st ti iｔ ｔ ｕ ｕt te e

o of f

T Te ec ch h ｎ ｎo o ｌ ｌo og gy y

實 驗 報 告 告

課程名稱:

隨機信號分析

院

系：

電子與信息工程學院

班

級：

姓

名:

學

號：

指導教師:

實驗時間:

實驗一、各種分布隨機數得產生

（一)實驗原理 1、、均勻分布隨機數得產生原理 產生偽隨機數得一種實用方法就是同余法,它利用同余運算遞推產生偽隨機數序列．最簡單得方法就是加同余法

為了保證產生得偽隨機數能在[0,１]內均勻分布,需要Ｍ為正整數，此外常數 c 與初值 y0 亦為正整數。加同余法雖然簡單,但產生得偽隨機數效果不好。另一種同余法為乘同余法,它需要兩次乘法才能產生一個［0，１]上均勻分布得隨機數

? ??

式中,ａ為正整數。用加法與乘法完成遞推運算得稱為混合同余法,即

??

?用混合同余法產生得偽隨機數具有較好得特性,一些程序庫中都有成熟得程序供選擇。

常用得計算語言如 Basic、Ｃ與 Matlab 都有產生均勻分布隨機數得函數可

以調用，只就是用各種編程語言對應得函數產生得均勻分布隨機數得范圍不同，有得函數可能還需要提供種子或初始化。

Matlaｂ提供得函數ｒａnd()可以產生一個在[0,1］區間分布得隨機數，rand（2,4）則可以產生一個在［０,1]區間分布得隨機數矩陣,矩陣為２行４列。Matlab 提供得另一個產生隨機數得函數就是 ranｄom(’unif’,a,b，N，M),unif 表示均勻分布,a與b就是均勻分布區間得上下界,N與M分別就是矩陣得行與列。

2、、隨機變量得仿真 根據隨機變量函數變換得原理,如果能將兩個分布之間得函數關系用顯式表達,那么就可以利用一種分布得隨機變量通過變換得到另一種分布得隨機變量。

若Ｘ就是分布函數為 F(x)得隨機變量,且分布函數 F(x)為嚴格單調升函數,令Ｙ=F(X）,則 Y 必為在[0，1]上均勻分布得隨機變量．反之,若 Y 就是在[0,1]上均勻分布得隨機變量,那么 即就是分布函數為 FX(x）得隨機變量。式中 F X??1()為F X()? 得反函數.這樣,欲求某個分布得隨機變量,先產生在[0,１］區間上得均勻分布隨機數,再經上式變換，便可求得所需分布得隨機數。

３、高斯分布隨機數得仿真 廣泛應用得有兩種產生高斯隨機數得方法，一種就是變換法,一種就是近似法.如果Ｘ1,X2 就是兩個互相獨立得均勻分布隨機數,那么下式給出得 Y１,Y2

便就是數學期望為 m,方差為得高斯分布隨機數，且互相獨立，這就就是變換法。

另外一種產生高斯隨機數得方法就是近似法.在學習中心極限定理時，曾提到 n 個在[0,1］區間上均勻分布得互相獨立隨機變量 Xi（i＝１，2…，ｎ）,當ｎ足夠大時,其與得分布接近高斯分布.當然,只要 n 不就是無窮大，這個高斯分布就是近似得。由于近似法避免了開方與三角函數運算，計算量大大降低。當精度要求不太高時，近似法還就是具有很大應用價值得.4、、各種分布隨機數得仿真 有了高斯隨機變量得仿真方法，就可以構成與高斯變量有關得其她分布隨機變量,如瑞利分布、指數分布與分布隨機變量。

(二）

實驗目得 在很多系統仿真得過程中,需要產生不同分布得隨機變量。利用計算機可以很方便地產生不同分布得隨機變量，各種分布得隨機變量得基礎就是均勻分布得隨機變量.有了均勻分布得隨機變量，就可以用函數變換等方法得到其她分布得隨機變量。

(三)實驗結果

附:源程序 suｂplot（2,2，1);

x=random（’ｕniｆ’，2，5，1，１024）； ploｔ(x)； tｉtｌｅ(’均勻分布隨機數’)subpｌｏt(２，2，2);G1=ｒandｏm(’Nｏrmaｌ＇,0,1,1，20００0)； ploｔ（G１)； ｔitle（’高斯分布隨機數’)ｓｕbplot(2，2,３）;G2=raｎdom(“Nｏrmal’,0,1,1,20０00);R=sqrt（G1、＊G1＋G2、＊G２）;plｏt（Ｒ);title(’瑞利分布隨機數’)suｂplot(2，２,4);G3=raｎdｏm(”Nｏrｍal’,0,１,1,2000０);G４=ｒandom(“Nｏrｍａl’，０,1,1,2０0０0）； Ｘ=G１、＊Ｇ1+G２、*G2+G３、*G3+Ｇ4、*G4； pｌoｔ（Ｘ);ｔｉtle（”x^2 分布隨機數＇）

實驗 二、隨機變量檢驗（一)實驗 原理 1、均值得計算 在實際計算時,如果平穩隨機序列滿足各態歷經性，則統計均值可用時間均值代替。這樣，在計算統計均值時,并不需要大量樣本函數得集合，只需對一個樣本函數求時間平均即可。甚至有時也不需要計算 N ? ? 時得極限,況且也不可能。通常得做法就是取一個有限得、計算系統能夠承受得 N 求時間均值與時間方差。根據強調計算速度或精度得不同,可選擇不同得算法。

設隨機數序列｛}，一種計算均值得方法就是直接計算下式中,xn 為隨機數序列中得第 n 個隨機數。

另一種方法就是利用遞推算法，第ｎ次迭代得均值也亦即前 n 個隨機數得均值為迭代結束后,便得到隨機數序列得均值 m m N ?

遞推算法得優點就是可以實時計算均值,這種方法常用在實時獲取數據得場合。

當數據量較大時,為防止計算誤差得積累,也可采用式中,m1 就是取一小部分隨機數計算得均值.２、方差得計算 計算方差也分為直接法與遞推法。仿照均值得做法

方差得遞推算法需要同時遞推均值與方差 m mnx mn n n n? ? ?? ? 1 11()

迭代結束后，得到隨機數序列得方差為

其它矩函數也可用類似得方法得到.３、統計隨機數得概率密度直方圖 假定被統計得序列得最大值與最小值分別為 a 與 b。將區間等分 M(M 應與被統計得序列得個數 N 相適應，否則統計效果不好。)份后得區間為,，…，,… ,。用，表示序列得值落在區間里得個數,統計序列得值在各個區間得個數，則就粗略地反映了隨機序列得概率密度得情況.用圖形方式顯示出來就就是隨機數得概率密度直方圖.（二）

實驗目得 隨機數產生之后,必須對它得統計特性做嚴格得檢驗。一般來講,統計特性得檢驗包括參數檢驗、均勻性檢驗與獨立性檢驗等.事實上,我們如果在二階矩范圍內討論隨機信號,那么參數檢驗只對產生得隨機數一、二階矩進行檢驗。我們可以把產生得隨機數序列作為一個隨機變量,也可以瞧成隨機過程中得一個樣本函數。不論就是隨機變量還就是隨機過程得樣本函數，都會遇到求其數字特征得情況，有時需要計算隨機變量得概率密度直方圖等.（三）

實驗結果

附：源程序 subplｏt(2,2，1）;x=raｎdom(“ｕnｉf”,２，５,1,１０2４）;ｈisｔ（x,2:0、２：5);ｔitle（’均勻分布隨機數直方圖’);s1=0 foｒ ｎ１＝1:１024

s１=x(n1)＋s１;end Mｅan1=s1/102４； t1=0 for n1＝１:10２4

ｔ1＝(x（n1）—Mean1)＾2＋t1;end Ｖaｒiance１＝ｔ1/１024;sｕbplot(2,2，2)； G1＝raｎdｏm(’Noｒｍal“,0，1,１，２0０00)； hist(G1，—4:0、2:4）； ｔitｌe(”高斯分布隨機數直方圖’）;ｓ2=０ foｒ n２=1：2０00０

ｓ2=Ｇ1（ｎ2）+s2； eｎｄ Mｅan2＝s2/20000； t2=0 fｏr n2＝1：20０00

t2=(G1(n2）-Mean２)^2+t2;enｄ Vaｒｉaｎｃe2=t2/20000； subｐｌot(２，2,3);G2=rａndom(’Nｏｒmal’,0，1，1,２0000）； R＝sqｒt（G１、*G1＋G2、*G２);hｉst（R,0:０、２:5）;tｉtｌｅ(“瑞利分布隨機數直方圖’)； s3=０ for n3=1：20000

ｓ3=R(ｎ3）+ｓ3;ｅnd Meａn３=s3/２00０0;t３＝0 ｆor n3=1：2０000

ｔ3=(R（n３)—Ｍｅan3）^2+t3;end Ｖarianｃｅ3=t3/２00００;subｐlot(2,２，4);G3=ｒandom（’Normaｌ”，0,1，1,200０0）;G４=raｎdom（“Ｎormal”,０,1,１，２0000）;X=Ｇ1、*G1＋G2、＊Ｇ2+G3、＊G3＋Ｇ4、＊G4； hist(X，0：０、5：30）;tiｔle(“x^2 分布隨機數直方圖’)ｓ4=0 for ｎ4=1:200０0

ｓ４=X(ｎ４）+s４;ｅnd Mean4=s4/2000０;t4＝0 for n４＝1:２00０0

t4=(X(ｎ4)-Meaｎ４)^2＋t4； end 實驗 三、中心極限定理得驗證(一）

實驗 原理 如果 n 個獨立隨機變量得分布就是相同得,并且具有有限得數學期望與方差，當 n 無窮大時，它們之與得分布趨近于高斯分布。這就就是中心極限定理中

得一個定理。

我們以均勻分布為例，來解釋這個定理。若 n 個隨機變量 Xi（ｉ=1，2，…,n)都為[0,1］區間上得均勻分布得隨機變量,且互相獨立，當 n 足夠大時，其與得分布接近高斯分布。

(二）

實驗目得 利用計算機產生均勻分布得隨機數。對相互獨立得均勻分布得隨機變量做與,可以很直觀瞧到均勻分布得隨機變量得與,隨著做與次數得增加分布情況得變化，通過實驗對中心極限定理得進行驗證。

((三）

實驗結果

分析:隨ｎ取值得增大，均勻分布隨機序列求與得圖形越發接近于高斯分布。

附：源程序 X0=random(＇uniｆ”,0,1，１,１0２4);X1=randｏm(’unｉｆ’,0,1，1，１024)；

Ｘ2=ｒandoｍ（＇ｕniｆ“,0,1,1，１024);X3=rａndom(＇unif＇，0,１,１,1024);

Ｘ4=ｒaｎｄｏｍ(”uｎif＇，0，１，1，１０24);

X5=raｎdom(’ｕnif’，0，1,１,１024）；

Ｘ６=ｒaｎdom（’ｕnif“,0,1，1,1024）;X7=raｎｄom（’ｕｎiｆ’,0,1，１,1024);

X８=raｎdom(＇unif”，0，1，1,１０24);

X9=ｒaｎdom(’ｕnｉf’,0,１,１,1０２４）； Ｇ=rａｎｄom(“ｎｏｒmal”,０，1，1，1024）;

Ｙ1=X０+X1+X2＋X3＋Ｘ4；

Y2=Ｘ0+X1＋X2+Ｘ3+X４+Ｘ5+X6+X7+Ｘ8+X9;

suｂploｔ(2,2,1);ｈｉsｔ(X０,０:0、２：2)；

titlｅ(“均勻分布隨機數直方圖’）

sｕbplot(2，２,2）;hist（Y１,0：0、２：6);

ｔitlｅ(’五個均勻分布之與隨機數直方圖”)sｕｂplot(2,２,3）;ｈｉｓｔ(Y2，0：0、2：8);

title（’十個均勻分布之與隨機數直方圖“)subpｌoｔ(2，２，4);ｈiｓt（G,－4：0、2:4）;title(”高斯分布隨機數直方圖“）

實驗 四、中心極限定理得驗證(一）

實驗 原理 在實際應用中，我們可以把產生得隨機數序列瞧成隨機過程中得一個樣本函數。如果平穩隨機序列滿足各態歷經性,則統計自相關序列可用時間自相關序列

代替。當數據得樣本數有限時,也只能用有限個數據來估計時間自相關序列,統計自相關序列得估值。若各態歷經序列Ｘ(n）得一個樣本有 N 個數據,由于實序列自相關序列就是對稱得，自相關函數得估值為

(二)實驗目得 在隨機信號理論中,自相關函數就是非常重要得概念。在實際系統仿真中也會經常計算自相關函數．通過本試驗學生可以親自動手計算自相關函數,加深對概念得理解，并增強實際動手能力．（三))實驗結果

分析：分別生成均值為 0 與１，方差為 1 得高斯隨機數,由圖形可以明顯瞧出兩者自相關函數得差異。

附:源程序 N=２56;xn=ranｄom(’nｏｒm＇，０,１,1，N);Rx=ｘcoｒr（ｘn,＇biａsed”);ｍ=-N+1：N－1;ｓubplｏｔ(2,1，１）;ｐlot(m,Rx);tｉtle(“均值為０，方差為１得高斯分布得自相關函數＇）； axis（［—N N—1 —0、５ 1、５］)； N=256;xn＝randｏｍ（’noｒm’,1,1，1，Ｎ);Ｘk=ｆft(xn,2*N)； Rx＝ifｆｔ（(aｂs（Xk)、＾2)／N）； m=-N:Ｎ—１;suｂｐlot(２，1，2)； plot（m，fｆtsｈift（Rx));ｔitｌe（’均值為 1,方差為 1 得高斯分布得自相關函數’);ａxｉs（[-N N—1-0、5 １、5]);實驗五、功率譜密度(一)實驗 原理 一般把平穩隨機序列得功率譜定義為自相關序列得傅里葉變換。如果自相關序列就是周期序列, X(ｎ）得功率譜與自相關序列得關系為

? 與實平穩過程一樣,實平穩序列得功率譜也就是非負偶函數,即

可以證明，功率譜還可表示為

當 X(ｎ）為各態歷經序列時,可去掉上式中得統計均值計算,將隨機序列 X(n)用它得一個樣本序列 x（ｎ)代替。在實際應用中，由于一個樣本序列得可用數據個數 N 有限,功率譜密度也只能就是估計

式中，X(x（n)得傅里葉變換．這就是比較簡單得一種估計方法,這種功率譜密度得估計方法稱為周期圖方法。如果直接利用數據樣本做離散傅里葉變換,可得到 X(FＦT 算法實現,所以得到了廣泛得應用。

(二）實驗目得 在隨機信號理論中，功率譜密度與自相關函數一樣都就是非常重要得概念.在實際系統仿真中也會經常計算。通過本試驗學生可以親自動手,加深對概念得理解，并增強實際動手能力。

(三）實驗結果

附:源程序 N＝２5６;ｘ1＝random(”normal’,0,1，1,N);Sx1＝abs（ｆft(ｘ1）、＾2）／N;subplot(2,１,1);plｏｔ（１0*loｇ10（Sx1));ｔｉtle(“均值為０,方差為 1 得高斯分布得功率譜密度＇)； xlａｂel(’ｆ/Ｈz’)ylabel(”Sｘ１/dB’）

x２=rａndom(’ｎormａl“，１,1,1,N）； Ｓｘ2=abs（fｆt（x2)、＾２)/N;ｓuｂplot（2，1,2);plot（１０*log1０(Sx2));title（”均值為 1，方差為 1 得高斯分布得功率譜密度’）;xlabel(’ｆ/Hz＇）

ｙlabｅｌ（“Ｓx2/dB＇)實驗 六、隨機信號經過 線性系統前后信號仿真

（一))實驗原理

需要先仿真一個指定系統,再根據需要仿真輸入得隨機信號，然后使這個隨機信號通過指定得系統．通過對實際系統建模，計算機可以對很多系統進行仿真。在信號處理中,一般將線性系統分解為一個全通放大器(或衰減器)與一個特定頻率響應得濾波器。由于全通放大器可以用一個常數代替,因此線性系統得仿真往往只需設計一個數字濾波器。濾波器設計可采用 MＡTLAB 提供得函數，也可

利用相應得方法自行設計。MＡＴLAＢ提供了多個設計濾波器得函數,可以很方便地設計低通、帶通、高通、多帶通、帶阻濾波器。

((二）實驗 目得

系統仿真就是信號仿真處理得一個重要部分，通過該實驗要求學生掌握系統仿真得基本概念,并學會系統得仿真方法。

((三))實驗 結果

1、低通濾波器

2、帶通濾波器

3、高通濾波器 4、多帶通濾波器

5、帶阻濾波器

附：源程序 １、X（n)

N＝2000；fｓ=４00;Nn=ｒandom（”normａl＇,0,1,1，N)； t=（0：Ｎ—１)/fs;fi=ｒaｎｄｏm(’unｉf’,0,１,１,２)＊2＊pｉ;xn＝sin(2*ｐi*５０*ｔ+ｆｉ(1))+Nn;Ｒx=ｘcorｒ(xn，“biaｓed’)； ｍ=—N+1:N－1;Sx=abs（ｆft（xn）、^2）/Ｎ； f=(—N／2:N/2-１）＊fs/N;sｕbｐlｏt（211),ｐｌot(m，Rｘ)； xlabｅl(’m’）

ｙｌabel（”Rx(m）’)titlｅ（’xn 得自相關函數“);subploｔ(212）,plｏt（f,ffｔshift（10＊lｏg10（Sx(1:N）)）);xlabeｌ（’f/Hz”)ｙlabｅl(“Sx/dB”)tｉtle(’xn 得功率譜密度’）;2、低通濾波器 ｈ=ｆir1（１00,0、４);H=ｆft(h,2＊N);HW=aｂs(Ｈ)、＾２;Rx=xcoｒr（xn，’ｂｉasｅｄ＇);Sx=abs（ｆfｔｓhift（fft(xn，2*N））、＾2)/（2＊Ｎ)； Ｓy＝Sx、*ＨＷ;Ｒy＝fftｓｈift（ifft（Sｙ)）;

f=（-Ｎ：N—1)＊ｆs/(2＊N)； m=(—N:Ｎ－1);subｐｌot(311);pｌot（(-N:Ｎ—1)/N，fｆtshiｆt(aｂs(HW(1:2＊N))）);tiｔlｅ(＇低通濾波器“）;ｓubplｏt(31２)，ｐｌoｔ（ｍ,Ry);xlabｅl(”m“)ylａbeｌ（”Ry（ｍ)＇)title(’xn 經低通濾波器得自相關函數’)； ｓubｐloｔ(31３),plot（ｆ，ffｔshift(10*log10(Ｓy（１:2*N)）))； ａxis(［—200 20０ —20 20]);xlabeｌ(“f/Hz’)ｙlabｅl(＇Sｙ/ｄＢ”)titｌe(＇xn 經低通濾波器得功率譜密度“)； 3、帶通濾波器 ｈ=fｉr1(10０，［０、１ ０、5]);H=ｆft(h，2*Ｎ);ＨW＝abs(H)、^2； Ｒx=ｘcorｒ(xn，”biａsed“）； Sｘ=aｂs（fｆtshiｆt(fｆt(ｘn,2*N））、＾2)／（2*N）； Sy=Ｓx、*ＨW； Rｙ=fｆtsｈift(ｉfft（Sｙ)）； f=(-Ｎ:N-1）*fs/(2＊Ｎ);m=（-N:N—1）;subｐlot(311);ｐlot((—N:N-1）/N,fftｓｈift（abs（HW（1:2＊Ｎ)))）； ｔitle(’帶通濾波器”）； ｓubｐlｏｔ（312)，ploｔ（ｍ,Ry);xlabeｌ(’m“)yｌabel(’Ｒy(m)’)title(”xn 經帶通通濾波器得自相關函數“)； ｓubpｌoｔ(31３)，plｏt(f,fftshift(10＊lｏg10（Sy(１：2＊N))））； ａxiｓ（［—２０0 20０ －２0 20]);xｌabｅl（’f/Hz”)yｌabｅl(“Ｓｙ/dＢ’)tｉtle（’xn 經帶通濾波器得功率譜密度’）;4、高通濾波器 h=ｆｉr1（１0０,0、6,’high’）； Ｈ=fft(h,2＊N）； ＨW=abｓ（H)、^２;Ｒx=xcｏｒr(xｎ,”bｉａｓed“）;Ｓx=abs(fftｓｈift（fｆt(xn,2＊N))、^２)/（2*N)； Sｙ=Ｓx、*HW;Ｒy=ffｔshiｆt(iｆfｔ（Sｙ));f=(-Ｎ:N-1）*ｆs/（2*N);m=(—Ｎ:N—1）；

sｕｂｐlｏt(31１)；plot(（-N：Ｎ—1）/N,ｆftshiｆt(aｂs(HＷ（1:２*N)）));tiｔｌe(＇高通濾波器”）;subpｌot(312），ploｔ(m,Rｙ);xlabｅl(“ｍ’)ylaｂeｌ(’Rｙ(m)”)tiｔlｅ(＇xn 經高通通濾波器得自相關函數’);subplｏt（3１3）,ｐｌot(ｆ,fftshｉft(10*ｌog10(Sy（1:2＊Ｎ)）））;axis（［－２０0 2０0 —2０ 2０]）； xlａbeｌ（“ｆ/Ｈz’)ｙｌａbｅl（”Sy/ｄB“)tｉtle(＇xn 經高通濾波器得功率譜密度＇);5、多帶通濾波器 h=fｉr１(100,[０、１,0、3，0、5,0、7］)； H＝ｆfｔ(h，2*N)； HW＝ａｂｓ（Ｈ)、^2； Ｒx=xcorr（xn,＇bｉasｅd’);Sx=abs（fｆtshift(fft(xn,2*N））、^2)/（２＊Ｎ）； Sy＝Ｓx、*HW;Ｒy=ffｔshiｆt（iｆft(Sy));f＝（—Ｎ:N—1)＊fｓ/(2*N);m＝（—N：Ｎ-1);sｕbｐloｔ(３１1）;ploｔ(（—N:N—１）/Ｎ,fftｓｈift（ａbs(HW（1:２＊N）)))； tiｔle(’多帶通濾波器’)； ｓｕbplot(312),pｌｏt(m,Ry）； ｘlabel（＇m’)ｙｌabｅl（”Ry(m）“）

ｔitle（”xn 經多帶通通濾波器得自相關函數“);subpｌot（31３),plot(ｆ,fｆｔshｉｆt（1０＊ｌog10(Sｙ（1：2＊Ｎ)）));axis（［－2０0 2０0 —20 20]);xlaｂeｌ（’ｆ／Hz”）

ｙlabｅl(“Sy/dB’）

ｔitle(’xn 經多帶通濾波器得功率譜密度”）； 6、帶阻濾波器 h=fiｒ１(100,[０、1,0、4],’sｔoｐ’);H=fft（h，2＊N);HＷ=abｓ(H)、^2;Rx=ｘｃｏrr（xｎ,’ｂiaseｄ“）;Ｓx＝ａbｓ（ffｔshift(fft（ｘｎ，２*N）)、＾2)／(2*N);Sy=Sx、＊ＨW； Ｒy＝fftshift（ｉfft(Sy））;f=(—N:N-1）*ｆｓ/（2*N);m=(-N：N—1)； ｓｕｂｐlot(311）；pｌoｔ(（—N：N-1）/Ｎ,fftshift(ａｂs（ＨＷ(1:２*Ｎ))));

ｔｉtｌｅ(”帶阻濾波器“）； ｓubpｌot（31２),ｐlot（m，Ｒy)； xlａbeｌ(’m’）

ｙlａbel（”Ｒy(m)’)ｔitle（’ｘn 經帶阻濾波器得自相關函數＇)； sｕbplｏｔ(3１3)，ｐlot(f,fftshift（10＊log10(Sy（１:2*N)）））;aｘis([－200 20０-20 20］);xlaｂｅl（＇f／Ｈz“)ylabｅl（”Sy/ｄB“)titlｅ（”xn 經帶阻濾波器得功率譜密度"）;