序列相关
改革开放的历史意义-t3000
2023年3月18日发(作者：园林培训)第6章功率谱估计
离散随机序列的特征描述
平稳随机序列通过LTI系统
经典功率谱估计
现代功率谱估计
6.1离散随机序列的特征描述
随机过程的分布函数
随机信号的数字特征
平稳各态遍历随机信号的时域描述
平稳各态遍历随机信号的频域描述（功率谱密度）
一、随机过程的分布函数
{X[k],kZ}表示一个随机过程
一维分布函数
二维分布函数
N维分布函数
二、随机信号的数字特征
均值
方差
自相关函数
互相关函数
三、平稳各态遍历随机信号的时域描述
1平稳随机序列
指统计特性不随时间的平移而变化的那一类随机序列
严平稳随机序列：
宽平稳随机序列：
x
mkXE]}[{
][]}[][{nRnkXkXE
x
平稳随机信号自相关函数特性
（1）对称性
][][nRnR
xx
][][*nRnR
xx
),,;,,,(),,;,,,(
21212121
nknknkxxxFkkkxxxF
NNNN
)][(),(xkXPkxF
)][,][(),;,(
22112121
xkXxkXPkkxxF
),,;,,,(
2121NN
kkkxxxF
)][,][(
11NN
xkXxkXP
]}[{][kXEkm
x
][]}[{}])[][{(][2222kmkXEkmkXEk
xxx
]}[][{],[
2121
kXkXEkkR
x
]}[][{],[
2121
kYkXEkkR
xy
（2）极限值
0n
]}[{]0[2kXER
x
n2][
xx
mR
（3）不等式
][]0[nRR
xx
2.各态遍历随机信号
集平均等于时间平均
N
Nk
N
x
kx
N
kXEm][
12
1
lim]}[{
四、平稳各态遍历随机信号的频域描述
功率谱密度
维纳——辛钦公式
n
nj
xx
enRP)()(
dePnRnj
xx
)(
2
1
)(
当自相关函数绝对可积时，平稳随机信号的自相关函数和
功率谱密度是一对傅里叶变换对。
]),(
12
1
[lim)(
2
NF
N
EP
X
N
x
][][
12
1
lim]}[][{][nkxkx
N
nkXkXEnR
N
Nk
N
x
N
Nk
x
N
xx
mkx
N
mkXE222]][[
12
1
lim}]][{[
6.2平稳随机序列通过LTI离散时间系统
输出序列的均值
输出序列的自相关函数
输出序列的功率谱
输入/输出序列的互相关函数及互功率谱
平稳随机序列通过LTI系统
一、输出序列的均值
H(ej0)
二、输出序列的自相关函数
][][][nRnRnR
xhy
Ry[n]是系统单位脉冲响应h[k]的自相关函数Rh[n]与输入随机序列X[k]
的自相关函数Rx[n]的卷积.
系统单位脉冲响应h[k]是确定信号，其自相关函数定义为
]}[{DTFT]}[{DTFTnRnR
xh
]}[{][]}[{][nkxEnhkyEkm
n
y
n
x
nhm][
)(][0jeHmkm
xy
三、输出序列的功率谱
四、输入/输出序列的互相关函数及互功率谱互相关
][][]}[][{][nRnhnkXkYEnR
xyx
][*][]}[][{][nRnhknykxEnR
xxy
互功率谱
)()(]}[{DTFT)(
x
j
xyxy
PeHnRP
)()(*]}[{DTFT)(
x
j
yxyx
PeHnRP
[例]一离散时间平稳白噪声通过一阶IIR数字滤波器
求输出的自相关函数、平均功率和功率谱。
零均值白噪声的特征
0]}[{kXE2)(P
][d
π2
1
][2j2
π
π
nΩenRn
x
解：
az
az
zH
11
1
)(
][][kukhk
j
j
e
eH
1
1
)(
2
2
cos21
1
)(
jeH
（1）计算输出的自相关函数
]}[][{DTFT]}[{DTFT)(nRnRnRP
xhyy
)()()(
2
x
j
y
PeHP
]}[{DTFT]}[{DTFTnRnR
xh
2
)(]}*][{DTFTjeHnhnh
P
x
()
1]1[][][kykxky
][][][nkhkhnR
k
h
0
1
0
1
2
2
0
n
n
n
knk
nk
n
knk
k
21
n
][][][nRnRnR
xhy
][2nR
h
2
2
1
n
（2）输出平均功率
]0[][
12
1
lim2
y
N
Nk
N
Rky
N
2
2
1
（3）输出功率谱
)()()(
2
x
j
y
PeHP
2
2
cos21
6.3经典功率谱估计
谱估计的质量
相关法(间接法)
周期图法(直接法)
周期图法的改进
利用MATLAB实现功率谱估计
一、谱估计的质量
1．估计量的偏差
}
ˆ
{}
ˆ
{biaE
2．估计量的方差
})}
ˆ
{
ˆ
({}
ˆ
var{2EE
的一致估计为则称若
ˆ
,0}
ˆ
var{lim,0}
ˆ
{bia
N
3．估计量的均方差
}
ˆ
{bia}
ˆ
var{})
ˆ
({}
ˆ
{MSE2E
二、相关法(间接法)进行功率谱估计
相关法的理论基础
自相关函数估计的计算
相关法进行功率谱估计
功率谱估计的质量
1.维纳—辛钦定理
)(][
x
F
x
PnR
计算方法：
（1）由随机序列一个样本的N个观测值计算自相关函数的估计
（2）对][
ˆ
nR
x
进行DTFT即得该随机序列的功率谱估计)(
ˆ
x
P
2.自相关函数估计的计算
X[k]是宽平稳各态遍历随机信号，x[k]是其一个样本
N
Nk
N
x
nkxkx
N
nR][][
12
1
lim][
已知x[k]的N个观测值x[0],x[1],,x[N-1]，则自相关函数
的估计为
][][
1
][
ˆ
1
0
nkxkx
N
nR
N
k
x
][*][
1
nxnx
N
1)1(NnN
[例]已知平稳各态遍历的实随机序列X[k]的单一样本的N个观
测值为x[k]={1,0，-1}，试计算该随机序列的自相关函数估计。
解:
][*][
1
][
ˆ
nxnx
N
nR
x
}1,0,2,0,1{
3
1
利用MATLAB计算相关函数的估计
1.利用conv函数计算
2.利用数字处理工具箱中提供的函数xcorr
xcorr(x,y);%随机序列X和Y的互相关
xcorr(x);%随机序列X的自互相关
利用DFT计算自相关函数的估计
1)对x[k]补零形成L点序列
)12]([NLkx
L
2)]}[{DFT][kxmX
LL
101
1
0
1
10
1
1
1
000
0
3)}][{IDFT
1
][
ˆ2mX
N
nR
Lx
3.相关法进行功率谱估计
][
ˆ
),(
ˆ
DTFT,DFT
][
ˆ
][mPPnRkx
xxx
估计
nN
k
xnkxkx
N
nR
1
0
][][
1
][
1,][
ˆ
)(
ˆ
NLenRPjn
L
Ln
xx
1][][
1
][
1
0
Nnnkxkx
N
nR
nN
k
x
10Nn
][][
1
][
ˆ
1
0
nkxkx
N
nR
nN
k
x
0)1(nN
nN
l
N
nk
x
nlxlx
N
nkxkx
N
nR
1
0
1
][][
1
][][
1`
][
ˆ
[例]已知实平稳随机序列X[k]单一样本的N个观测值为
x[k]={1,0，-1}，
试利用相关法估计其功率谱。
解：X[k]的自相关函数估计值为
}1,0,2,0,1{
3
1
][*][
1
][
ˆ
nxnx
N
nR
x
0N1
x[k]
nN1n
x[k+n],
0n
N1n
x[k+n],0n
n
对][
ˆ
nR
x
进行傅里叶变换得X[k]的功率谱估计
]}[
ˆ
{DTFT)(nRP
x
x
}2{
3
1
22jjee)2cos1(
3
2
4.相关法功率谱估计的质量
功率谱估计的质量与自相关函数估计的
质量密切相关
][]}[{bianR
N
n
nR
xx
])[][][(
1
]}[
ˆ
var{2nrRnrRrR
N
nR
r
N，偏差、方差趋于零，是一致估计。
N固定时，nN，偏差、方差较大
三、周期图法(直接法)进行功率谱估计
周期图法功率谱估计的计算
周期图法功率谱估计的质量
1.周期图法功率谱估计的计算
已知:
其它0
1,,1,0][
][
Nkkx
kx
N
方法基础:
][*][
1
][
ˆ
nxnx
N
nR
NNx
由维纳—辛钦定理
]}[*][{DTFT
1
)(
ˆ
nxnx
N
P
NNx
)()(
1
j
N
j
N
eXeX
N
2
)(
1
j
N
eX
N
周期图法功率谱估计的步骤
2
DTFT)(
1
)()(][
j
N
x
j
NN
eX
N
PeXkx功率谱估计
2
DFT][
1
][][][mX
N
mPmXkx
N
x
NN
功率谱估计
其中
1
0
][]}[{DTFT)(
N
k
kj
NN
j
N
ekxkxeX
1
0
2
][]}[{DFT][
N
k
mk
N
j
NNN
ekxkxmX
[例]已知实平稳随机序列X[k]单一样本的N个观测值为
x[k]={1,0，-1}，试利用周期图法估计其功率谱。
解：对x[k]进行离散时间傅里叶变换
2
1
0
1][)(j
N
k
kjj
N
eekxeX
功率谱估计为：
)()(
1
)(
1
)(j*j
2
jeXeX
N
eX
N
I
NNNN
)1)(1(
3
1
22jjee
平稳高斯白噪声功率谱估计结果(周期图法)
2.周期图法功率谱估计的质量
均值nj
N
Nn
N
enR
N
nN
IE
][)}({
1
)1(
N，E{IN(W)}=Px(W)}，渐进无偏估计
方差}
sin
)sin(
1{)}(var{
2
4
N
N
I
N
N增加，方差不减小，不是一致估计
四、周期图法的改进
问题的提出
平滑周期图(Blackman-Tukey法)
平均周期图法(Bartlett法)
重叠平均周期图法(Welch法)
1.问题的提出
周期图法进行功率谱估计,方差不随N的增加减小.
如何提高谱估计质量?减小方差方法:
1.对自相关函数估计值加窗
2.将N个观测值分段,计算各段的周期图,再取平均
2.平滑周期图(Blackman-Tukey法)
对自相关函数估计值加窗,将误差较大的估计值截去
jn
x
N
Nn
M
enRnwP][
ˆ
][)(
1
)1(
窗函数w[n](M Mnnw Mnnwnw wnw 0][ ][][ 1]0[][0 B-T法进行功率谱估计的主要步骤 (1)利用观测数据估计自相关序列。 (2)对自相关函数估计值加窗。 (3)计算加窗后自相关函数的DTFT。 优点：PM(W)波动比IN(W)小，可证是一致估计 缺点：降低了频率分辨率 3.平均周期图法(Welch-Bartlett法) 将随机序列X[k]的N个观测值分成A段 1,1,0;1,1,0];[][MkAikiMxkxi 第i段序列的周期图为 2 )( 1 )(j M i M eX M I 平均周期图 )( 1 )( 1 0 i M A i A M I A P 平均周期图法估计质量 )}(var{ 1 )}(var{i M A M I A P A,方差为零，是一致估计 因为 1 )1( ][)]([bia N Nn nj N enR N n I 所以)]([bia)]([bia N i M II 平均周期图方差减小的代价之一是偏差增大 4.重叠平均周期图法(Welch法) 平均周期图法优点：减小方差 缺点：增加估计的偏差，降低了谱的分辨率 原因：分段即加窗，段越多，窗越短，主瓣宽度越大 解决方法：将各段数据有一定程度的重叠 平稳高斯白噪声功率谱估计结果(Welch法) 五、利用MATLAB进行非参数功率谱估计 周期图法 [Pxx,F]=PERIODOGRAM(X,WINDOW,NFFT,Fs) X：进行功率谱估计的输入有限长序列； WINDOW：指定窗函数，默认值为矩形窗(boxcar)； NFFT：DFT的点数，NFFT>X，默认值为256； Fs：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； F：Pxx值所对应的频率点 Welch-Bartlett平均周期图法的MATLAB实现 [Pxx,F]=PSD(X,NFFT,Fs,WINDOW,NOVERLAP) X,NFFT,Fs用法同periodogram函数； WINDOW：指定窗函数，默认值为hanning窗； NOVERLAP指定分段重叠的样本数。 如果使用boxcar窗，且NOVERLAP=0，则可得到 0123 -40 -20 0 20 M=256 Frequency P o w e r S p e c t r al ( d B) 0123 -40 -20 0 20 M=128 Frequency P o w e r S p e c t r al ( d B) 0123 -40 -20 0 20 M=64 Frequency P o w e r S p e c t r al ( d B) 0123 -40 -20 0 20 M=32 Frequency P o w e r S p e c t r al ( d B) 0N-1 i=1i=3 i=2 i=A MM M M Bartlett法的平均周期图。 如果NOVERLAP=L/2，则可得到重叠50%的Welch 法平均周期图。 [Pxx,F]=PWELCH(X,WINDOW,NOVERLAP,NFFT,Fs) [例]一序列含有白噪声和两个频率间隔很近的余弦 信号，设 ][)32.0cos()3.0cos(][kwkkkx 分别采用周期图法和Welch法估计该序列的功率谱。 %PowerSpectralEstimation:Periodogram N=512;Nfft=1024;Fs=2*pi; n=0:N-1; xn=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); XF=fft(xn,Nfft); Pxx=abs(XF).^2/length(n); index=0:round(Nfft/2-1); f=index*Fs/Nfft; plot(f,10*log(Pxx(index+1))),grid 或直接采用periodogram函数 %PowerSpectralDensityUsingWelchAlgorithm N=512;Nfft=1024;Fs=2*pi; n=0:N-1; xn=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); L=input('L=') window=boxcar(L); noverlap=L/2; [Pxx2f]=psd(xn,Nfft,Fs,window,noverlap); plot(f,10*log(Pxx2)),grid 周期图法进行功率谱估计的结果 00.511.522.53 -60 -40 -20 0 20 40 60 Frequency(Hz) P o w e r S p e c t r al D e n si t y ( d B/ H z ) PSDEstimation:Periodogram Welch法，M=64的谱估计结果 Welch法，M=128的谱估计结果 6.4现代谱估计简介 问题提出 平稳随机信号的参数模型 AR模型参数与自相关函数的关系 AR模型参数与线性预测滤波器的关系 Y-W方程的L-D递推算法 伯格(Burg)递推算法 利用MATLAB进行AR模型功率谱估计 一、问题提出 经典法存在问题: 1.方差性能不好，不是Px(W)的一致估计 2.平滑周期图和平均周期图改善了周期图的方差性能，但却降低了谱分 辨率和增大了偏差。 3.可能使短序列的功率谱估计出现错误的结果 出现问题的原因: 00.511.522.53 -60 -40 -20 0 20 40 60 Frequency P o w e r S p e c t r al ( d B) WelchPSDEstimation,M=64 00.511.522.53 -60 -40 -20 0 20 40 60 Frequency P o w e r S p e c t r al ( d B) WelchPSDEstimation,M=128 将观测数据以外的数据一律视为零，与实际不符。 参数模型法的基本思想 根据所研究信号的先验知识，对观测数据以外的数据作出某 种比较合理的假设。 方法： 选择一个好的模型，在输入是冲激函数或白噪声的情况下，使 其输出等于所研究的信号，至少也是对该信号的一个良好近似。 利用已知的自相关函数或数据求模型的参数。 )利用求出的模型参数或数据估计该信号的功率谱。 二、平稳随机信号的参数模型 AR模型 )( 1 1 1 )( 0 zA za zH p n n n MA模型 q l l l zbzH 1 1)( ARMA模型 )( )( 1 )( 0 0 zA zB za zb zH p n n n q l l l 若输入白噪声的功率谱 2)( e P 则输出序列的功率谱为 2 2 2 )()()()(j e j x eHPeHP 若能确定模型中各参数an和bl就可以求得功率谱Px() 三、AR模型参数与自相关函数的关系 ][][][ 1 kwnkyaky n p n Yule-Walker（Y-W）方程 0 0 0 1 ]0[]2[]1[][ ]2[]0[]1[]2[ ]1[]1[]0[]1[ ][]2[]1[]0[2 2 1 pyyyy yyyy yyyy yyyy a a a RpRpRpR pRRRR pRRRR pRRRR 若已知Ry[n]，由Y-W方程解出各参数a1,a2,…,ap，则可 由AR模型参数获得功率谱Py()的估计值。 四、AR模型参数与线性预测滤波器的关系 前向线性预测滤波器 y[k]的预测值][ ˆ ky由其过去值y[k-1],y[k-2],,y[k-p]的线性加权 得到。 前向预测误差 ][)(][][ ˆ ][][ 1 nkynakykykyke p p n f p 前向预测误差滤波器系统函数 后向线性预测滤波器 由y[k],y[k-1],,y[k-p+1]p个数据预测数据y[k-p] 后向预测误差 ][ ˆ ][][pkypkykeb p ][)(][ 1 pnkynapky p p n 前向预测误差滤波器系统函数 )(])(1[)(1 1 zAzznazzApn p p n pb 五、Y-W方程的L-D递推算法 (1)计算自相关函数的估计值 (2)由自相关函数的估计值，递推2 21 ,,,, pp aaa ]0[ ]1[ )1( 1 y y R R a 2 11 2 1 )1](0[aR y 2 1 1 1 1 ][)(][ )( p yp p n y p npRnapR pa )1,,2,1()()()()( 11 pnnpapanana pppp 2 1 22])(1[ ppp pa A(z)y[k] ][kef p Ab(z)y[kp] ][keb p (3)求出功率谱估计 2 1 2 1 )( ˆ p n nj n AR ea P [例]利用L-D算法对x[k]=cos(0.3p)+cos(0.32p)+[k]进行 谱估计 频率分辨率比N=128点高 六、伯格(Burg)递推算法 L-D算法缺点： 在计算相关函数估计时，对N个观测数据以外的数据作 零的假设，故谱估计误差较大。 伯格(Burg)递推算法基本思想： 直接从观测的数据利用线性预测器的前向和后向预测的 总均方误差之和为最小的准则来估计反射系数，进而通过 L-D算法的递推公式求出AR模型优化的参数。 伯格(Burg)递推算法步骤 (1)确定初始条件 ][][][ 00 kykekebf][ 11 0 2 0 ky N N k (2)从p=1开始迭代计算:计算AR模型参数 }]1[][{ ]1[][2 2 1 2 1 1 11 1 keke keke K b p f p N pk b p f p N pk p )1,,2,1()()()( 11 pnnpaKnana pppp 递推p阶均方误差2 1 22)1( ppp K (3)递推高一阶前、后向预测误差 ]1[][][ 11 keKkekeb pp f p f p ][]1[][ 11 keKkekef pp b p b p (4)若阶数小于p，则阶数加1，回到步骤(2)进行下一次 迭代，直到达到预定阶数p。 估计功率谱 2 1 2 1 )( ˆ p n nj n AR ea P 七、计算AR模型参数的MATLAB函数 A=LEVINSON(R,ORDER) ORDER:AR模型的阶数；R:观测序列的自相关函数； A:白噪声序列的方差和AR模型参数。 [A,E,K]=ARYULE(X,ORDER) X为观测序列；E为预测误差；K为反射系数。 [A,E,K]=ARBURG(X,ORDER) Pxx=PYULEAR(X,ORDER) Pxx=PBURG(X,ORDER) 直接绘制功率谱估计的曲线 PYULEAR(X,ORDER,NFFT,Fs) PBURG(X,ORDER,NFFT,Fs) NFFT:为DFT点数，默认值为256； Fs:绘制功率谱曲线的抽样频率，默认值为1。 [例]利用Burg法对x[k]=cos(0.3p)+cos(0.32p)+[k]进行谱估计 N=512;Nfft=1024;Fs=2*pi; n=0:N-1; xn=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); order=50;figure(1) pburg(xn,order,Nfft,Fs) title('BurgAlgorithm,p=50') order=80;figure(2) pburg(xn,order,Nfft,Fs) title('BurgAlgorithm,p=80') AR模型阶数p=50的谱估计结果(Burg法) AR模型阶数p=80的谱估计结果(Burg法) 00.511.522.53 -60 -40 -20 0 20 40 60 Frequency(Hz) P o w e r S p e c t r al D e n si t y ( d B/ H z ) PSDUsingBurgAlgorithm,p=50 00.511.522.53 -60 -40 -20 0 20 40 60 Frequency(Hz) P o w e r S p e c t r al D e n si t y ( d B/ H z ) PSDUsingBurgAlgorithm,p=80