matlab求均值-方差

1、實(shí)驗(yàn)報(bào)告隨機(jī)信號(hào)的數(shù)字特征分析、方差、均方值等;一、實(shí)驗(yàn)?zāi)康? .了解隨機(jī)信號(hào)自身的特性，包括均值(數(shù)學(xué)期望)2 .掌握隨機(jī)信號(hào)的分析方法；二、實(shí)驗(yàn)原理1 .均值測量方法均值亂表示集合平均值或數(shù)學(xué)期望值?；陔S機(jī)過程的各態(tài)歷經(jīng)性,最常用的方法是取N個(gè)樣本數(shù)據(jù)并簡單地進(jìn)行平均，即r?x1N1Ni0Xdi其中，樣本信號(hào)的采樣數(shù)據(jù)記為2 .均方誤差的測量方法XdiX(iT,XTs為采樣間隔。隨機(jī)序列的均方誤差定義為:E(X2)1NOlimx2(n)NNii')3 .方差測量方法如果信號(hào)的均值是已知的，則其方差估計(jì)設(shè)計(jì)為1N1?XV(Xdimx)2Ni0它是無偏的與漸進(jìn)一致的。三、實(shí)驗(yàn)內(nèi)容利用

2、MATLABH勺偽隨機(jī)序列產(chǎn)生函數(shù)randn()產(chǎn)生多段1000點(diǎn)的序歹【,編制一個(gè)程序，計(jì)算隨機(jī)信號(hào)的數(shù)字特征，包括均值、方差、均方值、最后把計(jì)算結(jié)果平均，繪制數(shù)字特征圖形。源程序如下:clearall;clc;%產(chǎn)生50個(gè)1000以內(nèi)點(diǎn)的偽隨機(jī)序列x=randn(50,1000);%計(jì)算隨機(jī)產(chǎn)生的50個(gè)點(diǎn)序列的均值，方差，均方average=zeros(1,50);variance=zeros(1,50);square=zeros(1,50);%計(jì)算均值fori=1:50forj=1:1000average(i)=average(i)+x(i,j);endaverage(i)=averag

3、e(i)/1000;end%計(jì)算方差fori=1:50forj=1:1000variance(i)=variance(i)+(x(i,j)-average(i).A2;endvariance(i)=variance(i)/1000;end%計(jì)算均方值fori=1:50forj=1:1000square(i)=square(i)+x(i,j).A2;endsquare(i)=square(i)/1000;endEX=sum(average)/50;DX=sum(variance)/50;RMS=sum(square)/50;plot(average);title('50個(gè)隨機(jī)序列的均值&

4、#39;);figure;plot(variance);title('50個(gè)隨機(jī)序列的方差');figure;plot(square);title('50個(gè)隨機(jī)序列的均方值');四、實(shí)驗(yàn)結(jié)果及分析砌個(gè)隨機(jī)序列的均方值由上結(jié)果可知：將圖中的計(jì)算結(jié)果平均后，得到的結(jié)果為：產(chǎn)生的50個(gè)點(diǎn)的隨機(jī)序列均值的平均值為：EX=0.0090197;產(chǎn)生的50個(gè)點(diǎn)的隨機(jī)序列方差的平均值為DX=1.0078;產(chǎn)生的50個(gè)點(diǎn)的隨機(jī)序列均方值的平均值為RMS=1.0087由上面所得到的圖形可以看出50個(gè)點(diǎn)的偽隨機(jī)序列的均值都在0附近,方差以及均方差都在1附近，將這些均值平均后得出的均值

5、也是在0值附近，方差在1附近，與統(tǒng)計(jì)的結(jié)果相符合。實(shí)驗(yàn)二 數(shù)字相關(guān)和數(shù)字卷積程序、實(shí)驗(yàn)?zāi)康氖煜?shù)字相關(guān)和數(shù)字卷積運(yùn)算。二、實(shí)驗(yàn)原理1. 線性以及循環(huán)相關(guān)的原理1.1 線性相關(guān)的原理假定x1(n)是列長為N的有限長序列，x2(n)是列長為M的有限長序歹兩者的線性相關(guān)的結(jié)果為：y(n)x1(m)x2(mn)m1.2 循環(huán)相關(guān)的原理假定x1(n)是列長為N的有限長序列，x2(n)是列長為M的有限長序歹兩者循環(huán)相關(guān)的結(jié)果為：N1y(n)x1(m)x2(mn)NRN(n)m02. 線性以及循環(huán)卷積的原理2.1 線性卷積的原理假定x1(n)是列長為N的有限長序列，x2(n)是列長為M的有限長序歹1,兩者的

6、線性卷積的結(jié)果為：y(n)x1(n)*x2(n)x1(m)x2(nm)m2.2 循環(huán)卷積的原理循環(huán)卷積的矩陣表示形式如下所示：其中x和H是兩個(gè)輸入的序列，是循環(huán)卷積得到的實(shí)驗(yàn)結(jié)果。x x(0),x(1),.,x(N1)TyHx其中，yy(0),y(1),.,y(N1)Th(0)h(N1)h(1)h(1)h(0)h(2)Hh(N1)h(N2)h(0)三、實(shí)驗(yàn)內(nèi)容的程編寫函數(shù)實(shí)現(xiàn)兩個(gè)隨機(jī)序列的線性、循環(huán)相關(guān)和線性、循環(huán)卷積序:源程序如下：兩個(gè)序列線性相關(guān)的函數(shù)：clearallclcx=ones(1,8);h=ones(1,10);nx=length(x);nh=length(h);n=nx+nh

7、-1;fori=nh+1:nh(i)=0;endfori=nx+1:nx(i)=0;endfori=1:nforj=1:nH(i,j)=h(mod(i+j-2,n)+1);endend隨機(jī)序列1');隨機(jī)序列2');線性相關(guān)結(jié)果');y=H*x'subplot(3,1,1);stem(x);title('subplot(3,1,2);stem(h);title('subplot(3,1,3);stem(y);title('兩個(gè)序列循環(huán)相關(guān)的函數(shù)：clearallclcx=ones(1,8);h=ones(1,10);nx=length(x

8、);nh=length(h);n=nx;if(nx>nh)fori=nh+1:nh(i)=0;endendif(nx<nh)n=nh;fori=nx+1:nx(i)=0;endendfori=1:nforj=1:nH(i,j)=h(mod(i+j-2,n)+1);endend隨機(jī)序列1');隨機(jī)序列2');循環(huán)相關(guān)結(jié)果');y=H*x'subplot(3,1,1);stem(x);title('subplot(3,1,2);stem(h);title('subplot(3,1,3);stem(y);title('兩個(gè)序列線性卷

9、積的函數(shù)：clearallclcx=ones(1,8);h=ones(1,10);nx=length(x);nh=length(h);n=nx+nh-1;fori=nx+1:nx(i)=0;endfori=nh+1:nh(i)=0;endfori=1:nforj=1:nH(i,j)=h(mod(i+n-j,n)+1);endendy=H*x'隨機(jī)序列1');隨機(jī)序列2');線性卷積結(jié)果');subplot(3,1,1);stem(x);title('subplot(3,1,2);stem(h);title('subplot(3,1,3);stem

10、(y);title('兩個(gè)序列循環(huán)卷積的函數(shù)：clearallclcx=ones(1,8);h=ones(1,10);n=15;nx=length(x);nh=length(h);if(n<nx|n<nh)fprintf('輸入圓周卷積的點(diǎn)數(shù)不正確');breakendfork=nh+1:nh(k)=0;endfork=nx+1:nx(k)=0;endfork=1:nforl=1:nH(k,l)=h(mod(k+n-l,n)+1);endendy=H*x'隨機(jī)序列1');隨機(jī)序列2');subplot(3,1,1);stem(x);t

11、itle('subplot(3,1,2);stem(h);title(subplot(3,1,3);stem(y);title(循環(huán)卷積結(jié)果');四、實(shí)驗(yàn)結(jié)果及分析1.線性相關(guān)實(shí)現(xiàn)的程序及結(jié)果y=88865432112345672 .循環(huán)相關(guān)實(shí)現(xiàn)的程序及結(jié)果y=88888888883 .線性卷積實(shí)現(xiàn)的程序及結(jié)果y=123456788876543214.循環(huán)卷積實(shí)現(xiàn)的程序及結(jié)果當(dāng)n=15時(shí)y=333456788876543當(dāng)n=17時(shí)y=12345678887654由上圖可知：15點(diǎn)循環(huán)卷積結(jié)果與線性卷積的結(jié)果是不一致的，但是17點(diǎn)循環(huán)卷積結(jié)果與線性卷積的結(jié)果是一致的。實(shí)驗(yàn)三維納

12、-霍夫方程的求解一、實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)使用Matlab實(shí)現(xiàn)W-H程序的編寫。二、實(shí)驗(yàn)原理一個(gè)線性系統(tǒng)，如果它的單位樣本響應(yīng)為h(n),當(dāng)輸入一個(gè)隨機(jī)信號(hào)x(n):x(n)s(n)(n)其中s(n)表示信號(hào)，表示噪聲，則輸出y(n)為y(n)h(m)x(nm)m我們希望x(n)通過線性系統(tǒng)h(n)后得到的y(n)盡量接近于s(n),因此稱y(n)為s(n)的估計(jì)值，用S(n)表示，即y(n)S(n)維納濾波的標(biāo)準(zhǔn)方程如果我們以片？分別表示信號(hào)的真值與估計(jì)值，而用e(n)表示它們之問的誤差e(n)s(n)s(n)目標(biāo)：均方誤差E|e(n)2min(MMSE準(zhǔn)則)y(n)s(n)h(m)x(nm)m0上式

13、可看成輸出等于現(xiàn)在和過去各輸入的加權(quán)之和Mn)hx，其中i1im1或mi1hih(i1)h(m)為x(ni1)x(nm)現(xiàn)在的問題是需要求得使Es§2最小的(h),為此，將這式對(h)求偏導(dǎo)，并令其結(jié)果等于0,得2Eshxixj0i1j1于正Es(n)h0Dt(m)x(nm)x(nk)0k0pm0Ee(n)xj0這樣就得到維納濾波的標(biāo)準(zhǔn)方程sx(k)hopt(m)xx(km),k0m0FIR維納濾波器設(shè)h(n)是一個(gè)因果序列可以用有限長(長度為N)的序列去逼近它，有上述得到W-H方程的矩陣形式為：i1hlx1x1h2乂伐LhNxNx1si2hlx2x1h2x2x2LhNx?xnx?s

14、MMiNhlxn*h2xnx2LhNxnxnxNs即：xxh其中xs，八1xx自相關(guān)矩陣稱，xs為x三ss的互相關(guān)矩h1xsxx1取2Lxxnhh2x?sx2x1x2x2Lx2xNMxsMxxMMMhNxNs2NX1xNX2LxNX這樣得到W-H方程的解為：hh0pt1xxxs三、實(shí)驗(yàn)內(nèi)容編寫函數(shù)解W-H方程，尋找最優(yōu)的濾波器，并檢驗(yàn)該程序的準(zhǔn)確性。源程序如下：clearall;clc;%俞入信號(hào)A=1;f=1000;fs=10A5;t=(0:999);Mlag=100;x=A*cos(2*pi*f*t/fs);%合正弦波信號(hào)加入信噪比為20dB的高斯白噪聲xn=awgn(x,5);figur

15、e;subplot(2,2,1)plot(t,xn)title('輸入信號(hào)圖像')%計(jì)算輸入信號(hào)自相關(guān)函數(shù)Rxn=xcorr(xn,Mlag,'biased');subplot(2,2,2)plot(-Mlag:Mlag),Rxn)title('輸入信號(hào)自相關(guān)函數(shù)')%維納濾波N=100;Rxnx=xcorr(xn,x,Mlag,'biased');rxnx=zeros(N,1);rxnx(:)=Rxnx(101:101+N-1);Rxx=zeros(N,N);Rxx=diag(Rxn(101)*ones(1,N);fori=2:

16、Nc=Rxn(101+i)*ones(1,N+1-i);Rxx=Rxx+diag(c,i-1)+diag(c,-i+1);endRxx;h=zeros(N,1);h=inv(Rxx)*rxnx;yn=filter(h,1,xn);subplot(2,2,3)plot(yn);title('經(jīng)過維納濾波器后信號(hào)信號(hào)');Ryn=xcorr(yn,Mlag,'biased');subplot(2,2,4);plot(-Mlag:Mlag),Ryn);');title('經(jīng)過維納濾波器后信號(hào)自相關(guān)函數(shù)四、實(shí)驗(yàn)結(jié)果輸入信號(hào)自相關(guān)函數(shù)經(jīng)過維納濾波器后信號(hào)信

17、號(hào)輕過維納濾波器后信號(hào)自相關(guān)函數(shù)從圖中可以看出，濾波后得到的正弦信號(hào)仍然有一定的誤差，但是輸入信號(hào)的自相關(guān)函數(shù)在0點(diǎn)出有明顯的噪聲成分，通過維納濾波以后得到的信號(hào)的自相關(guān)函數(shù)在0點(diǎn)處已經(jīng)的噪聲給消除了很多。實(shí)驗(yàn)四Yule-Walker方程的求解、實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)使用Matlab實(shí)現(xiàn)Y-W程序的編寫二、實(shí)驗(yàn)原理1.AR模型的Yule-Walker方程AR模型，又稱為自回歸模型，是一個(gè)全極點(diǎn)的模型，可用如下差分方程來表示：px(n)akx(nk)(n)(1)k1其中(n)是均值為零、方差為的白噪聲序列，p是AR模型的階數(shù)a(k),k=1,2,p是p階AR模型的參數(shù)。AR模型系統(tǒng)H(z)的轉(zhuǎn)移函數(shù)為：H

18、(z)(2)1k1akzk1從而得到AR模型的功率譜估計(jì)的計(jì)算公式:巳()由上式可以看出 數(shù)和白噪聲序列的方差xx(m)xx(m)2A(ej )d pj k1ake Jk 1,要利用AR模型進(jìn)行功率譜估計(jì)(3)必須得到模型參式(4) xx(0) xx(1)Mxx( p)將(pakk 1pakk 11)式變形有:xx(m k)(4)的矩陣形式為:xx( 1) xx(0)Mxx(p 1)xx (k)xxxx，( M2)1)xx ( p)xx( (p-1)M1a1M(5)xx(p 2)xx(0)ap 0式(4)和(5)是AR模型的ARYule-Walker方程。BurgE最小。利用Burg法求解AR

19、模型參數(shù)的步驟:第一步：由初始條件eo(n) x(n), b0(n)x(n)根據(jù)公式(7)求出反射系數(shù)K1:KpN 12ep1(n)bp(n 1)n pN1e2 1(n) bp 1(n 1) n p第二步：根據(jù)序列x (n)自相關(guān)函數(shù)x(0)N 1N x(n)N n 02,求出階次m=12.burg算法求解方法算法是使序列x(n)的前后向預(yù)測誤差功率之和：(6)2ep1(n)bp1(n1)Kpb；1(n1)Kpep1(n)0時(shí)的AR模型參數(shù)a(1,1)=k1與前后向預(yù)測誤差功率之和。第三步：由式(8)求出前向預(yù)測誤差e1(n)與后向預(yù)測誤差U(n),然后由式(7)估計(jì)出反射系數(shù)k2;ep(n)

20、bp(n)ep 1(n) Kpbp 1(n 1) bp 1(n 1) Kpep 1(n)(8)第四步：由( 數(shù) a(2,1)和 2；Kpapi2 p9) Levinsion遞推關(guān)系,求出階次 m=2時(shí)的AR模型參appap 1,i(1appappap 1,p2) 2 ) p 1,xx(0)(9)第五步：重復(fù)上述過程，直到階次m=p,這樣就求出了所有階次的AR,避開了序列的自相關(guān)模型參數(shù)。Burg算法的遞推過程是建立在數(shù)據(jù)序列基礎(chǔ)上函數(shù)的估計(jì)，所以與自相關(guān)法相比，具有較好的頻率分辨率。三、實(shí)驗(yàn)內(nèi)容已知觀測信號(hào)，編寫函數(shù)解Y-W方程，尋找參數(shù)系數(shù)，并檢驗(yàn)該程序的準(zhǔn)確性和掌握用法。源程序如下：cle

21、arall;x=cos(0:0.1:50);N=length(x);%周期圖法FFTX=fft(x);POW=(abs(FFTX)A2)/N;subplot(1,2,1),plot(1/N:2*pi/N/2/pi:0.5,POW(1:N/2);title('周期圖求取功率譜');xlabel('f/Hz');%Burg法白AR譜估計(jì)p=8;a=zeros(p,p);%初始化前向和后向誤差以及Te=x;b=x;sigma=0;fori=1:Nsigma=sigma+x(i).A2;endsigma=sigma/N;forn=1:psum1=0;sum2=0;for

22、j=n+1:Nsum1=sum1+2*e(j)*b(j-1);sum2=sum2+e(j).A2+b(j-1).A2;enda(n,n)=-sum1/sum2;sigma=sigma*(1-abs(a(n,n).A2);ifn>=2fori=1:n-1a(n,i)=a(n-1,i)+a(n,n)*a(n-1,n-i);endendforj=n+1:Nc(j)=e(j)+a(n,n)*b(j-1);d(j)=b(j-1)+a(n,n)*e(j);ende=c;b=d;end%計(jì)算并輸出功率譜form=1:Nsum=0;forn=1:psum=a(p,n)*exp(-sqrt(-1)*2*p

23、i*n*m/N)+sum;endPOW2(m)=sigma/(abs(1+sum).A2);endsubplot(1,2,2),plot(1/N:2*pi/N/2/pi:0.5,POW2(1:N/2);title('AR模型譜估計(jì)法求取功率譜');xlabel('f/Hz');四.實(shí)驗(yàn)結(jié)果及分析140周期圖求取功率譜120100804030D.20萬0 80100050000.20,4060 8M模型譜估計(jì)法求取功率諾 1500f/Hz0實(shí)驗(yàn)中輸入的信號(hào)為余弦信號(hào)，理想情況下其功率譜是在余弦信號(hào)頻率上的一個(gè)沖擊函數(shù)。從實(shí)驗(yàn)的結(jié)果圖可以看出，用AR模型估計(jì)的功率譜

24、同用周期圖法估計(jì)的功率譜一樣，說明了用該方法所計(jì)算的功率譜的準(zhǔn)確性。實(shí)驗(yàn)五自適應(yīng)噪聲抵消算法的軟件設(shè)計(jì)與實(shí)現(xiàn)一、實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)使用MATLAB編寫LMS自適應(yīng)濾波器，以及如何在生物醫(yī)學(xué)信號(hào)中進(jìn)行應(yīng)用。二、實(shí)驗(yàn)原理1、自適應(yīng)干擾抵消的原理Nf (n)自適應(yīng)處理器，圖1自適應(yīng)干擾抵消原理圖圖1所示的是自適應(yīng)干擾抵消器的基本結(jié)構(gòu)。期望信號(hào)d(n)是信號(hào)與噪聲之和，即d(n)=x(n)+N(n),自適應(yīng)處理器的輸入是與N(n)相關(guān)的另一個(gè)噪聲N'(n)。當(dāng)x(n)和N(n)不相關(guān)時(shí)，自適應(yīng)處理器將調(diào)整自己的參數(shù)，使y(n)成為N(n)的最佳估計(jì)N(n)。這樣，e(n)將逼近信號(hào)x(n),且其均方差Ee2(n)為最小。噪聲N(n)就得到了一定程度的抵消。2、LMS自適應(yīng)濾波算法圖2單輸入自適應(yīng)線性組合器LMS算法使用的準(zhǔn)則是使濾波器的期望輸出值和實(shí)際輸出值之間的均方誤最小化的準(zhǔn)則，即使用均方誤差來做性能指標(biāo)。自適應(yīng)濾波的結(jié)果如圖2所示。各符號(hào)的意義是：x(n)輸入信號(hào)，y(n)為濾波器的輸出，d(n)為y(n)想要趨近的理想信號(hào)，d(n)是已知的，e(n)為誤差信號(hào)。濾波器均方誤差可表示為：(n)Ee2(n)設(shè)自適應(yīng)濾波器的輸入矢量為：X(n)