matlab求均值方差Word格式.docx
《matlab求均值方差Word格式.docx》由会员分享，可在线阅读，更多相关《matlab求均值方差Word格式.docx（30页珍藏版）》请在冰豆网上搜索。
10０0
variaｎce（ｉ）＝vaｒiance（i）+（ｘ（i,j）—ａvｅrａｇｅ（i））。
^２；
end
variａnｃe（ｉ）=variancｅ（i）/1０0０；
ｅｎd
％计算均方值
for ｉ＝1:
50
ｆｏrj＝１:
１0０0
sqｕａrｅ（i）=sqｕａｒｅ（i）+x（i,ｊ）。
^2;
end
squａｒｅ（ｉ）＝sｑuａrｅ（i）/1000；
end
EX＝sum（avｅrａge）/5０；
DＸ=sum（variaｎcｅ）/５0;
RMS=sum（squaｒe）／5０；
plot（average）;
titｌe（’50个随机序列的均值'
）;
fｉｇuｒe;
ｐlｏt（varianｃe）；
title（'
50个随机序列的方差'
figurｅ;
plｏt（sqｕare）;
titｌe（’50个随机序列的均方值'
四、实验结果及分析
由上结果可知：
将图中的计算结果平均后，得到的结果为：
产生的50个点的随机序列均值的平均值为:
ＥX=０。
00９0197;
产生的５0个点的随机序列方差的平均值为ＤX=１。
0078;
产生的50个点的随机序列均方值的平均值为RMS=１.0０87．
由上面所得到的图形可以看出50个点的伪随机序列的均值都在０附近，方差以及均方差都在1附近，将这些均值平均后得出的均值也是在0值附近,方差在1附近，与统计的结果相符合。
ﻬ
实验二数字相关和数字卷积程序
熟悉数字相关和数字卷积运算。
１。
线性以及循环相关的原理
1线性相关的原理
假定x1（n）是列长为Ｎ的有限长序列,x2（ｎ）是列长为Ｍ的有限长序列，两者的线性相关的结果为:
2循环相关的原理
假定x1（n）是列长为Ｎ的有限长序列，x２（n）是列长为M的有限长序列,两者循环相关的结果为:
线性以及循环卷积的原理
2.1 线性卷积的原理
假定x1（n）是列长为N的有限长序列，x2（n）是列长为M的有限长序列,两者的线性卷积的结果为：
2．2循环卷积的原理
循环卷积的矩阵表示形式如下所示:
其中x和H是两个输入的序列,y是循环卷积得到的实验结果。
其中,
编写函数实现两个随机序列的线性、循环相关和线性、循环卷积的程序:
源程序如下：
两个序列线性相关的函数：
clｅaｒall
cｌc
x=ｏｎes（1,8）;
h＝ones（1,１0）;
nx ＝lenｇｔh（x）；
nｈ=ｌenｇｔｈ（h）;
n=nｘ + nｈ-1；
for i=nh＋１:
n
h（i）＝0;
fｏｒi＝nx＋1：
x（i）＝0；
ｅnd
for ｉ=1:
ｆoｒj=1:
H（i，j）=h（mod（i+ｊ—２,n）+1）;
ｅnd
y=H＊ x’;
subpｌoｔ（3，1，1）;
stem（x）；
title（’随机序列１'
）；
subpｌot（3,１,2）;
sｔem（ｈ）;
tiｔｌe（’随机序列2’）;
subｐloｔ（3,１,３）;
steｍ（y）；
ｔiｔle（'
线性相关结果'
两个序列循环相关的函数：
ｃlearａｌl
clc
ｘ=oneｓ（１,8）;
h＝onｅs（1,10）；
nｘ=lｅngth（x）;
ｎh=lｅｎgtｈ（ｈ）；
n =nｘ；
if（ｎx〉nｈ）
foｒi = nh＋1:
ｈ（i）= 0;
if（ｎx<
nｈ）
n =nh；
for i＝nｘ+１:
ｘ（i） = ０；
enｄ
foｒｉ＝１：
foｒj=1:
Ｈ（ｉ,j）=h（mｏｄ（i+ｊ—2，ｎ）＋１）;
eｎd
y=H* x'
；
subpｌot（3,1， 1）;
ｓtem（ｘ）;
tｉtlｅ（’随机序列１'
subpｌot（3,１，2）;
sｔｅm（h）；
title（’随机序列2’）；
subpｌｏｔ（3,1,３）；
ｓｔem（y）;
循环相关结果＇）;
两个序列线性卷积的函数:
cleａraｌｌ
x=ｏnes（１,8）;
ｈ＝oｎｅｓ（1，１0）;
nx =leｎgth（x）；
nｈ ＝leｎgth（h）;
n＝nｘ+nｈ—１;
ｆor ｉ=nx＋1:
x（i）= ０;
ｆori=nh＋１:
ｈ（i） = 0；
ｆori＝1:
ｎ
forｊ=1：
H（ｉ,j）＝h（mod（i+n-j，ｎ）+1）;
y= Ｈ*x’;
subｐlot（3， １， 1）;
sｔeｍ（x）；
titlｅ（’随机序列1＇）;
suｂpｌｏt（3,1，2）;
steｍ（h）；
tｉｔlｅ（'
随机序列2'
ｓubpｌot（3,1,3）；
stem（ｙ）;
线性卷积结果’）；
两个序列循环卷积的函数:
cｌearalｌ
x=ones（１，8）；
h＝oneｓ（1,１0）；
n=15；
nx＝ｌengｔh（x）;
nh＝ｌengｔh（h）;
ｉｆ（n〈nx||n＜nh）
fｐrｉntｆ（'
输入圆周卷积的点数不正确’）；
brｅak
foｒk=nh+1:
h（k） = ０;
eｎd
fork=nx+1:
x（k）= 0;
foｒ k=１：
fｏrｌ=1:
H（k,l）=ｈ（ｍod（k+n—l，n）＋１）；
enｄ
y=H*x’；
subpｌｏt（3,1,１）;
steｍ（ｘ）;
随机序列1'
sｕｂｐloｔ（３,１,２）;
sｔem（h）;
tｉtle（’随机序列2'
sｕbplot（3，1, 3）;
stｅm（y）;
ｔｉｔle（’循环卷积结果'
1.线性相关实现的程序及结果
y={８8865 4321 12 3 45 ６ 7｝
循环相关实现的程序及结果
y=｛8 8 8 8 888 88 ８}
3．线性卷积实现的程序及结果
y＝{ 1２ 34 56 ７ ８ 8 87 ６ 5 43 ２ 1}
4.循环卷积实现的程序及结果
当n＝1５时
y=｛3 33４5６７8 ８8 7 6 ５ 4 ３｝
当n＝17时
y=｛1 2 ３ 4５６ ７ 8８ 8 7 6 54 ３2 1}
由上图可知:
15点循环卷积结果与线性卷积的结果是不一致的，但是17点循环卷积结果与线性卷积的结果是一致的.
实验三维纳－霍夫方程的求解
学习使用Matlaｂ实现W-H程序的编写。
一个线性系统,如果它的单位样本响应为ｈ（n）,当输入一个随机信号x（n）:
其中s（ｎ）表示信号,表示噪声，则输出y（n）为
我们希望ｘ（ｎ）通过线性系统h（n）后得到的y（n）尽量接近于s（ｎ），因此称ｙ（n）为s（n）的估计值，用
表示，即
维纳滤波的标准方程
如果我们以
分别表示信号的真值与估计值,而用e（ｎ）表示它们之间的误差
目标:
均方误差
ｍin（MMＳE准则）
上式可看成输出等于现在和过去各输入的加权之和
其中
现在的问题是需要求得使
最小的
，为此,将这式对
求偏导,并令其结果等于0,得
于是
这样就得到维纳滤波的标准方程
ＦIR维纳滤波器
设ｈ（n）是一个因果序列可以用有限长（长度为N）的序列去逼近它,有上述得到Ｗ-H方程的矩阵形式为:
即:
自相关矩阵称，
为x与s的互相关矩阵
这样得到W－H方程的解为:
编写函数解W—Ｈ方程,寻找最优的滤波器，并检验该程序的准确性.
cｌearalｌ;
cｌc;
%输入信号
A=1；
f=10０0；
fs=１0^5；
t= （0:
９99）;
Mlａｇ = 1０0;
ｘ=A*ｃos（２*pｉ*f*t/ｆｓ）;
％给正弦波信号加入信噪比为20dB的高斯白噪声
xｎ=awgn（x，5）；
figure;
subplｏt（２,2,1）
ｐlｏｔ（t,xｎ）
输入信号图像’）
％计算输入信号自相关函数
Ｒxn=xcorr（ｘn,Ｍlag，'
bｉaｓed'
sｕｂｐlot（２，2，２）
ｐloｔ（（—Mlag:
Mlag）,Rxn）
ｔitle（＇输入信号自相关函数’）
%维纳滤波
N=１00;
Rxｎｘ =ｘcorr（xｎ,x,Mlag,’biａsed’）;
ｒxnx= ｚｅros（N,1）;
ｒｘnx（:
） =Ｒxnx（１01:
1０１＋N—1）;
Rxｘ =ｚｅｒos（Ｎ,Ｎ）;
Rxx =diａg（Rxn（101）＊ones（1,N））；
ｆor i=2:
Ｎ
ｃ=Rxn（10１+i）*ｏnes（１,N+1—i）；
Rxx=Ｒxx+dｉag（c，ｉ-1）+diａg（c,—ｉ+1）；
Rxx；
h=ｚeｒos（N，1）；
h＝ｉnｖ（Rxx）*ｒｘnx;
yｎ= filteｒ（ｈ，1,ｘn）；
subplot（2,2,3）
plot（yn）;
ｔｉtlｅ（’经过维纳滤波器后信号信号’）；
Ryn=ｘcorr（ｙn,Mlag,’biａsed'
subｐloｔ（2，２,4）;
plｏt（（—Mｌａg:
Mlag），Ryn）；
title（’经过维纳滤波器后信号自相关函数’）；
四、实验结果
从图中可以看出，滤波后得到的正弦信号仍然有一定的误差，但是输入信号的自相关函数在0点出有明显的噪声成分，通过维纳滤波以后得到的信号的自相关函数在0点处已经的噪声给消除了很多。
实验四
Yｕle-Wａlｋer方程的求解
学习使用Matlab实现Y—W程序的编写.
AR模型的
Yule-Ｗaｌker方程
AR模型,又称为自回归模型，是一个全极点的模型，可用如下差分方程来表示:
（1）
其中
是均值为零、方差为的白噪声序列, p是AR模型的阶数， ａ（ k）， k＝ 1，2,,p 是ｐ阶AＲ模型的参数。
AR 模型系统H（z） 的转移函数为：
（2）
从而得到AR模型的功率谱估计的计算公式：
（3）
由上式可以看出,要利用ＡR模型进行功率谱估计，必须得到模型参数和白噪声序列的方差。
将
（1）式变形有:
（4）
式（4）的矩阵形式为:
（5）
式（ 4）和（5）是AR模型的ARYｕｌe-Wａlker方程．
2.burg算法求解方法
Bｕｒg算法是使序列x（n）的前后向预测误差功率之和：
（6）
最小。
利用Burg法求解ＡＲ模型参数的步骤:
第一步：
由初始条件
根据公式（７）求出反射系数
:
（7）
第二步：
根据序列x（ｎ）自相关函数
，求出阶次m=１时的AR模型参数ａ（1，1）=k1与前后向预测误差功率之和。
第三步：
由式（8）求出前向预测误差
与后向预测误差
，然后由式（7）估计出反射系数ｋ2；
（8）
第四步：
由（9）Levinｓion递推关系,求出阶次ｍ=2时的ＡR模型参数a（2,1）和
;
（9）
第五步:
重复上述过程,直到阶次m＝p，这样就求出了所有阶次的AＲ模型参数。
Burg算法的递推过程是建立在数据序列基础上,避开了序列的自相关函数的估计,所以与自相关法相比, 具有较好的频率分辨率。
已知观测信号,编写函数解Y-W方程,寻找参数系数,并检验该程序的准确性和掌握用法。
clearall；
x ＝ ｃｏs（０：
0。
1:
5０）;
N＝lｅngｔh（x）；
%周期图法
FFＴＸ=fｆt（ｘ）；
POW =（aｂs（FFTX）.＾2）／N;
ｓubｐｌoｔ（1,2，1），plot（1/N：
2*pi/N/２/ｐi:
0。
5，PＯW（１:
N／２））;
ｔiｔle（’周期图求取功率谱’）;
xｌaｂel（’ｆ/Hz’）；
%Ｂｕｒg法的AR谱估计
p=8;
a=zｅros（ｐ,p）;
%初始化前向和后向误差以及γ
e=ｘ；
b=x;
sｉｇma＝0;
ｆorｉ＝1：
N
sigｍａ=sｉgma+x（i）。
ｓｉgma=sigｍa/N;
fｏrn=1:
p
sum１=０；
sum2=0;
forj＝n+1：
sum１=sum1+２*e（ｊ）*b（j－1）；
sｕm2=ｓum2＋ｅ（j）.＾2+ｂ（ｊ-1）.^2;
eｎｄ
a（n,n）=—suｍ1／sum2;
sigma=ｓigma*（1-aｂs（a（n,ｎ））．^2）;
ｉfn>
=2
fori＝1:
ｎ—１
a（n,i）＝a（ｎ—1,i）+a（n,ｎ）*ａ（ｎ-１，n-i）；
enｄ
ｆoｒｊ=n＋1:
c（j）=ｅ（j）+a（n,ｎ）*b（j—1）；
d（ｊ）=b（j—１）+a（n,n）*ｅ（j）;
end
e=c;
b=d;
％计算并输出功率谱
form=1:
ｓum=0；
fｏrｎ＝１:
sum＝ a（p，n）*eｘｐ（—sqrｔ（-1）*2*pi＊n*m/N）＋ｓum；
eｎd
PＯW２（m）=ｓｉgｍa/（aｂs（１+sum）.^2）;
sｕbpｌot（１,2，2）,plot（1／N：
2*ｐi/N/2/pi:
０。
5，PＯW2（1：
N/2））；
tｉtlｅ（＇AＲ模型谱估计法求取功率谱＇）；
xlabｅl（’f/Hz’）；
四．实验结果及分析
实验中输入的信号为余弦信号,理想情况下其功率谱是在余弦信号频率上的一个冲击函数．从实验的结果图可以看出,用AR模型估计的功率谱同用周期图法估计的功率谱一样,说明了用该方法所计算的功率谱的准确性．
实验五 自适应噪声抵消算法的软件设计与实现
学习使用MATLAB编写LＭS自适应滤波器，以及如何在生物医学信号中进行应用。
1、自适应干扰抵消的原理
图1自适应干扰抵消原理图
图1所示的是自适应干扰抵消器的基本结构。
期望信号d（n）是信号与噪声之和,即d（ｎ）＝x（n） +N（n）,自适应处理器的输入是与N（ｎ）相关的另一个噪声N＇（n）.当ｘ（n）和N（n）不相关时，自适应处理器将调整自己的参数,使y（n）成为N（n）的最佳估计
。
这样，e（n）将逼近信号ｘ（n）,且其均方差E［e2（ｎ）]为最小．噪声N（n）就得到了一定程度的抵消.
2、LMS自适应滤波算法
图２单输入自适应线性组合器
LMＳ算法使用的准则是使滤波器的期望输出值和实际输出值之间的均方误最小化的准则,即使用均方误差来做性能指标.
自适应滤波的结果如图2所示。
各符号的意义是：
x（n）输入信号,ｙ（n）为滤波器的输出,d（n）为y（n）想要趋近的理想信号，d（ｎ）是已知的,e（n）为误差信号.滤波器均方误差可表示为:
设自适应滤波器的输入矢量为:
加权矢量（即滤波器参数矢量）为：
滤波器的输出为：
误差信号为期望输出ｄ（n）与滤波器实际输出之间的误差，即
（1）
ＬMS算法是取单个误差样本的平方e2（n）的梯度作为均方误差梯度的估计,由式
（1）可得梯度矢量的估计为：
由此得到一个新的权矢量递推公式,即ＬMS算法递推公式为:
已知观测信号，编写LＭS滤波器,对该信号进行滤波处理,检验该程序的准确性和掌握ＭATLAB自带函数的用法。
clｅar all；
d=ｓiｎ（0:
Ｍ =leｎgｔh（d）;
％设定滤波器的长度为15
N=15；
noise＝ ０.2＊sin（5０＊（0：
0.１:
5０））;
%噪声信号
x＝d+noｉse;
ｕ=0。
0１；
%初始化滤波器的参数为全0的矩阵
w=zeros（1,N）；
for ｎ=N：
Ｍ
ｙ（n）=0;
fｏri＝1：
ｙ（ｎ）=ｙ（ｎ）+x（ｎ—i+1）*w（i）；
e（n）=d（ｎ）—ｙ（n）；
ｆor i=1:
Ｎ ％进行滤波器参数的调整
ｗ（i）=w（i）+２*u*e（n）*x（n—i＋1）；
subpｌot（2,１，1）；
plot（y）;
titlｅ（’滤波处理后得到的信号’）
suｂpｌot（2，1,2）;
plot（x）;
tiｔlｅ（’观测到的信号'
实验中输入的原始信号为正弦信号，噪声为在正弦信号的基础上加入一个５0倍频率的噪声信号,观测到的信号是这两个信号的叠加,利用LMS算法进行滤波处理,从从实验结果中可以看出,在LMS的加权参数稳定之后的输出结果比较理想。
升级会员 