数字信号处理II
——随机信号的功率谱估计方法
一、实验目的
1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。 2．观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳
定性、主瓣宽度和旁瓣效应的影响。 3．学习使用FFT 提高谱估计的运算速度。 4．体会非参数化功率谱估计方法的优缺点。
二、实验原理与方法
假设信号()x n 为平稳随机过程，其自相关序列定义为:
}*
()()()m E x n x n m φ+
(0.1)
其中{}E ?表示取数学期望，{}*
?表示取共轭。根据定义，()x n 的功率谱密度()P w 与自相关序列()m φ存在如下关系：
j m
m P m e
ωωφ+∞
-=-∞
=
∑
(0.2)
1()()2j m
m P e
d π
ωπ
φωωπ
-=
?
(0.3)
然而，实际中我们很难得到准确的自相关序列()m φ，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。目前常用的线性谱估计方法有两种：自相关函数法和周期图方法，本实验将对这两种方法分别予以讨论。
1．自相关函数法
假设已知随机信号()x n 的N 个观测样本，则其自相关序列可以用下式进行估计:
||1
*
1?()()()
||1||
N m n m x n x n m m N N m φ--==
+≤--∑
(0.4)
当仅使用长度为2M-1的自相关序列时，对其进行傅立叶变换即可得到功率谱估计如下:
1
1
??()()M j m
m M P
m e
ωωφ--=-+=∑
(0.5)
其中M 为加窗长度，R e ()c
M
W m 为矩形窗函数，定义如下: R e 1,||()0,
||c
M
m M W
m m M
≥? (0.6)
因此，?()P
w 在一定程度上可以看作是“真正的功率谱()P w ”与窗函数傅立叶变换的卷积。
矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣，旁瓣效应使那些处于旁瓣附近功率较小的频率分量被淹没掉。为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为：
0.540.46cos ,||()0,||H am
M m m M W m M
m M
π?
+
(0.7)
这种窗函数可以有效的抑制旁瓣，但此时主瓣宽度增大，从而降低了谱估计的分辨率。这种主瓣和旁瓣之间的矛盾在非参数化功率谱估计方法中是无法解决的。
2．周期图方法
假设已知随机信号()x n 的N 个样本，利用周期图方法，信号()x n 的功率谱估计为:
1
2
1?()|()|N j m
n P
x n e
N
ωω--==∑
(0.8)
利用上述方法得到的谱估计方差与信号()x n 的功率谱平方2[()]P w 成正比，为了减小它的方差，可以将信号序列进行分段处理，然后再求各分段结果的平均，这就是平均周期图方法，即Bartlett 方法,它特别适用于FFT 直接计算功率谱估值。 （1）Bartlett 平均周期图方法
将一个随机序列()x n (0)n N ≤≤分成K 段，每段长度为L ，各段之间互不重叠，因而N
L K
=。可以想到，第i 段的信号序列可表示为：
()((1))
0n<,1i x n x n i L L i K
=+-≤≤≤ (0.9)
每一段的周期图又可写成：
120
1?()|()|1L j m i i n I x n e i K
L ωω--==≤≤∑ (0.10)
于是，功率谱估计定义为：
1
1??()()K
i
i P
I K
ωω==∑
(0.11)
对于固定的记录长度来讲，分段数K 增大可使谱估计的方差减小，但是由于L 的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率也是矛盾的。
除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以得到较平滑的谱估计，当然，相应的分辨率也有所下降。
（2）平滑平均周期图方法
这是一种改进的Bartlett 周期图方法，将长度为 N 的平稳随机信号序列()x n 分成K 段，每段长度为L ，即/L N K
=。这里在计算周期图之前，先用窗函数
L W (n)给每段序列(n)i x 加权，K
个修正的周期图定义为：
2
-1
-=0
1?()=
(n)W
(n)e
1L W
j m
i i
L
n I x i K LU
ωω≤≤∑ (0.12)
其中U 表示窗函数序列每个样本的平均能量，即：
2
-1
=0
1=
W
(n)
L L
n U L
∑ (0.13)
在这种情况下，功率谱估计可按下面表达式给出：
=1
1??
()=()
K
W i
i P
I K
ωω∑ (0.14)
本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计：
(n+)
=1
(n)=+(n)s
i i N j i i x a e
w ω?∑ (0.15)
其中N s 为正弦个数，i ω，i ?和i a 分别为第i 个正弦信号的数字频率、相位和
幅度，i ?随机分布在[0,2)π之间，(n)w 为零均值，
方差等于2
w σ的复高斯白噪声。 三、实验内容及步骤
1、仔细阅读有关功率谱估计的内容，按照图1给出的框图编制自相关函数法谱
估计程序。运行程序，输入2
=100,=10,=1,=0.6,=0,=0,s i i w
N M N ωπ?σ选择矩形窗。观察谱峰位置是否正确（注意：由于窗效应可能引起谱估计的非正定）。
图 1实验流程框图
图2 M=10时，矩形窗的谱估计
从图中可以看出谱峰位置在0.6*pi处，这个位置是正确的。
2、观察并记录参数变化对谱估计性能的影响。
(1)改变M = 5，其它输入同步骤1，观察功率谱估计的主瓣宽度和旁瓣大小随自
相关序列长度的变化情况。
图3 M=5时，矩形窗的谱估计
观察M=10和M=5的仿真结果，可以看出M越小，则主瓣宽度越大，分辨率越低，谱峰高度也降低了；旁瓣越少，旁瓣最大幅度也减小；主瓣与旁瓣的幅值比、宽度比基本没有改变。
(2) 选择窗函数为Hamming 窗，其它输入同步骤1，观察不同的窗函数对谱估计性能的影响。（一般来讲，与矩形窗相比，加Hamming 窗，主瓣宽度增加一倍，但旁瓣减小许多）。
图4 M=10时， hamming 窗的谱估计
与矩形窗相比，选用汉明窗，功率谱的主瓣宽度增加了近一倍，旁瓣也减少了许多，分辨率降低了。
(3) 改变1=/4?π，其它输入同步骤 1，观察初始相位的变化对谱估计性能的影响。
图5 4/1π?=，矩形窗的谱估计
由图可知改变初始相位，对频谱估计的结果没有影响。
(4) 改变2=1
σ，其它输入同步骤1，观察信噪比的变化对谱估计性能的影响。
w
σ=，矩形窗的谱估计
图621
w
当噪声增加，信噪比减小时，会使谱估计得到的结果的精确性减小，而且增加的噪声越大，对谱估计影响就越大，
(5) 改变2
=10,=1
Nσ，其它输入同步骤1，结合（4）的内容，观察数据长
w
度及信噪比对谱估计性能的影响。
图7 N=10，2
1w σ=，矩形窗的谱估计
由图知，在N 较大时，改变信噪比对谱估计影响比较小；而当N 较小时，改变信噪比对谱估计的影响较大。在上图中我们可以看出此时频谱峰值已偏离了0.6pi ，同时旁瓣的幅值增大了许多。
3、运行程序，输入2
121212=100,=2,=0.6,=0.8,==0,==1,=0,s w
N N a a ωπωπ??σ选择矩形窗，调整自相关序列长度M ，使得两个正弦频率分量临界分辨出来，记录此时M 值，并绘出功率谱图。同样，在加 Hamming 窗的情况下，记录使两个正弦频率分量临界分辨出来的M 值，并绘出此时的功率谱图。
图8 当M=7时，矩形窗的谱估计
图9 当M=8时，矩形窗的谱估计
由上面两图可知，在M=8时，加矩形窗可以将两个正弦频率分量临界分辨出来。
图10 当M=9时，汉明窗的谱估计
图11 当M=10时，汉明窗的谱估计
由上两图可知，在M=10时，加汉明窗可以将两个正弦频率临界分辨出来。
4、运行自相关函数法谱估计程序，输入2
=100,=0,=1,M=10,
N N 选择矩
s w
形窗，观察利用自相关函数法得到白噪声信号谱估计。改变M =3、20，观察M 的变化对白噪声谱估计的影响。
图12 当M=3时，加矩形窗的白噪声谱估计
图13 当M=10时，加矩形窗的白噪声谱估计
图14 当M=20时，加矩形窗的白噪声谱估计
由上面三幅图知，M 越大，加矩形窗的谱估计各个分量变得越清晰，更容易分辨。
5．根据图1给出的框图，编制周期图法谱估计程序，子程序FFT 可以直接调用Matlab 中的函数。运行程序，输入=100,=128,=10,F N N K 选择窗函数为矩
形窗，输入2
111=1,=0.6,=0,=1,=0,s w
N a ωπ?σ观察谱峰位置是否正确，并与步骤 1 结果比较。
图15 N/K=10时的周期图法谱估计
此图峰值在0.6*pi 处，位置正确，并且由对比可知，主瓣峰值明显下降了，而且其宽充增加，旁瓣数目有所减少。
6．利用周期图法重复 2，3，4 的内容，这里=/K L N 相当于自相关函数法中的M ，观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。
（1）重复步骤2(1)，选矩形窗可得如下图所示的周期图法谱估计
图16 N/K=5时，矩形窗的Bartlett周期图谱估计
对比每2(1)步中结果可知，用Bartlett周期图谱估计的主瓣在N/K减小时，主瓣宽度变大，分辨率降低，峰值下降，而且旁瓣数目减少，其规律和自相关函数法中减小M时，规律一致。
（2）重复步骤2(2)，选Hamming窗可得如下图所示的周期图法谱估计
图17 N/K=5时，汉明窗的Bartlett周期图谱估计
对比第2(2)步中结果可知与矩形窗相比，汉明窗知主瓣宽度加倍，旁瓣减小许多。Bartlett周期图谱估计和自相关谱出现的规律相似。
（3）重复步骤2(3)得到如下图所示的周期图法谱估计
图18 4/1π?=，N/K=10时矩形窗的Bartlett 周期图谱估计
对比每2(3)可知，φ变化对谱估计没有太大影响。 （4）重复步骤2(4)得如下图所示的周期图法谱估计
图19 取12
=w σ，N/K=10时矩形窗的Bartlett 周期图谱估计
对比第2(4)步可知，增大高斯白噪声的方差时，对谱估计的峰值没有影响，但却使旁瓣增大。
（5）重复步骤2(5)得如下图所示的周期图法谱估计
图20 N=50，12
=w σ时矩形窗的Bartlett 周期图谱估计
由图可知，当N 减小，同时增大噪声方差时，矩形窗的Bartlett 周期图谱估计的
主瓣宽度增加。
（5）用周期图法重复步骤3，输入100,2,N N S ==10.6,ωπ=20.8,
ωπ=12120,1a a φφ====，02
=w σ，加矩形窗，调整段数K ，其仿真结果如下图所示；
加Hamming 窗时，仿真如下图所示
图21 K=20时，矩形窗的Bartlett 周期图谱估计
图22 K=25时，矩形窗的Bartlett周期图谱估计
由上两幅图可知，加矩形窗时，K=20时可以将两个正弦频率分量临界分辨出来，但峰值位置有所变化。
图23 K=20时，汉明窗的Bartlett周期图谱估计
图24 K=25时，汉明窗的Bartlett 周期图谱估计
同样由上面两图可知，加汉明窗时，当K=20时，可以将两个正弦频率分量临界分辨出来。同样其频率中心也发生了变化。
（6）利用周期图方法重复步骤4，加矩形窗，2
100,0,1w
N N S σ===，当N/K=10、4、20，仿真结果分别如下图所示
图25N/K=4时，矩形窗的Bartlett 周期图谱估计
图26N/K=10时，矩形窗的Bartlett周期图谱估计
图27N/K=20时，矩形窗的Bartlett周期图谱估计
由上图可知，N/K越大，谱估计的各个分量变得越清晰，更容易分辨。
综上可右，对比自相关函数法的谱估计，由步骤5、6的仿真结果，我们发现周期图法谱估计减小了峰值的幅度，加宽主瓣宽度，但是周期图使方差减小，从而得到一致的谱估计。周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计。
四、思考题
1. 证明式（1.5）可以表示为：
-1
-=0
???()=2R eal[(m )e ]-(0),M j m m P ωωφφ∑其中Re al[]?表示取实部。
证明：1
1
1
1
????()()()()M M j m
j m
j m
m M m M m P
m e
m e
m e
ωωωωφφφ-----=-+=-+===
+
∑
∑
∑
k=-m 由得：
-1
1
k 00
???()(-)()M M j k
j m
m P
k e m e
ωωωφφ--===+∑∑?-(0)φ
由自相关函数的对称性得：
-1
1
k 00
???()()()M M j k
j m
m P
k e m e
ωωωφφ--===+∑∑?-(0)φ
即是：1
??()()(+)M j k
j m
m P
m e
e
ωωωφ--==∑?-(0)φ
即得：-1
-=0
???()=2R eal[(m )e ]-(0)M j m m P ωωφφ∑
2. 已知实信号(n)x 的 N 个样本(0),(1),...,(1),x x x N -，可以定义()A P ω的估计如下：
1
1
??()(),N j m
A
m N P m e
ωωφ--=-+=∑
其中||1
*
1?()()()N m i m x i x i m N
φ
--==+∑
试证明：2
1
1?()()N j n
A
n P x n e
N
ωω--==∑。
证明：1||1
1
*10
1
1??()()()()N m N N j m
j m
A
m N i m N P m e
x i x i m e
N
ωωωφ------=-+==-+==
+=
∑
∑
∑
1||
1
*
()
1
1[()]()N m N j i j m i i m N x i e
x i m e
N
ωω-----+==-++∑
∑
2
1
1
*
11[()]()()N N j n j n
j n
n n x n e
x n e
x n e
N
N
ωωω-----===
=
∑∑
五、实验程序
1、自相关函数法（信号不为零）
clc;
clear all; n=100;
m=10;
ns=2;
w1=0.6*pi;
w2=0.8*pi;
fi1=0;
fi2=0;
a1=1;
a2=1;
sigma2=0;
l=256;
sn1=zeros(n,1);
sn2=zeros(n,1);
xn=zeros();
sn1=a1*exp(i*(w1*(1:n)'+fi1));
sn2=a2*exp(i*(w2*(1:n)'+fi2));
wn=sqrt(sigma2).*(randn(n,1));
xn=sn1+sn2+wn;
xxn=xcorr(xn,'unbiased'); %产生信号long=2*m-1; %矩形窗Rectangle=window(@rectwin,long);
p1=xxn(n-m+1:n+m-1).*Rectangle;
p2=fft(p1,l);
p3=abs(p2(1:l/2));
long1=2*m-1; %汉明窗Hammwin=window(@hamming,long1);
p11=xxn(n-m+1:n+m-1).*Hammwin;
p22=fft(p11,l);
p33=abs(p22(1:l/2));
figure(1);
stem((0:l/2-1)/(l/2),p3,'.');
title('矩形窗功率谱估计')
grid on;
figure(2)
stem((0:l/2-1)/(l/2),p33,'.');
title('汉明窗功率谱估计')
grid on;
2、自相关函数法（白噪声）
clear all;
n=100;
m=20;
w=0.6*pi;
sigma2=1;
l=256;
现代信号处理作业 实验题目： 设信号)()8.0cos(25.0)47.0cos()35.0cos()(321n v n n n n x ++++++=θπθπθπ，其中321,,θθθ是[]ππ,-内的独立随机变量，v(n)是单位高斯白噪声。 1.利用周期图法对序列进行功率谱估计。数据窗采用汉明窗。 2.利用BT 法对序列进行功率谱估计，自相关函数的最大相关长度为M=64,128,256,512采用BARTLETT 窗。 3.利用Welch 法对序列进行功率谱估计，50%重叠，采用汉明窗，L=256,128,64。 4.利用Burg 法对序列进行AR 模型功率谱估计，阶数分别为10，13. 要求每个实验都取1024个点，fft 作为谱估计，取50个样本序列的算术平均，画出平均的功率谱图。 实验原理： 1）。周期图法： 又称间接法，它把随机信号的N 个观察值x N (n)直接进行傅里叶变换，得到X N (e jw ),然后取其幅值的平方，再除以N ，作为对x （n ）真实功率谱的估计。 2^ )(1)(jw e X N w P N per = ， 其中∑-=-=1 )()(N n jwn N jw N e n x e X 2）。BT 法： 对于N 个观察值x(0),x(1),。。。,x(N-1),令x N (n)=a(n)x(n)。计算r x （m ）为
∑--=-≤+= m N n N N x N m m n x n x N m r 10 1),()(1 )(,计算其傅里叶变换 ∑-=--≤= M M m jwm x BT N M e m r m v w P 1 ,)()()(^ ^ ，作为观察值的功率谱的估计。 其中v(m)是平滑窗。 3)。Welch 法： 假定观察数据是x(n),n=0,1,2...,N-1，现将其分段，每段长度为M,段与段之间的重叠为M-K,第i 个数据段经加窗后可表示为 1,...,1,0 )()()(-=+=M i iK n x n a n x i M 其中K 为一整数，L 为分段数，该数据段的周期图为 2)(1)(^w X MU w P i M i per =，其中∑-=-=1 0)()(M n j w n i M i M e n x w X 。由此得到平均周期图为 ∑-==10 ^_ )(1)(L i i per w P L w P 。其中归一化U 取∑-== 10 2 )(1M n n a M U 。 4）。Burg 法： 在约束条件下，使得)(2 1^^^ b f ρρρ+=极小化，其中，约束条件是它所得到的 各阶模型解要求满足Levison 递归关系。 仿真结果： 1.周期图法
功率谱估计方法的比较 摘要： 本文归纳了信号处理中关键的一种分析方法, 即谱估计方法。概述了频谱估计中的周期图法、修正的协方差法和伯格递推法的原理，并且对此三种方法通过仿真做出了对比。 关键词：功率谱估计；AR 模型；参数 引言： 谱估计是指用已观测到的一定数量的样本数据估计一个平稳随机信号的谱。由于谱中包含了信号的很多频率信息，所以分析谱、对谱进行估计是信号处理的重要容。谱估计技术发展 渊源很长，它的应用领域十分广泛，遍及雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域，其容、方法都在不断更新，是一个具有强大生命力的研究领域。谱估计的理论和方法是伴随着随机信号统计量及其谱的发展而发展起来的，最早的谱估计方法是建 立在基于二阶统计量, 即自相关函数的功率谱估计的方法上。功率谱估计的方法经历了经典谱估计法和现代谱估计法两个研究历程，在过去及现在相当长一段时间里，功率谱估计一直占据着谱估计理论里的核心位置。经典谱估计也成为线性谱估计，包括BT 法、周期图法。现代谱估计法也称为非线性普估计，包括自相关法、修正的协方差法、伯格（Burg ）递推法、特征分解法等等。 原理： 经典谱估计方法计算简单，其主要特点是谱估计与任何模型参数无关，是一类非参数化的方法。它的主要问题是：由于假定信号的自相关函数在数据的观测区间以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，经典法的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。现代谱估计方法使用参数化的模型，他们统称为参数化功率谱估计，由于这类方法能够给出比经典法高得多的频率分辨率，故又称为高分辨率方法。下面分别介绍周期图法、修正的协方差法和伯格递推法。修正的协方差法和伯格递推法采用的模型均为AR 模型。 （1）周期图法 周期图法是先估计自相关函数, 然后进行傅里叶变换得到功率谱。假设随机信号x(n)只观测到一段样本数据，n=0, 1, 2, …, N-1。根据这一段样本数据估计自相关函数，如公式（1） 对（1）式进行傅里叶变换得到（2）式。 ∑--=+=1||0 *) ()(1 )(?m N n xx m n x n x N m r
第3章 平稳随机过程的谱分析 付里叶变换是处理确定性信号的有效工具，它信号的频域内分析处理信号，常常使分析工作大为简化。 对于随机信号，是否也可以应用频域分析方法？付里叶变换是否可引入随机信号中？ 3.1 随机过程的谱分析 3.1.1 回顾：确定性信号的谱分析 )(t f 是非周期实函数， )(t f 的付里叶变换存在的充要条件是： 1．)(t f 在),(∞-∞上满足狄利赫利条件； 2．)(t f 绝对可积： +∞
3.1.2 随机过程的功率谱密度 一、样本函数的平均功率 问题1：由于付里叶变换是针对确定性函数进行的，在处理随机过程)(t X 时，取 )(t X 的一个样本函数)(t x （在曲线族中取某一曲线）来进行付里叶分 析。 问题2：随机过程)(t X 的样本函数)(t x 一般不满足付里叶变换的条件，它的总能 量是无限的，需考虑平均功率。 若随机过程)(t X 的样本函数)(t x 满足 +∞<=? -∞→T T T dt t x T W 2 )(21 lim W 称为样本函数)(t x 的平均功率。 对于平稳过程，其样本函数的平均功率是有限的。 二、截取函数 对于)(t X 的一个样本函数)(t x ，在)(t x 中截取长为T 2的一段，记为)(t x T ， 它满足： ???? ?≥<=T t T t t x t x T 0 ) ()( 称)(t x T 为)(t x 的截取函数。 三、截取函数的付里叶变换 0>T ，取定后，)(t x T 的付里叶变换一定存在： ??--+∞ ∞--==T T t j t j T T dt e t x dt e t x X ωωω)()()( 其付里叶逆变换为： ? +∞ ∞ -= ωωπ ωd e X t x t j T T )(21 )( 其帕塞瓦（Parseval ）等式为 ? ? ? +∞ ∞ --+∞ ∞ -= =ωωπ d X dt t x dt t x T T T T 2 2 2 )(21 )()(
1 随机信号的经典谱估计方法 估计功率谱密度的平滑周期图是一种计算简单的经典方法。它的主要特点是与任 何模型参数无关，是一类非参数化方法[4]。它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。本章主要介绍了周期图法、相关法谱估计（BT ）、巴特利特(Bartlett)平均周期图的方法和Welch 法这四种方法。 2.1 周期图法 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。周期图法[5]包含了下列两条假设： 1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段 )(n x N 来估计该随机序列的功率谱。这当然必然带来误差。 2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数)(m R x 时将 )(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外 x(n)是全零序列，这种处 理方法显然与周期图法不一样。 但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。通过比较我们发现：如果相关法中M=N ，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N 时相关法的特例。因此相关法和周期图法可结合使用。 2.2 相关法谱估计（BT ）法
功率谱估计及其MATLAB仿真 詹红艳 （201121070630控制理论与控制工程） 摘要:从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。 关键词:功率谱估计;周期图法;AR参数法;Matlab Power Spectrum Density Estimation and the simulation in Matlab Zhan Hongyan （201121070630Control theory and control engineering） Abstract:Mainly introduces the principles of classical PSD estimation and modern PSD estimation,discusses the characteristics of the methods of realization in Matlab.Moreover,It gives an example of each part in realization using Matlab functions. Keywords:PSDPstimation,Periodogram method,AR Parameter method,Matlab 1引言 现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。 功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。 Matlab是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件,人称矩 阵实验室,它集数值分析、矩阵运算、信号处理和图形显示于一体,构成了一个方便的、界面友好的用户环境,成为目前极为流行的工程数学分析软件。也为数字信号处理进行理论学习、工程设计分析提供了相当便捷的途径。本文的仿真实验中,全部在Matlab6.5环境下调试通过;随机序列由频率不同的正弦信号加高斯白噪声组成。 2经典功率谱估计 经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1: Fs=500;%采样频率 n=0:1/Fs:1;
随机信号 利用经典谱估计法估计信号的功率谱
作业综述： 给出一段信号“asd.wav”,利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数。采用MATLAB语言，利用MATLAB语言强大的数据处理和数据可视化能力，通过GUI的对话框模板，使操作更为简便！在一个GUI界面中，同时呈现出不同方法产生出的功率谱。 这里给出了几种不同的方法：BT法，周期图法，平均法以及Welch法。把几种不同方法所得到的功率谱都呈现在一个界面中，便于对几种不同方法得到的功率谱作对比。 一．题目要求 给出一段信号及采样率，利用经典谱估计法估计出信号的功率谱。 二．基本原理及方法 经典谱估计的方法，实质上依赖于传统的傅里叶变换法。它是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有BT法，周期图法，平均法以及Welch法。 1. BT法(Blackman-Tukey) ●理论基础： (1)随机序列的维纳-辛钦定理 由于随机序列{X(n)}的自相关函数Rx(m)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为，则可将随机序列的自相关函数用连续时间函数表示为 等式两边取傅里叶变换，则随机序列的功率谱密度 (2)谱估计 BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值。即 其中可有式得到。 2. 周期图法 ●理论基础： 周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。在前面我们已知，各态历经的连续随机过程的功率谱密度满足
式中 是连续随机过程第i 个样本的截取函数 的频谱。对应在随机序列中则有 由于随机序列中观测数据 仅在 的点上存在，则 的N 点离散傅里叶变换为： 因此有随机信号的观测数据 的功率谱估计值（称“周期图”）如下： 由于上式中的离散傅里叶变换可以用快速傅里叶变换计算，因此就可以估计出功率 谱。 3.平均法： 理论基础： 平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的，有一个定理：如果 ， ， ， 是不相关的随机变量，且都有个均值 及其方差 ，则可以证明它们的算术平均的均值为 ，方差为 。 由定理可见：具有 个独立同分布随机变量平均的方差，是单个随机变量方差的 ， 当 时，方差 ，可以达到一致估计的目的。因此，将 个独立的估计量经过算术 平均后得到的估计量的方差也是原估计量方差的 。 平均图法即是将数据 ， ， 分段求周期图法后再平均。例如，给定N=1000个数据样本（平均法适用于数据量大的场合），则可以将它分成10个长度为100的小段，分别计算每一段的周期图 ()()2 1001100,100(1) 1 ,1,2,```,10100 l j l n l G w X e l ω-=-= =∑ 然后将这10个周期图加以平均得谱估计值： ()() 10 100100,1 110l l G w G w ==∑ 由于这10小段的周期图取决于同一个过程，因而其均值相同。若这10个小段的周期图是统计独立的，则这10个小段平均之后的方差却是单段方差的 。
功率谱估计性能分析及Matlab 仿真 1 引言 随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。 信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。 2 经典功率谱估计 经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。 周期图法( Periodogram ) Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。 取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换 1 ()()N j j n N n X e x n e ω ω---==∑ 然后进行谱估计
信号的功率谱分析 1、功率谱密度函数的定义 对于随机信号)(t x ，由于其任一样本函数都是时间的无限的函数，一般不能满足傅里叶变换的存在条件（即积分?∞ ∞-dt t x )(必须收敛）。如果将样本函数取在一个有限区间]2 ,2[T T -内，如图所示，令在该区间以外的0)(=t x ，则积分?∞ ∞-dt t x )(收敛，满足傅里叶变换条件，变换后用功率谱密度函数表示。 2、功率谱密度函数（又称功率谱）的物理意义 是在频域中对信号能量或功率分布情况的描述。功率谱表示振动能量在频率 域的分解，其应用十分广泛。功率谱的横坐标是频率，纵坐标是实部、虚部的模 的平方。 功率谱密度函数作为随机信号在频域内描述的函数。对于随机信号而言，它 不存在频谱函数，只存在功率谱密度函数（功率大小在频谱中反映为频谱的面 积）。时域中的相关分析为在噪声背景下提取有用信息提供了途径。功率谱分析 则从频域提供相关技术所能提供的信息，它是研究平稳随机过程的重要方法。 3．功率谱密度函数的应用 (1)结构各阶固有频率的测定 工程结构特别是大型结构(如高层楼房、桥梁、 高塔和重要机械设备等)要防止共振引起的破坏，需要测定其固有频率。如果对 结构加以激励(或以大地的脉动信号作为激励信号)，即可测定结构的响应(振动信 号)，再对响应信号作自功率谱分析，便可由谱图中谱峰确定结构的各阶固有频
率。 (2)利用功率谱的数学特点求取信号传递系统的频率响应函数。 (3)作为工业设备工作状况的分析和故障诊断的依据 根据功率谱图的变化，可以判断机器设备的运转是否正常。同时．还可根据机器设备正常工作和不正常工作时，振动加速度信号的功率谱的差别，查找不正常工作时，功率谱图中额外谱峰产生的原因以及排除故障的方法。 自功率谱密度函数 定义及其物理意义 假如)(t x 是零均值的随机过程，即0=x μ（如果原随机过程是非零均值的，可以进行适当处理使其均值为零）又假设)(t x 中没有周期分量，那么当∞→τ， 0)(→τx R 。这样，自相关函数)(τx R 可满足傅里叶变化的条件∞?∞ ∞- ττd R x )(。 )(τx R 傅里叶变换)(f S x ，ττπτd e R f S j x x 2)()(-∞∞-? =和逆变换 df e f S R j x x πττ2)()(-∞ ∞ -?=定义)(f S x 为)(t x 的自功率谱密度函数，简称自谱或自功率谱。 出此可见， )(f S x 曲线下和频率轴所包围的面积就是信号的平均功率，)(f S x 就是信号的功率密度沿频率轴的分布，故称)(f S x 为自功率谱密度函数。自功率谱密度函数和幅值谱的关系为 2)(1lim )(f X T f S x T x →=利用这一种关系，就可以通过直接对时域信号作傅里叶变换来计算功率谱。自相关分析可以有效地检测出信号中有无周期成分。自功率谱密度也能用来检测信号中的周期成分。周期信号的频谱是脉冲函数，在某特定频率上的能量是无限的。但是在实际处理时，用矩形窗函数对信号进行截断，这相当于在频域用矩形窗函数的频谱sinc 函数和周期信号的频谱——δ函数实行卷积，因此截断后的周期函数的频谱已不再是脉冲函数，原来为无限大的谱线高度变成有限长，谱
随机信号的功率谱密度估计和相关函数
随机信号的功率谱密度估计和相关函数 1．实验目的 了解估计功率谱密度的几种方法，掌握功率谱密度估计在随机信号处理中的作用。 ⒉实验原理 随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。功率谱密度简称为功率谱，是自相关函数的傅里叶变换。对功率谱密度的估计又称功率谱估计。 1.线性估计法（有偏估计）：线性估计方法是有偏的谱估计方法，谱分辨率随数据长度的增加而提高。包括自相关估计、自协方差法、周期图法。 2.非线性估计（无偏估计）：非线性估计方法大多是无偏的谱估计方法，可以获得高的谱分辨率。包括最大似然法、最大熵法 ⒊实验任务与要求 1. 所有功能均用matlab仿真。 2. 输入信号为：方波信号+n(t)，方波信号信号基频1KHz，幅值为1v，n(t)为白噪声。 3. 编写自相关估计法、自协方差法、周期图法、最大似然法、最大熵法的matlab 程序。正确的运行程序。 4. 必须用图示法来表示仿真的结果。对几种功率谱估计的方法进行比较分析，发现它们各自有什么特点？。 5. 按要求写实验报告。 4.Matlab程序如下： 生成输入信号： clear; fs=1024;%设采样频率为1024 n=0:1/fs:1; N=length(n); W=2000*pi;%因方波频率F=1000HZ所以角频率W=2000pi X1n=square(W*n);%方波信号 X2n=randn(1,N);%白噪声信号 xn=X1n+X2n; %产生含有噪声的信号序列XN subplot(3,1,1) plot(n,xn); xlabel('n') ylabel(…输入信号?) %绘输入信号图
1. （10分）随机变量12,X X 彼此独立，且特征函数分别为12(),()v v φφ，求下列随机变量的特征函数： （1） 122X X X =+ （2）12536X X X =++ 解：（1）() 121222()jv X X jvX jv X jvX X v E e E e E e e φ+???? ??===?????? ? 12 21212()(2)jvX jv X X X E e E e v v φφ????=????和独立 （2）() 1212536536()jv X X jv X jv X jv X v E e E e e e φ++???? ==????? ? 12536 12jv X jv X jv X X E e E e E e ?????? ??????和独立 6 12(5)(3)jv e v v φφ= 2. （10分）取值()1,1-+，概率[0.4,0.6]的独立()半随机二进制传输信号()X t ，时隙长度为T ，问： （1） 信号的均值函数()E X t ????； （2） 信号的自相关函数(),X R t t τ+； （3） 信号的一维概率密度函数();X f x t 。 解：（1）()10.410.60.2E X t =-?+?=???? (2) 当,t t τ+在同一个时隙时： []222(,)()()[()]10.6(1)0.41X R t t E X t X t E X t ττ+=+==?+-?= 当,t t τ+不在同一个时隙时： [][][](,)()()()()0.20.20.04 X R t t E X t X t E X t E X t τττ+=+=+=?= （3）()()();0.610.41X f x t x x δδ=-++ 3. （10分）随机信号0()sin()X t t ω=+Θ，()()0cos Y t t ω=+Θ，其中0 ω为常数，Θ为在[]-,ππ上均匀分布的随机变量。
频谱分析（也称频率分析），是对动态信号在频率域内进行分析，分析的结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。 功率谱 频谱和功率谱有什么区别与联系？ 谱是个很不严格的东西，常常指信号的Fourier变换， 是一个时间平均（time average）概念 功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。有两个重要区别： 1。功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列） 2。功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。 功率谱是个什么概念？它有单位吗? 随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。 功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。 一是用相关函数的傅立叶变换来定义谱密度；二是用随机过程的有限时间傅立叶变换来定义谱密度；三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零，这样才能保证相关函数在时差趋向于无穷时衰减，所以lonelystar说的不全对，光靠相关函数解决不了许多问题，要求太严格了；对于第二种方式，虽然一个平稳随机过程在无限时间上不能进行傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换（FFT）估计谱密度的依据；第三种方式是根据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进行重构，在依靠正交性来建立的。
功率谱估计浅谈 摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计 前言 功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。 周期图、自相关法及其改进方法(Welch)为经典(非参数)谱估计方法, 其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题,特别是在数据记录很短的情况下,这一问题尤其突出。以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性,可实现信号检测或信噪分离,对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义,并广泛应用于通信、自动控制、地球物理等领域。在现代参数法功率谱估计方法中,比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关,对短数据具有更强的适应性,从而弥补了经典谱估计法的不足,但其也有一些自身的缺陷。 下面就给出这两类谱估计的简单原理介绍与方法实现。 经典谱估计法 经典法是基于传统的傅里叶变换。本文主要介绍一种方法：周期图法。 周期图法 由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。下面讨论离散随机信号序列的功率谱问题。 连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即：
随机信号分析题目及答 案 HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】
1. （10分）随机变量12,X X 彼此独立，且特征函数分别为12(),()v v φφ，求下列随机变量的特征函数： （1） 122X X X =+ （2）12536X X X =++ 解：（1） ()121222()jv X X jvX jv X jvX X v E e E e E e e φ+??????===??????? （2） () 121 2 536536 ()jv X X jv X jv X jv X v E e E e e e φ++????==????? ? 2. （10分）取值()1,1-+，概率[0.4,0.6]的独立 ()半随机二进制传输信号()X t ，时隙长度为T ，问： （1） 信号的均值函数()E X t ????； （2） 信号的自相关函数(),X R t t τ+； （3） 信号的一维概率密度函数();X f x t 。 解：（1）()10.410.60.2E X t =-?+?=???? (2) 当,t t τ+在同一个时隙时： 当,t t τ+不在同一个时隙时：
（3）()()();0.610.41X f x t x x δδ=-++ 3. （10分）随机信号0()sin()X t t ω=+Θ，()()0cos Y t t ω=+Θ，其中0 ω为常数，Θ为在[]-,ππ上均匀分布的随机变量。 （1） 试判断()X t 和()Y t 在同一时刻和不同时刻的独立性、相关性及正交性； （2） 试判断()X t 和()Y t 是否联合广义平稳。 解： （1） 由于X (t )和Y(t )包含同一随机变量 θ，因此非独立。 根据题意有 1 2f ()θπ = 。 []001 sin()02E[X(t )]E t sin(w t )d π π ωθθπ - =+Θ= +=?， 由于0XY XY R (t,t )C (t,t )==，X (t )和Y(t )在 同一时刻正交、线性无关。 除()012w t t k π-=±外的其他不同时刻 12120XY XY R (t ,t )C (t ,t )=≠，所以1X (t )和2Y(t )非正交 且线性相关。
Classical Power Spectrum Estimation Abstract With the increasing need of spectrum, various computational methods and algorithms have been proposed in the literature. Keeping these views and facts of spectrum shaping capability by FRFT based windows we have proposed a closed form solution for Bartlett window in fractional domain. This may be useful for analysis of different upcoming generations of mobile communication in a better way which are based on OFDM technique. Moreover, it is useful for real-time processing of non-stationary signals. As per our best knowledge the closed form solution mentioned in this paper have not been reported in the literature till date.This paper focuses on classical period spectral estimation and moderu spectral estimation based on Burg algorithm. By comparing various algorithms in computational complexity and resolution, Burg algorithm was used to signal processing finally. Experimental and simulation results indicated that digital signal processing system would meet system requirements for measurement accuracy. Keywords periodogram spectral estimation ; Burg algorithm I. INTRODUCTION When we expand the frequency response of any digital filter by means of Fourier series, we get impulse response of the digital filter in the form of coefficients of the Fourier series. But the resultant filter is unrealizable and also its impulse response in infinite in duration. If we directly truncate this series to a finite number of points we have to face with well known Gibbs phenomenon, so we modify the Fourier coefficients by
功率谱估计性能分析及其MATLAB实现 一、经典功率谱估计分类简介 1．间接法 根据维纳-辛钦定理，1958年Blackman和Turkey给出了这一方法的具体实现，即先由N个观察值，估计出自相关函数，求自相关函数傅里叶变换，以此变换结果作为对功率谱的估计。 2．直接法 直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N 个观察值直接进行傅里叶变换，得到，然后取其幅值的平方，再除以N，作为对功率谱的估计。 3．改进的周期图法 将N点的观察值分成L个数据段，每段的数据为M，然后计算L个数据段的周期图的平均 ，作为功率谱的估计，以此来改善用N点观察数据直接计算的周期图的方差特性。根据分段方法的不同，又可以分为Welch法和Bartlett法。 Welch法 所分的数据段可以互相重叠，选用的数据窗可以是任意窗。 Bartlett法 所分的数据段互不重叠，选用的数据窗是矩形窗。
二、经典功率谱估计的性能比较 1．仿真结果 为了比较经典功率谱估计的性能，本文采用的信号是高斯白噪声加两个正弦信号，采样率Fs=1000Hz，两个正弦信号的频率分别为f1=200Hz，f2=210Hz。所用数据长度N=400. 仿真结果如下： Figure1(a)示出了待估计信号的时域波形；
Figure2(b)示出了用该数据段直接求出的周期图，所用的数据窗为矩形窗； Figure2(c)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为矩形窗，长度M=128，数据没有加窗； Figure2(d)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为Hamming 窗，长度M=64，数据没有加窗； Figure2(e)是用Welch平均法求出的功率谱曲线，每段数据的长度为64点，重叠32点，使用的Hamming窗； Figure2(f)是用Welch平均法求出的功率谱曲线，每段数据的长度为100点，重叠48点，使用的Hamming窗； 2．性能比较 1)直接法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现 虚假谱峰； 2)间接法由于使用了平滑窗对直接法估计的功率谱进行了平滑，因此方差性能比直接 法好，功率谱比直接法估计的要平滑，但其分辨率比直接法低。 3)Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱 也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。 3．关于经典功率谱估计的总结 1)功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而 仍是目前较常用的谱估计方法。 2)谱的分辨率较低，它正比于2π/N，N是所使用的数据长度。 3)方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。 4)周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图 的方差性能，但往往又减小了分辨率和增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。 三、AR模型功率谱估计 1．A R模型功率谱估计简介 AR模型功率谱估计是现代谱估计中最常用的一种方法，这是因为AR模型参数的精确估计可以用解一组线性方程（Yule-Walker方程）的方法求得。其核心思想是：将信号看成是一个p 阶AR过程，通过建立Yule-Walker方程求解AR模型的参数，从而得到功率谱的估计。 由于已知的仅仅是长度有限的观测数据，因此AR模型参数的求得，通常是首先通过某种算法求得自相关函数的估计值，进而求得AR模型参数的估计值。常用的几种AR模型参数提取方法有： 1)自相关法 假定观测数据区间之外的数据为0，在均方误差意义下使得数据的前向预测误差最小。
经典功率谱估计 直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); window=boxcar(length(xn)); %矩形窗 nfft=1024; [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx)); 间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); 改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。 1. Bartlett法 Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。 Matlab代码示例： clear； Fs=1000; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
《现代信号处理》 姓名：李建强 学号：201512172087 专业：电子科学与技术 作业内容：在MATLAB平台上对一个特定的平稳随机信号进行经典功率谱估计和现代功率谱估计的比较 一、前言 功率谱估计是信息学科中的研究热点,在过去的30多年里取得了飞速的发展。在许多工程应用中，它能给出被分析对象的能量随频率的分布情况。平滑周期图是一种计算简单的经典方法，它的主要特点是与任何模型参数无关，但估计出来的功率谱很难与信号的真是功率谱相匹配。与周期图方法不同，现代谱估计主要是针对经典谱估计(周期图和自相关法)的分辨率低和方差性能不好的问题而提出的。其使用参数化的模型，能够给出比周期图方法高得多的频率分辨率。其内容极其丰富，涉及的学科和领域也相当广泛，按是否有参数大致可分为参数模型估计和非参数模型估计，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等;后者有最小方差方法、多分量的MUSIC方法等。 二、总体概述 本次实验分别使用经典的功率谱估计（如周期图法）与AR模型法对某一特定的平稳随机信号进行其功率谱估计，由图像得到信号的频率。利用MATLAB平台，直观形象地观察并比较二者估计效果的区别，以便于加深对功率谱估计的理解和掌握。 三、具体的实现步骤 1、经典法功率谱估计 周期图法又称直接法，它是从随机信号x(n)中截取N长的一段，把它视为能量有限的
真实功率谱的估计的一个抽样。 1.1、实现步骤 （1）、模拟系统输出参数x(n)=A*sin(2πf1*n)+B*sin(2πf2*n)，包括序列长度N（128或512或1024，加性高斯白噪声（AGWN）功率一定，设置A，B，f1，f2，n的值。 （2）、应用周期图法（不加窗）对信号的功率谱密度进行估计，使用直接法在MATLAB 平台上进行编程实现。 （3）、输出相应波形图，进行观察，记录。 1.2 MATLAB源代码实现 clear all; %清除工作空间所有之前的变量 close all; %关闭之前的所有的figure clc; %清除命令行之前所有的文字 n=1:1:128; %设定采样点n=1-128 f1=0.2; %设定f1频率的值0.2 f2=0.213; %设定f2频率的值0.213 A=1; %取定第一个正弦函数的振幅 B=1; %取定第一个正弦函数的振幅 a=0; %设定相位为0 x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a); %定义x1函数，不添加高斯白噪声x2=awgn(x1,3); %在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0; %定义临时值，并规定初始值为0 temp=fft(x2,128); %对x2做快速傅里叶变换 pw1=abs(temp).^2/128; %对temp做经典功率估计