matlab求均值,方差

实验报告

随机信号的数字特征分析

一、 实验目的

1.了解随机信号自身的特性，包括均值（数学期望）、方差、均方值等； 2. 掌握随机信号的分析方法；

二、实验原理

1.均值测量方法

ˆx表示集合平均值或数学期望值。基于随机过程的各态历经性，均值m最常用的方法是取N个样本数据并简单地进行平均，即

1ˆxmN2.均方误差的测量方法

随机序列的均方误差定义为：

Xd[i] i0N1其中，样本信号的采样数据记为Xd[i]X(iT,)，Ts为采样间隔。

1N2E(X)limxi(n)

NNi123.方差测量方法

如果信号的均值是已知的，则其方差估计设计为

1ˆN2X(Xd[i]mx)2 i0N1它是无偏的与渐进一致的。

三、实验内容

利用MATLAB中的伪随机序列产生函数randn()产生多段1000点的序列，编制一个程序，计算随机信号的数字特征，包括均值、方差、均方值、最后把计算结果平均，绘制数字特征图形。

源程序如下： clear all;

文档

实用

clc;

%产生50个1000以内点的伪随机序列 x=randn(50,1000);

%计算随机产生的50个点序列的均值，方差，均方 average=zeros(1,50); variance=zeros(1,50); square=zeros(1,50);

%计算均值 for i=1:50

for j=1:1000

average(i)=average(i)+x(i,j); end

average(i)=average(i)/1000;

end

%计算方差 for i=1:50

for j=1:1000

variance(i)=variance(i)+(x(i,j)-average(i)).^2; end

variance(i)=variance(i)/1000; end

%计算均方值 for i=1:50

for j=1:1000

square(i)=square(i)+x(i,j).^2; end

square(i)=square(i)/1000; end

EX=sum(average)/50; DX=sum(variance)/50; RMS=sum(square)/50;

plot(average);

title('50个随机序列的均值'); figure;

文档

实用

plot(variance);

title('50个随机序列的方差'); figure;

plot(square);

title('50个随机序列的均方值');

四、实验结果及分析

文档

实用

文档

实用

由上结果可知：

将图中的计算结果平均后，得到的结果为：产生的50个点的随机序列均值的平均值为：EX=0.0090197；产生的50个点的随机序列方差的平均值为DX=1.0078；产生的50个点的随机序列均方值的平均值为RMS=1.0087。

由上面所得到的图形可以看出50个点的伪随机序列的均值都在0附近，方差以及均方差都在1附近，将这些均值平均后得出的均值也是在0值附近，方差在1附近，与统计的结果相符合。

文档

实用

实验二 数字相关和数字卷积程序

一、实验目的

熟悉数字相关和数字卷积运算。

二、实验原理

1.线性以及循环相关的原理 1.1 线性相关的原理

假定x1(n)是列长为N的有限长序列，x2(n)是列长为M的有限长序列，两者的线性相关的结果为：

y(n)1.2 循环相关的原理

x1(m)x2(mn) m假定x1(n)是列长为N的有限长序列，x2(n)是列长为M的有限长序列，两者循环相关的结果为：

y(n)x1(m)x2((mn))NRN(n)

m0N12.线性以及循环卷积的原理 2.1 线性卷积的原理

假定x1(n)是列长为N的有限长序列，x2(n)是列长为M的有限长序列，两者的线性卷积的结果为：

y(n)x1(n)*x2(n)mx(m)x(nm)

122.2 循环卷积的原理

循环卷积的矩阵表示形式如下所示：其中x和H是两个输入的序列，y是循环卷积得到的实验结果。

yHx

其中，y[y(0),y(1),...,y(N1)]T，x[x(0),x(1),...,x(N1)]T

文档

实用

h(N1)h(1)h(0)h(1)h(0)h(2) Hh(N1)h(N2)h(0)

三、实验内容

编写函数实现两个随机序列的线性、循环相关和线性、循环卷积的程序:

源程序如下：

两个序列线性相关的函数： clear all clc

x=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; end

for i=nx+1:n x(i) = 0; end

for i=1:n

for j=1:n

H(i,j) = h(mod(i+j-2,n)+1); end end

y = H * x';

subplot(3, 1, 1);stem(x);title('随机序列1'); subplot(3, 1, 2);stem(h);title('随机序列2'); subplot(3, 1, 3);stem(y);title('线性相关结果');

两个序列循环相关的函数： clear all clc

文档

实用

x=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; end end

if (nxfor i=nx+1:n x(i) = 0;

end end

for i=1:n

for j=1:n

H(i,j) = h(mod(i+j-2,n)+1); end

end

y = H * x';

subplot(3, 1, 1);stem(x);title('随机序列1'); subplot(3, 1, 2);stem(h);title('随机序列2'); subplot(3, 1, 3);stem(y);title('循环相关结果');

两个序列线性卷积的函数： clear all clc

x=ones(1,8); h=ones(1,10); nx = length(x); nh = length(h); n = nx + nh - 1; for i = nx+1:n x(i) = 0; end

for i=nh+1:n

文档

实用

h(i) = 0;

end

for i=1:n

for j=1:n

H(i,j) = h(mod(i+n-j,n)+1); end end

y = H * x';

subplot(3, 1, 1);stem(x);title('随机序列1'); subplot(3, 1, 2);stem(h);title('随机序列2'); subplot(3, 1, 3);stem(y);title('线性卷积结果');

两个序列循环卷积的函数： clear all clc

x=ones(1,8); h=ones(1,10); n=15;

nx = length(x); nh = length(h);

if (nfprintf('输入圆周卷积的点数不正确'); break end

for k=nh+1:n h(k) = 0; end

for k=nx+1:n x(k) = 0; end

for k=1:n

for l=1:n

H(k,l)=h(mod(k+n-l,n)+1); end end

y = H*x';

subplot(3, 1, 1);stem(x);title('随机序列1'); subplot(3, 1, 2);stem(h);title('随机序列2');

文档

实用

subplot(3, 1, 3);stem(y);title('循环卷积结果');

四、实验结果及分析

1.线性相关实现的程序及结果

y={ 8 8 8 6 5 4 3 2 1 1 2 3 4 5 6 7}

2.循环相关实现的程序及结果

y={ 8 8 8 8 8 8 8 8 8 8}

文档

实用

3.线性卷积实现的程序及结果

y={ 1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1}

文档

实用

4.循环卷积实现的程序及结果 当n=15时

y={ 3 3 3 4 5 6 7 8 8 8 7 6 5 4 3}

当n=17时

y={ 1 2 3 4 5 6 7 8 8 8 7 6 5 4 3 2 1}

文档

实用

由上图可知：

15点循环卷积结果与线性卷积的结果是不一致的，但是17点循环卷积结果与线性卷积的结果是一致的。

文档

实用

实验三 维纳－霍夫方程的求解

一、实验目的

学习使用Matlab实现W-H程序的编写。

二、实验原理

一个线性系统，如果它的单位样本响应为h(n)，当输入一个随机信号x(n) ：

x(n)s(n)(n)

其中s(n)表示信号，表示噪声，则输出y(n)为

y(n)h(m)x(nm)

m我们希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因

ˆ(n) ˆ(n)表示，即 y(n)s此称y(n)为s(n)的估计值，用s维纳滤波的标准方程

ˆ分别表示信号的真值与估计值，而用e(n)表示它们之如果我们以s与s间的误差

ˆ(n) e(n)s(n)s目标：均方误差 E{e(n)2}min (MMSE准则)

ˆ(n)h(m)x(nm) y(n)sm0ˆ(n)hixi，上式可看成输出等于现在和过去各输入的加权之和 s其中

i1im1 或 mi1

hih(i1)h(m)xix(ni1)x(nm)2最小的h，ˆ现在的问题是需要求得使Ess为此，将这式对hi求i偏导，并令其结果等于0，得

2Eshixixj0于是Es(n)h(m)x(nm)x(nk)i1 j10 k0 0ptm0Ee(n)x0j这样就得到维纳滤波的标准方程

文档

实用

sx(k)hopt(m)xx(km), k0

m0FIR维纳滤波器

设h(n)是一个因果序列可以用有限长(长度为N)的序列去逼近它 ，有上述得到W-H方程的矩阵形式为：

i1 h1x1x1h2x1x2i2 h1x2x1h2x2x2hNx2xNx2s

iN h1xNx1h2xNx2hNxNxNxNs即：xxhxs，其中xx自相关矩阵称，xs为x与s的互相关矩阵

hNx1xNx1sx1sx1x1 x1x2 x1xNh1  x2x2x2xNh2x2sx2x1 hxsxx   xsX X XhNxN2xNN2N1N这样得到W-H方程的解为：hh0ptxxxs

1三、实验内容

编写函数解W-H方程，寻找最优的滤波器，并检验该程序的准确性。

源程序如下： clear all; clc;

%输入信号 A = 1; f = 1000;

fs = 10^5; t = (0:999); Mlag = 100;

x = A*cos(2*pi*f*t/fs);

%给正弦波信号加入信噪比为20dB的高斯白噪声 xn = awgn(x,5); figure;

subplot(2,2,1)

文档

实用

plot(t,xn)

title('输入信号图像') %计算输入信号自相关函数

Rxn=xcorr(xn,Mlag,'biased'); subplot(2,2,2)

plot((-Mlag:Mlag),Rxn)

title('输入信号自相关函数') %维纳滤波 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));

for i = 2:N

c = Rxn(101+i)*ones(1,N+1-i);

Rxx = Rxx+diag(c,i-1)+diag(c,-i+1); end Rxx;

h = zeros(N,1);

h = inv(Rxx)*rxnx; yn = filter(h,1,xn); subplot(2,2,3) plot(yn);

title('经过维纳滤波器后信号信号'); Ryn=xcorr(yn,Mlag,'biased'); subplot(2,2,4);

plot((-Mlag:Mlag),Ryn);

title('经过维纳滤波器后信号自相关函数');

四、实验结果

文档

实用

从图中可以看出，滤波后得到的正弦信号仍然有一定的误差，但是输入信号的自相关函数在0点出有明显的噪声成分，通过维纳滤波以后得到的信号的自相关函数在0点处已经的噪声给消除了很多。

文档

实用

实验四 Yule－Walker方程的求解

一、实验目的

学习使用Matlab实现Y-W程序的编写。

二、实验原理

1.AR模型的 Yule－Walker方程

AR模型,又称为自回归模型,是一个全极点的模型,可用如下差分方程来表示：

x(n)akx(nk)(n) （1）

k1p其中(n)是均值为零、方差为的白噪声序列, p是AR模型的阶数, a( k), k= 1, 2, , p 是p阶AR模型的参数。AR 模型系统H ( z ) 的转移函数为：

H(z)11akzkk1p （2）

从而得到AR模型的功率谱估计的计算公式：

Pxx()2A(e)j221akejkk1p2 （3）

由上式可以看出, 要利用AR模型进行功率谱估计, 必须得到模型参数和白噪声序列的方差。将（1）式变形有：

pxx(m)akxx(mk)0 , m0k1 （4） p(m)a(k)2, m0xxkxxk1式（4）的矩阵形式为：

2xx(0) xx(1) xx(2) xx(p)1(1) (0) (1) ((p-1))a 0xxxxxxxx1 （5）    (p) (p1) (p2) (0)xxap 0xxxxxx式( 4)和( 5) 是AR模型的AR Yule-Walker方程。

文档

实用

2.burg算法求解方法

Burg算法是使序列x(n)的前后向预测误差功率之和：

22E2e(n)b(n1)Kb(n1)Keppp1p1p1p1(n)0 （6） 最小。

利用Burg法求解AR模型参数的步骤：

第一步：由初始条件e0(n)x(n)，b0(n)x(n)根据公式（7）求出反射系数

K1：

2ep1(n)bp1(n1)N1KpNn1p （7）

22e(n)bp1p1(n1)np1第二步：根据序列x（n）自相关函数x(0)Nn0N1x(n)，求出阶次m=1

2时的AR模型参数a(1,1)=k1与前后向预测误差功率之和。

第三步：由式(8)求出前向预测误差e1(n)与后向预测误差b1(n), 然后由式(7)估计出反射系数k2；

ep(n)ep1(n)Kpbp1(n1)bp(n)bp1(n1)Kpep1(n)数a(2,1)和2；

（8）

第四步：由（9）Levinsion递推关系, 求出阶次m=2时的AR模型参

Kpappapiap1,iappap1,pi222p(1app)p, 10xx(0)2 （9）

第五步: 重复上述过程, 直到阶次m=p, 这样就求出了所有阶次的AR模型参数。

Burg算法的递推过程是建立在数据序列基础上, 避开了序列的自相关函数的估计, 所以与自相关法相比, 具有较好的频率分辨率。

三、实验内容

已知观测信号，编写函数解Y-W方程，寻找参数系数，并检验该程序的准确性和掌握用法。

文档

实用

源程序如下：

clear all;

x = cos(0:0.1:50); N = 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谱估计 p=8;

a=zeros(p,p);

%初始化前向和后向误差以及γ e=x; b=x;

sigma=0; for i=1:N

sigma=sigma+x(i).^2; end

sigma=sigma/N; for n=1:p sum1=0; sum2=0;

for j=n+1:N

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

%计算并输出功率谱 for m=1:N sum=0; for n=1:p

sum = a(p,n)*exp(-sqrt(-1)*2*pi*n*m/N) + sum; end

POW2(m)=sigma/(abs(1 + sum).^2);

end

subplot(1,2,2),plot(1/N : 2*pi/N/2/pi : 0.5, POW2(1:N/2)); title('AR模型谱估计法求取功率谱'); xlabel('f/Hz');

四．实验结果及分析

实验中输入的信号为余弦信号，理想情况下其功率谱是在余弦信号频率上的一个冲击函数。从实验的结果图可以看出，用AR模型估计的功率

文档

实用

谱同用周期图法估计的功率谱一样，说明了用该方法所计算的功率谱的准确性。

文档

实用

实验五 自适应噪声抵消算法的软件设计与实现

一、实验目的

学习使用MATLAB编写LMS自适应滤波器，以及如何在生物医学信号中进行应用。

二、实验原理

1、自适应干扰抵消的原理

图1 自适应干扰抵消原理图

图1所示的是自适应干扰抵消器的基本结构。期望信号d(n)是信号与噪声之和，即d(n)= x(n) +N(n)，自适应处理器的输入是与N(n)相关的另一个噪声N’(n)。当x(n) 和N(n)不相关时，自适应处理器将调整自己的参数，使y(n)成为N(n)的最佳估计N(n)。这样，e(n)将逼近信号x(n)，且其均方差E[e2(n)]为最小。噪声N(n)就得到了一定程度的抵消。 2、LMS自适应滤波算法

图2 单输入自适应线性组合器

LMS算法使用的准则是使滤波器的期望输出值和实际输出值之间的均方误最小化的准则，即使用均方误差来做性能指标。

文档

实用

自适应滤波的结果如图2所示。各符号的意义是：x(n)输入信号，y(n)为滤波器的输出，d(n)为y(n)想要趋近的理想信号，d(n) 是已知的，e(n)为误差信号。滤波器均方误差可表示为：(n)E[e2(n)]

设自适应滤波器的输入矢量为：X(n)[x(n),x(n1),......,x(nN1)] 加权矢量（即滤波器参数矢量）为：W(n)[w1(n),w2(n),...,xN(n)]T

滤波器的输出为：y(n)WT(n)X(n)XT(n)W(n)

误差信号为期望输出d(n)与滤波器实际输出之间的误差，即 e(n)d(n)y(n)d(n)WT(n)X(n) （1）

LMS算法是取单个误差样本的平方e2(n)的梯度作为均方误差梯度的估计，由式（1）可得梯度矢量的估计为：

e2(n)e(n)(n)2e(n)2e(n)x(n)

ww由此得到一个新的权矢量递推公式，即LMS算法递推公式为：

w(n1)w(n)u((n))w(n)2ue(n)x(n)

三、实验内容

已知观测信号，编写LMS滤波器，对该信号进行滤波处理，检验该程序的准确性和掌握MATLAB自带函数的用法。

源程序如下：

clear all;

d = sin(0:0.1:50); M = length(d);

% 设定滤波器的长度为15 N=15;

noise = 0.2*sin(50*(0:0.1:50));%噪声信号 x = d + noise； u=0.01;

%初始化滤波器的参数为全0的矩阵 w=zeros(1,N); for n=N:M y(n)=0; for i=1:N

文档

实用

y(n)=y(n)+x(n-i+1)*w(i);

end

e(n)=d(n)-y(n);

for i=1:N %进行滤波器参数的调整 w(i)=w(i)+2*u*e(n)*x(n-i+1); end end

subplot(2,1,1); plot(y);

title('滤波处理后得到的信号') subplot(2,1,2); plot(x);

title('观测到的信号');

四、实验结果及分析

实验中输入的原始信号为正弦信号，噪声为在正弦信号的基础上加入一个50倍频率的噪声信号，观测到的信号是这两个信号的叠加，利用LMS算法进行滤波处理，从从实验结果中可以看出，在LMS的加权参数稳定之后的输出结果比较理想。

文档

实用

文档