matlab求均值,方差

1、隨機(jī)過(guò)程及在生物醫(yī)學(xué)工程中的應(yīng)用實(shí)驗(yàn)報(bào)告實(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.均值測(cè)量方法均值表示集合平均值或數(shù)學(xué)期望值?；陔S機(jī)過(guò)程的各態(tài)歷經(jīng)性，最常用的方法是取N個(gè)樣本數(shù)據(jù)并簡(jiǎn)單地進(jìn)行平均，即其中，樣本信號(hào)的采樣數(shù)據(jù)記為，為采樣間隔。2.均方誤差的測(cè)量方法隨機(jī)序列的均方誤差定義為：3.方差測(cè)量方法如果信號(hào)的均值是已知的，則其方差估計(jì)設(shè)計(jì)為它是無(wú)偏的與漸進(jìn)一致的。三、實(shí)驗(yàn)內(nèi)容利用MATLAB中的偽隨機(jī)序列產(chǎn)生函數(shù)randn()產(chǎn)生多段1000點(diǎn)的序列，編制一個(gè)程序，計(jì)算隨機(jī)信

2、號(hào)的數(shù)字特征，包括均值、方差、均方值、最后把計(jì)算結(jié)果平均，繪制數(shù)字特征圖形。源程序如下：clear all;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ì)算均值for i=1:50 for j=1:1000 average(i)=average(i)+x(i,j); end average(i)=average(i)/1000;end%計(jì)算方差for i=1:50 for j=1:1000 v

3、ariance(i)=variance(i)+(x(i,j)-average(i).2; end variance(i)=variance(i)/1000;end%計(jì)算均方值for i=1:50 for j=1:1000 square(i)=square(i)+x(i,j).2; end square(i)=square(i)/1000;endEX=sum(average)/50;DX=sum(variance)/50;RMS=sum(square)/50;plot(average);title('50個(gè)隨機(jī)序列的均值');figure;plot(variance);title

4、('50個(gè)隨機(jī)序列的方差');figure;plot(square);title('50個(gè)隨機(jī)序列的均方值');四、實(shí)驗(yàn)結(jié)果及分析由上結(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附近，將這些均值平均后得出的均值也是在0值附近，方差在1附近，與統(tǒng)計(jì)的結(jié)果相符合。實(shí)驗(yàn)二 數(shù)字相關(guān)和數(shù)字卷積程序一、

5、實(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)是列長(zhǎng)為N的有限長(zhǎng)序列，x2(n)是列長(zhǎng)為M的有限長(zhǎng)序列，兩者的線性相關(guān)的結(jié)果為：1.2 循環(huán)相關(guān)的原理假定x1(n)是列長(zhǎng)為N的有限長(zhǎng)序列，x2(n)是列長(zhǎng)為M的有限長(zhǎng)序列，兩者循環(huán)相關(guān)的結(jié)果為：2.線性以及循環(huán)卷積的原理2.1 線性卷積的原理假定x1(n)是列長(zhǎng)為N的有限長(zhǎng)序列，x2(n)是列長(zhǎng)為M的有限長(zhǎng)序列，兩者的線性卷積的結(jié)果為：2.2 循環(huán)卷積的原理循環(huán)卷積的矩陣表示形式如下所示：其中x和H是兩個(gè)輸入的序列，y是循環(huán)卷積得到的實(shí)驗(yàn)結(jié)果。其中，三、實(shí)驗(yàn)內(nèi)容編寫函數(shù)實(shí)現(xiàn)兩個(gè)隨機(jī)

6、序列的線性、循環(huán)相關(guān)和線性、循環(huán)卷積的程序:源程序如下：兩個(gè)序列線性相關(guān)的函數(shù)：clear allclcx=ones(1,8);h=ones(1,10);nx = length(x);nh = length(h);n = nx + nh - 1;for i = nh+1:n h(i) = 0;endfor i=nx+1:n x(i) = 0;endfor i=1:n for j=1:n H(i,j) = h(mod(i+j-2,n)+1); endendy = H * x' subplot(3, 1, 1);stem(x);title('隨機(jī)序列1');subplot(

7、3, 1, 2);stem(h);title('隨機(jī)序列2');subplot(3, 1, 3);stem(y);title('線性相關(guān)結(jié)果');兩個(gè)序列循環(huán)相關(guān)的函數(shù)：clear allclcx=ones(1,8);h=ones(1,10);nx = length(x);nh = length(h);n = nx;if (nx>nh) for i = nh+1:n h(i) = 0; endendif (nx<nh) n = nh; for i=nx+1:n x(i) = 0; endendfor i=1:n for j=1:n H(i,j) =

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

9、= 0;endfor i=nh+1:n h(i) = 0;endfor i=1:n for j=1:n H(i,j) = h(mod(i+n-j,n)+1); endendy = H * x'subplot(3, 1, 1);stem(x);title('隨機(jī)序列1');subplot(3, 1, 2);stem(h);title('隨機(jī)序列2');subplot(3, 1, 3);stem(y);title('線性卷積結(jié)果');兩個(gè)序列循環(huán)卷積的函數(shù)：clear allclcx=ones(1,8);h=ones(1,10);n=15;n

10、x = length(x);nh = length(h);if (n<nx|n<nh) fprintf('輸入圓周卷積的點(diǎn)數(shù)不正確'); breakendfor k=nh+1:n h(k) = 0;endfor k=nx+1:n x(k) = 0;endfor k=1:n for l=1:n H(k,l)=h(mod(k+n-l,n)+1); endendy = H*x'subplot(3, 1, 1);stem(x);title('隨機(jī)序列1');subplot(3, 1, 2);stem(h);title('隨機(jī)序列2')

11、;subplot(3, 1, 3);stem(y);title('循環(huán)卷積結(jié)果');四、實(shí)驗(yàn)結(jié)果及分析1.線性相關(guān)實(shí)現(xiàn)的程序及結(jié)果y= 8 8 8 6 5 4 3 2 1 1 2 3 4 5 6 72.循環(huán)相關(guān)實(shí)現(xiàn)的程序及結(jié)果y= 8 8 8 8 8 8 8 8 8 83.線性卷積實(shí)現(xiàn)的程序及結(jié)果y= 1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 14.循環(huán)卷積實(shí)現(xiàn)的程序及結(jié)果當(dāng)n=15時(shí)y= 3 3 3 4 5 6 7 8 8 8 7 6 5 4 3當(dāng)n=17時(shí)y= 1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1由上圖可知：15點(diǎn)循環(huán)卷積

12、結(jié)果與線性卷積的結(jié)果是不一致的，但是17點(diǎn)循環(huán)卷積結(jié)果與線性卷積的結(jié)果是一致的。實(shí)驗(yàn)三 維納霍夫方程的求解一、實(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) ：其中s(n)表示信號(hào)，表示噪聲，則輸出y(n)為 我們希望x(n)通過(guò)線性系統(tǒng)h(n)后得到的y(n)盡量接近于s(n)，因此稱y(n)為s(n)的估計(jì)值，用表示，即 維納濾波的標(biāo)準(zhǔn)方程 如果我們以分別表示信號(hào)的真值與估計(jì)值，而用e(n)表示它們之間的誤差目標(biāo)：均方誤差 min (MMSE準(zhǔn)則) 上式可看成輸出等于現(xiàn)在和過(guò)去各輸入的加權(quán)之和 ，其中現(xiàn)

13、在的問(wèn)題是需要求得使最小的，為此，將這式對(duì)求偏導(dǎo)，并令其結(jié)果等于0，得于是這樣就得到維納濾波的標(biāo)準(zhǔn)方程 FIR維納濾波器設(shè)h(n)是一個(gè)因果序列可以用有限長(zhǎng)(長(zhǎng)度為N)的序列去逼近它 ，有上述得到W-H方程的矩陣形式為：即：，其中自相關(guān)矩陣稱，為x與s的互相關(guān)矩陣 這樣得到W-H方程的解為：三、實(shí)驗(yàn)內(nèi)容編寫函數(shù)解W-H方程，尋找最優(yōu)的濾波器，并檢驗(yàn)該程序的準(zhǔn)確性。源程序如下：clear all;clc;%輸入信號(hào)A = 1;f = 1000;fs = 105; t = (0:999); Mlag = 100; x = A*cos(2*pi*f*t/fs); %給正弦波信號(hào)加入信噪比為20dB的

14、高斯白噪聲xn = awgn(x,5); figure;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

15、(N,N); Rxx = diag(Rxn(101)*ones(1,N);for i = 2:N c = 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)過(guò)維納濾波器后信號(hào)信號(hào)');Ryn=xcorr(yn,Mlag,'biased'); subplot(2,2,4);plot(-Mlag:Mlag),Ry

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

17、 p 是p階AR模型的參數(shù)。AR 模型系統(tǒng)H ( z ) 的轉(zhuǎn)移函數(shù)為： （2）從而得到AR模型的功率譜估計(jì)的計(jì)算公式： （3）由上式可以看出, 要利用AR模型進(jìn)行功率譜估計(jì), 必須得到模型參數(shù)和白噪聲序列的方差。將（1）式變形有： （4）式（4）的矩陣形式為： （5）式( 4)和( 5) 是AR模型的AR Yule-Walker方程。2.burg算法求解方法Burg算法是使序列x(n)的前后向預(yù)測(cè)誤差功率之和： （6）最小。利用Burg法求解AR模型參數(shù)的步驟：第一步：由初始條件，根據(jù)公式（7）求出反射系數(shù)： （7）第二步：根據(jù)序列x（n）自相關(guān)函數(shù)，求出階次m=1時(shí)的AR模型參數(shù)a(1,1

18、)=k1與前后向預(yù)測(cè)誤差功率之和。第三步：由式(8)求出前向預(yù)測(cè)誤差與后向預(yù)測(cè)誤差, 然后由式(7)估計(jì)出反射系數(shù)k2； （8）第四步：由（9）Levinsion遞推關(guān)系, 求出階次m=2時(shí)的AR模型參數(shù)a(2,1)和； （9）第五步: 重復(fù)上述過(guò)程, 直到階次m=p, 這樣就求出了所有階次的AR模型參數(shù)。Burg算法的遞推過(guò)程是建立在數(shù)據(jù)序列基礎(chǔ)上, 避開了序列的自相關(guān)函數(shù)的估計(jì), 所以與自相關(guān)法相比, 具有較好的頻率分辨率。三、實(shí)驗(yàn)內(nèi)容已知觀測(cè)信號(hào)，編寫函數(shù)解Y-W方程，尋找參數(shù)系數(shù)，并檢驗(yàn)該程序的準(zhǔn)確性和掌握用法。源程序如下：clear all;x = cos(0:0.1:50);N =

19、 length(x);%周期圖法FFTX = fft(x);POW = (abs(FFTX).2)/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);%初始化前向和后向誤差以及e=x;b=x;sigma=0;for i=1:N sigma=sigma+x(i).2;endsigma=sigma/N;for n=1:p sum1=0; sum2=0; for j=n+1:N

20、sum1=sum1+2*e(j)*b(j-1); sum2=sum2+e(j).2+b(j-1).2; end a(n,n)=-sum1/sum2; sigma=sigma*(1-abs(a(n,n).2); if n>=2 for i=1:n-1 a(n,i)=a(n-1,i)+a(n,n)*a(n-1,n-i); end end for j=n+1:N c(j)=e(j)+a(n,n)*b(j-1); d(j)=b(j-1)+a(n,n)*e(j); end e=c; b=d;end%計(jì)算并輸出功率譜for m=1:N sum=0; for n=1:p sum = a(p,n)*ex

21、p(-sqrt(-1)*2*pi*n*m/N) + sum; end POW2(m)=sigma/(abs(1 + sum).2);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é)果及分析實(shí)驗(yàn)中輸入的信號(hào)為余弦信號(hào)，理想情況下其功率譜是在余弦信號(hào)頻率上的一個(gè)沖擊函數(shù)。從實(shí)驗(yàn)的結(jié)果圖可以看出，用AR模型估計(jì)的功率譜同用周期圖法估計(jì)的功率譜一樣，說(shuō)明了用該方法所計(jì)算的功率譜的準(zhǔn)確性。實(shí)驗(yàn)五 自適應(yīng)噪聲抵消算法的軟

22、件設(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)干擾抵消的原理圖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ì)。這樣，e(n)將逼近信號(hào)x(n)，且其均方差Ee2(n)為最小。噪聲N(n)就得到了一定程度的抵消。2、LMS自適應(yīng)濾波算法圖2 單輸入自適應(yīng)線性組合器LMS算法使用的準(zhǔn)則是使濾波器的期望輸出值和實(shí)際輸出值之間的均方誤最小化的準(zhǔn)則，即使用均方誤差來(lái)做性能指標(biāo)。自適應(yīng)濾波的結(jié)果如圖2所示。各符號(hào)的意義是：x(n)輸入信號(hào)，y(n)為濾波器的輸出，d(n)為y(n)想要趨