_周期图法的功率谱估计
随机信号的功率谱估计方法
数字信号处理II——随机信号的功率谱估计方法一、实验目的1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。
2．观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。
3．学习使用FFT 提高谱估计的运算速度。
4．体会非参数化功率谱估计方法的优缺点。
二、实验原理与方法假设信号()x n 为平稳随机过程，其自相关序列定义为:{}*()()()m E x n x n m φ+@(0.1)其中{}E •表示取数学期望，{}*•表示取共轭。
根据定义，()x n 的功率谱密度()P w 与自相关序列()m φ存在如下关系：()()j mm P m eωωφ+∞-=-∞=∑ (0.2)1()()2j m m P e d πωπφωωπ-=⎰ (0.3)然而，实际中我们很难得到准确的自相关序列()m φ，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。
目前常用的线性谱估计方法有两种：自相关函数法和周期图方法，本实验将对这两种方法分别予以讨论。
1．自相关函数法假设已知随机信号()x n 的N 个观测样本，则其自相关序列可以用下式进行估计:||1*01ˆ()()()||1||N m n m x n x n m m N N m φ--==+≤--∑ (0.4)当仅使用长度为2M-1的自相关序列时，对其进行傅立叶变换即可得到功率谱估计如下:11ˆˆ()()M j m m M Pm e ωωφ--=-+=∑(0.5)其中M 为加窗长度，Re ()cMW m 为矩形窗函数，定义如下: Re 1,||()0,||cMm M Wm m M <⎧=⎨≥⎩(0.6)因此， ˆ()Pw 在一定程度上可以看作是“真正的功率谱()P w ”与窗函数傅立叶变换的卷积。
矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣，旁瓣效应使那些处于旁瓣附近功率较小的频率分量被淹没掉。
经典功率谱估计方法实现问题的研究
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 ）法这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。
功率谱估计方法的比较
功率谱估计方法的比较1.周期图法周期图法是最简单直观的功率谱估计方法之一，通过将信号分成多个长为N的区间，计算每个区间内信号的一维傅里叶变换，然后将这些变换结果平方并取平均得到功率谱。
该方法简单快速，但由于其需要使用多个区间的数据进行平均，因此对信号长度有较高的要求，且在信号存在非平稳性时，该方法不适用。
2.自相关法自相关法是一种经典的功率谱估计方法，通过计算信号的自相关函数来估计功率谱。
具体步骤是将信号与其自身的延迟序列进行点乘，并取平均得到自相关函数。
然后对自相关函数进行傅里叶变换，得到功率谱估计值。
该方法计算简单，但精度一般，且在信号长度较长时计算复杂度较高。
3.傅里叶变换法傅里叶变换法是一种经典的功率谱估计方法，通过对信号直接进行傅里叶变换得到功率谱。
该方法计算简单，精确度高，但对信号的长度存在要求，较长的信号长度能提供更高的分辨率。
此外，傅里叶变换法只适用于周期性信号。
4.平均周期图法平均周期图法是一种对周期图法的改进。
它将信号分为多段，并对每一段进行周期图计算，然后将计算结果平均得到平均周期图。
与周期图法相比，平均周期图法可以降低误差，提高估计精度。
然而，该方法仍然对信号长度有一定要求，并且计算复杂度较高。
5.移动平均法移动平均法是一种基于滑动窗口的功率谱估计方法，其基本思想是通过对信号进行多次滑动窗口处理，将窗口内信号的傅里叶变换结果平方并取平均得到功率谱估计值。
该方法在计算复杂度上较低，适用于非平稳信号的功率谱估计。
但是，由于窗口大小的选择存在权衡，需要根据实际情况进行合理设置。
总结起来，各种功率谱估计方法各有优劣。
周期图法和自相关法计算简单，但方法的精度较低，受信号长度限制且无法处理非平稳信号。
傅里叶变换法具有较高的计算精度，但对信号的长度和周期性要求较高。
平均周期图法和移动平均法对周期图法进行了改进，在精度上有所提高，但计算复杂度较高。
因此，在实际应用中，需要根据具体的信号特点和处理要求选取合适的功率谱估计方法。
第三章功率谱估计_1_
1 N −1 ˆ 样本的自相关函数: Rx (k ) = ∑ x(n) x* (n + k ) k = 0,1,..., M 1 N n =0 周期图间接法： 功率谱：Px (ω ) =
k =− M
M<N
∑
M
ˆ Rx (k )e− jωk
周期图的不足之处
当观测数据数量N趋于无穷时，估计方差仍不趋于零，是非一 致估计。 分辨率和窗长N有关，是低分辨率估计，无法完成功率谱的高 分辨。
周期图法
数据窗
有偏估计,平滑性差 加窗函数
Px (ω ) =
1 N
N −1 k =0
x ( n ) c ( n ) e − jnT ω ∑
n=0
N −1
2
谱窗
功率谱曲线平滑, 但分辨率下降
Px (ω ) = ∑ R x ( k ) w ( k ) e − jkT ω
要提高分辨率，使用参数化的谱估计！ 经典谱估计：使用FFT的谱估计 现代谱估计：参数化谱估计
周期图及其改进方法
由N个离散随机数据样本x(0), x(1),..., x( N − 1), 估计信号的功率谱。 频谱： X N (ω ) = ∑ x(n)e− jωn
n =0 N −1
周期图直接法： 功率谱：Px (ω ) = 1 1 2 X N (ω ) = N N x(n)e− jωn ∑
n =0 N −1 2
第3章 平稳过程的线性模型
3.2 平稳随机信号通过线性系统
y (n) = x(n) ∗ h(n) =
m = −∞
∑
∞
x(m )h(n − m )
如果x(n)为确定性信号
Y (e ) = X (e ) H (e )
北京交通大学(数字信号处理研究生课程)ch6_2经典功率谱估计
j2
)
2 3
1 cos( 2 )
周期图法(直接法)进行功率谱估计
2. 周期图法功率谱估计的质量
N ( N 1)
Nn N
R[n]e jn
N，E{IN()}= Px()，渐进无偏估计 方差
sin( N) var{I N ()} {1 } N sin
功率谱估计
周期图法功率谱估计的步骤
x N [k ] X N (e
L点DFT
DTFT
j
2 1 j ) P x ( ) X N (e ) N 功率谱估计
x N [ k ] X N [ m ] P x [ m ]
功率谱估计
1 N
X N [m]
N 1 n 1
ˆ ( ) P x
L
jn ˆ ˆ R[0] 2 Re{ R[n]e }
功率谱估计
例：相关法计算功率谱估计 已知实平稳随机序列X[k]单一样本的N个观测值为 x[k]={1, 0，1}，试利用相关法估计其功率谱。 解：X[k]的自相关函数估计值为 1 ˆ [ n ] {1, 0, 2, 0, 1} R x 3
解:
1 ˆ R x [ n] N
1 3
N 1 k 0
x[k ]x[k n]
x[ n ] x[ n 2]
x[ n ] {1, 0, 1}
x[ n 2] {1, 0, 1,}
ˆ [ n ] {1, 0, 2, 0, 1} R x 3
1
利用MATLAB计算自相关函数的估计
N 1 Rx [n] lim x[k ]x[k n] N 2 N 1 k N
[matlab实现经典功率谱估计]matlab功率谱估计
![[matlab实现经典功率谱估计]matlab功率谱估计](https://img.taocdn.com/s3/m/acd3e43931b765ce04081407.png)
[matlab实现经典功率谱估计]matlab功率谱估计1、直接法：直接法又称周期图法，它是把随机序列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));2、间接法：间接法先由序列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);3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.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)); nfft=1024;window=boxcar(length(n)); %矩形窗noverlap=0; %数据无重叠p=0.9; %置信概率[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot_Pxxc=10*log10(Pxxc(index+1));figure(1)plot(k,plot_Pxx);pause;figure(2)plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);3.2、Welch法Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。
窗口效应
一种改进的窗处理是韦尔奇（Welch）法，因为窗在两 边渐变为零，所以这种方法降低了由于重叠导致的段间 统计相依的效应，而且，一个合适的非矩形窗可以减小 “旁瓣效应”，即“谱泄漏”。 经常采用的窗函数有矩形窗、汉宁窗(Hanning)、海明窗 (Hamming)、凯瑟窗(Kaiser)和切比雪夫窗(Chebwin)等， 其中矩形窗、汉宁窗以及海明窗都是广义余弦窗的特殊 形式，它们可以看作是频率为0、 ／(N一1)和 ／(N一1)的余弦曲线的线性组合，这里N为窗的长度，
巴特利特平均周期图法的实现过程就是把数据分成互 不重叠的L 段，每段有M个样本，即
其中 令
, 是长度为M的矩形窗。 ,则第i段的周期图可有下面公式求得：
由上式可求得整个序列的谱估计：
巴特利特平均周期图法的仿真
在巴特利特平均周期图 法中，对于固定长度的 N，分段K数目愈大， 方差愈小，但M愈小， 会使偏倚增大，谱分辨 率变得愈差。 因此在应用巴特利特平 均周期图法时，需要在 偏倚和估计方差之间进 行权衡，根据实际要求 选择适当的M和N。
在MALAB中可以采用下面的命令来生成：
加窗后的窗口效应
加矩形窗 0
-5
-10
-15
-20
-25
-30
-35
0
50
100
150
200
250
300
矩形窗处理的谱估 计的主瓣较窄，分 辨率最好，同时其 主瓣附近的衰减比 其他两个窗函数小， 旁瓣部分泄漏比较 严重，所以其方差 特性最差，噪声水 平较高。
例如测量震动物体的自震频率时，则可以选用主瓣宽度比 较窄的Rectangular窗；当要求谱泄漏小时，则应选用旁 瓣幅度较小的窗。
功率谱估计
W(n)为零均值方差为1的AWGN，n=1,2,3……,128
1.1周期图法：
我们知道随机信号的功率谱和自相关函数是一对傅式变换对：
而自相关函数定义为：
对于平稳随机过程，并由功率谱的偶函数特性得：
实际得到的随机信号只能是它的一个样本的片断，因此只能用有限长的样本序列来估计功率谱，这相当于用一个有限宽度（N）的窗函数 去乘样本序列，于是有（用离散频率K代替ω）：
title('周期图法');
xlabel('Hz');
ylabel('dB/Hz');
window1=hamming(128);
noverlap=20; %数据20%的重叠
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,'onesided');
plot_Pxx1=10*log10(Pxx1);
仿真结果：
2.现代功率谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大熵谱分析法（AR模型法）、Pisarenko谐波分解法、Prony提取计点法、Prony谱分解法以及Carpon最大似然法。其中AR模型应用较多，具有代表性。常用的模型有ARMA模型、AR模型、MA模型。
这就是用样本序列片断的DFT来估计功率谱的式子。由于加了矩形窗，使得这种直接的周期图估计平滑性、一致性和分辨率不能满足实际要求，因此有必要对上式作一些修改，这些修改主要有两种方法：
1.分段平均：即将长度为N的数据分成L段（允许有重叠），分别求出每一段的功率谱，然后即以平均。这样L个平均的方插笔每个随机变量的单独方差小L倍。
自功率谱估计的经典方法
5．自功率谱估计的经典方法 1） 周期图法（直接法）对于时间序列)(n x N ，其傅里叶变换（DTFT ——离散时间信号的傅里叶变换）为∑-=-=1)()(N n nj N j N en x e X ωω，⎰-=ππωωωπd e e X n x n j j N N )(21)(记为)()(ωj N D TFTN e X n x −−→←)(n x N 的离散傅里叶变换（DFT ）为∑-=-=102)()(N n kn Nj N N en x k X π，∑-==12)(1)(N k kn Nj NN e k X N n x π记为)()(k X n x N D FTN −−→←若)(n x N 是信号)(n x 在时间域截断的结果，即)()()(n d n x n x N N ⋅= (5-58)其中，)(n d N 是单边矩形窗，其表达式为⎩⎨⎧-≤≤=其它,010,1)(N n n d N 而)(n x 是确定性功率信号（或随机信号的一个样本序列），则根据第三章的讨论结果知，=)(ωj x e S 2,)(1)(limlim ωωj N N j x N N e X Ne P ∞→∞→= (5-59) 反映了信号)(n x 的平均功率在频域的分布情况，称为平均功率谱密度。
因此，估计量2,,)(1)()(ˆωωωj N j x N j PER x e X Ne P e S == (5-60) 为信号)(n x 的功率谱的一个估计。
此估计方法称为直接法或周期图法。
在)(ˆ,ωj PER x eS 的实际运算中采用DFT ，ω在单位园上均匀取值。
当取Nπω2=∆时，（5-60）改写为2,,)(1)()(ˆk X Nk P k S Nx N N PER x ==，1,,1,0-=N k (5-61) 其中，∑-=-=12)()(N n nk NjN N en x k X π，1,,1,0-=N k当取N22πω=∆时，需对)(n x N 补N 个零后再作DFT ，此时（5-60）改写为 22,22,)(1)()(ˆk X Nk P k S Nx N N PER x ==，12,,1,0-=N k (5-62) 其中，)(2k X N 参见(5-42)、(5-33)式。
周期图估计法
一种信号功率谱密度估计方法。
它的特点是：为得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度N，即(1)(2)由于序列x(n)的离散傅里叶变换X()具有周期性，因而这种功率谱也具有周期性，常称为周期图。
早期的统计学者曾利用这种方法从大量的数据中寻找隐藏的周期性的规律。
周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。
因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。
由于随机起伏大，使用周期图不能得到比较稳定的估值。
一些学者对此作了改进。
为了减小随机起伏，M.S.巴特利特提出平均周期图法，即先把信号序列分为若干段，对每段分别计算其周期图，然后取各个周期图的平均作为功率谱的估值。
平均周期图可以减小随机起伏，但是，如果信号序列不是足够长，由于每段序列长度变短，功率谱估值对不同频率成分的分辨能力也随之下降。
另一种改进方法是将周期图与一个适当的频域窗函数相褶积，从而对周期图产生平滑作用，以减小随机起伏。
加窗处理的结果虽然可以使随机起伏减小，但也会使周期图的分辨能力下降。
P.O.韦尔奇提出一种把加窗处理与平均处理结合起来的方法。
先把分段的数据乘以窗函数(进行加窗处理)，分别计算其周期图，然后进行平均。
韦尔奇方法是较常用的一种计算方法。
为了得到较好的功率谱估值，加窗和平均处理均应兼顾减小随机起伏和保证有足够的谱分辨率两个方面。
周期图法的优点是能应用离散傅里叶变换的快速算法来进行估值。
对利用式(1)、(2)得到的功率谱估值进行傅里叶反变换，可以得到信号的自相关函数估值。
这种方法适用于长信号序列的情况，在有足够的序列长度时，应用改进的周期图法，可以得到较好的功率谱估值，因而应用很广。
用matlab做经典功率谱估计
用matlab做经典功率谱估计经典功率谱估计1、直接法：直接法又称周期图法，它是把随机序列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));2、间接法：间接法先由序列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);3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N 太小，谱的分辨率又不好，因此需要改进。
3.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));nfft=1024;window=boxcar(length(n)); %矩形窗noverlap=0; %数据无重叠p=0.9; %置信概率[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot_Pxxc=10*log10(Pxxc(index+1));figure(1)plot(k,plot_Pxx);pause;figure(2)plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);3.2、Welch法Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。
功率谱图应用
1.基本方法周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。
假定有限长随机信号序列为x(n)。
它的Fourier变换和功率谱密度估计存在下面的关系：式中，N为随机信号序列x(n)的长度。
在离散的频率点f=kΔf,有：其中，FFT[x(n)]为对序列x(n)的Fourier变换，由于FFT[x(n)]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。
下面用例子说明如何采用这种方法进行功率谱用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。
为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch 法）等方法加以改进。
2. 分段平均周期图法（Bartlett法）将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。
对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。
平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。
对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。
这两种方法都称为平均周期图法，一般后者比前者好。
程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。
与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。
3.加窗平均周期图法加窗平均周期图法是对分段平均周期图法的改进。
在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。
由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。
周期图法功率谱估计------窗口效应
nfft=512;
•plot(f,plot_Pxx);
window=boxcar(100);%矩形窗
•title('加矩形窗');
window1=hamming(100);%海明窗
•figure(2)
window2=blackman(100);%blackman 窗 noverlap=20; %数据无重叠
周期图法功率谱估计------窗口效应
主要内容
概述 周期图法 窗函数法 谱估计法比较 总结
数字信号处理的两个主要分支：
数字滤波 频谱分析
对随机信号的频谱分析——功率谱估计
对确定信号，可以用傅立叶变换；而随机 信号无始无终具无限能量，不满足傅立叶变
换绝对可积的条件.
Hale Waihona Puke 概述 功率谱：随机信号的功率谱反映的是随机 信号的频率成分及各成分的相对强弱。
);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,ran ge);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,ran ge);
加矩形窗后的窗口效应
加矩形窗 0
-5
-10
-15
-20
-25
-30
-35
0
50
Matlab 代码:
Fs=600; n=0:1/Fs:1;
•plot_Pxx=10*log10(Pxx); •plot_Pxx1=10*log10(Pxx1);
xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size (n));
对功率谱估计常用方法的探讨及应用分析
对功率谱估计常用方法的探讨及应用分析功率谱估计是信号处理中常用的一种方法，它可以将信号的频率特性展示出来，对于信号的分析和处理具有重要意义。
常用的功率谱估计方法包括周期图法、解析法、Welch方法、Bartlett方法和Burg方法等。
本文将对这些方法进行探讨并分析其应用。
周期图法是一种基本的功率谱估计方法，它基于傅里叶变换的思想，通过将信号分解为不同频率的正弦波分量，然后计算每个分量的功率，从而得到信号的频谱特性。
该方法的优点是计算简单，但对于非平稳信号或信号中存在窗函数时会引入谱漏，导致估计结果不准确。
解析法是一种使用解析信号估计功率谱的方法。
解析信号是通过原始信号与希尔伯特变换得到的，它具有正频谱和负频谱的特点。
该方法的优点是可以避免频谱漏失的问题，但计算量较大。
应用方面，解析法常用于振动信号的分析和故障诊断中。
Welch方法是一种常用的频谱估计方法，它通过对信号进行分段处理，然后对每个片段进行傅里叶变换，最后将各个片段的功率谱进行平均得到最终的估计结果。
这样做的好处是可以减小谱漏的影响，并且可以根据需要进行频谱分辨率和频率平滑的调整。
Welch方法在信号处理中应用广泛，如语音和音频处理、通信系统等。
Bartlett方法是Welch方法的特例，它将信号分成互不重叠的窗函数片段，然后进行傅里叶变换并对功率谱进行平均。
这种方法的优点是计算简单，但对于非平稳信号可能会引入谱漏现象，导致估计结果不准确。
Bartlett方法在多传感器信号处理和谱估计的实时应用中常用。
Burg方法是一种利用自回归(AR)模型估计功率谱的方法。
AR模型假设信号的当前值与过去若干个值相关，通过建立AR模型并对其参数进行估计，可以得到信号的频谱特性。
该方法的优点是可以很好地处理非平稳信号，并且对信号中的噪声具有较好的抑制效果。
Burg方法在信号处理中广泛应用于信号的谱分析和预测等领域。
综上所述，功率谱估计方法在信号处理中具有重要的应用价值。
功率谱估计的经典方法
功率谱估计的经典方法周期图法是最早被提出的功率谱估计方法之一、它基于信号的周期性，将信号分解成一系列频率分量，然后计算每个频率分量的功率谱密度。
周期图法主要分为周期自相关法和周期平均法两种。
周期自相关法通过计算信号的自相关函数，然后进行傅里叶变换得到功率谱估计结果。
周期平均法则是通过对多个信号周期进行平均得到功率谱估计结果。
平均法是功率谱估计的另一种常用方法。
它通过对信号进行多次采样，然后计算采样信号的傅里叶变换得到频谱，再对多个频谱进行平均得到功率谱估计结果。
平均法的优点是抗噪声能力强，可以提高功率谱估计的准确性。
自相关法是一种基于信号自身特性的功率谱估计方法。
它通过计算信号的自相关函数，然后进行傅里叶变换得到功率谱估计结果。
自相关法的优点是计算简单，但是对信号的平稳性要求较高。
递归方法是一种实时性较好的功率谱估计方法。
它通过对信号进行递推计算，每次计算结果作为下一次计算的输入，以此来估计信号的功率谱。
递归方法通常会使用窗函数来平滑信号，减小频谱分辨率。
递归方法的优点是计算效率高，可以用于实时信号处理。
除了这些经典方法，还有一些其他的功率谱估计方法，如Yule-Walker方法、Burg方法、最大熵方法等。
每种方法都有其适用的场景和特点，选择合适的方法需要根据具体需求和信号特性进行判断。
在实际应用中，功率谱估计可以用于信号处理、通信系统设计、频谱分析等领域。
它可以帮助我们了解信号的频谱分布特性，对信号进行分析和处理，从而实现更好的信号传输和处理效果。
无论是音频信号、图像信号还是通信信号，功率谱估计都具有重要的意义。
因此，掌握功率谱估计的经典方法是进行信号处理和频谱分析的基础。
计算功率谱
计算功率谱clear;fs=500; %采样率df=1; %频率分辨率N=floor(fs/df)+1; %计算的序列点数t=0:1/fs:(N-1)/fs; %截取信号的时间段f=0:df:fs; %功率谱估计的频率分辨率和范围xt=sin(2*pi*50*5)+2*sin(2*pi*130*t)+randn(1,length(t));%截取时间段上的离散信号样本序列%%%%利用周期图法进行功率谱估计，但是其得出的功率谱很不平滑，相应的估计协方差比较大%增加采样点数也不能使周期图变得更加平滑，这是周期图法的缺点，在后面对其改进。
Px=abs(fft(xt)).^2/(N^2); %功率谱估计Pav_tm=sum(xt.^2)/N; %在时域计算信号功率Pav_fn=sum(Px); %通过功率谱计算信号功率figure(1)subplot(221), %作出功率谱密度图plot(f,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱（dB）');title('周期图法得出的功率谱估计');%%%%对周期图法进行改进的思想是将信号分段进行估计，然后再将这些估计结果进行平均，%从而减小估计的协方差，使估计功率谱图变得平滑。
%本程序是将以上501点的信号分为3段，分别作周期图法估计，然后平均。
Px=(abs(fft(xt(1:167))).^2+abs(fft(xt(168:334))).^2+...abs(fft(xt(335:501))).^2)/3/((N/3)^2);%分为3段，每段采样点数为N的1/3Pav_tm=sum(xt.^2)/N; %在时域计算信号功率Pav_fn=sum(Px); %通过功率谱计算信号功率subplot(222), %作出功率谱密度图plot(0:3:fs,10*log10(Px));xlabel('频率（Hz）');ylabel('功率谱（dB）');title('采用分段估计平均的方法降低估计协方差');%%%%增加分段数可以进一步降低估计的协方差，然而每段中的数据点太少，就会使估计的%频率分辨率下降很多。
周期图法估计功率谱
随机信号谱估计方法的Matlab实现摘要：功率谱估计是随机信号分析中的一个重要内容。
从介绍功率谱的估计原理入手分析经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。
经典功率谱估计的方差大、谱分辨率差，分辨率反比于有效信号的长度，但现代谱估计的分辨率不受此限制。
给出了功率谱估计的应用。
关键词：功率谱估计；周期图法；AR参数法；1 引言在一般工程实际中，随机信号通常是无限长的，例如，传感器的温漂，不可能得到无限长时间的无限个观察结果来获得完全准确的温漂情况，即随机信号总体的情况，一般只能在有限的时间内得到有限个结果，即有限个样本，根据经验来近似地估计总体的分布。
有时，甚至不需要知道随机信号总体地分布，而只需要知道其数字特征，如均值、方差、均方值、相关函数、功率谱的比较精确的情况即估计值。
功率谱估计(PSD)是用有限长的数据估计信号的功率谱，它对于认识一个随机信号或其他应用方面都是重要的，是数字信号处理的重要研究内容之一。
功率谱估计可以分为经典谱估计(非参数估计)和现代谱估计(参数估计)。
2 .平均周期图法和平滑平均周期图法对于周期图的功率谱估计, 当数据长度N 太大时, 谱曲线起伏加剧, 若N 太小, 谱的分辨率又不好,因此需要改进。
两种改进的估计法是平均周期图法和平滑平均周期图法。
(1)Bartlett 法:Bartlett 平均周期图的方法是将N 点的有限长序列x(n)分段求周期图再平均。
Matlab 代码示例1:fs=600;n=0:1/fs:1;xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n));nfft=512;window=hamming(nfft); %矩形窗noverlap=0;%数据无重叠p=0.9;%置信概率[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p);index=0:round(nfft/2- 1);k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot_Pxxc=10*log10(Pxxc(index+1));figure(1)plot(k,plot_Pxx);figure(2)plot(k,[plot_Pxx plot_Pxx- plot_Pxxcplot_Pxx+plot_Pxxc]);matlab调试图下图(2)Welch 法:Welch 法对Bartlett 法进行了两方面的修正, 一是选择适当的窗函数w(n), 并在周期图计算前直接加进去, 加窗的优点是无论什么样的窗函数均可使谱估计非负。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
nfft=1024； %傅里叶变换点数 Xk=fft（x，nfft） ； Pxx=abs（Xk）.^2/N； index=0：round（nfft/2-1） ； k=index*N/nfft；% 频率（Hz） plot_Pxx=10*log10（Pxx（index+1） ） ；
%功率谱 dB
2
基本的周期图法
论
著
与创新研发的投入，提高科技与创新投入的效益，优先保证 关键技术研发创新和重大产业化项目建设等的要素需求， 促 进经济发展从外延扩张型向内涵提高型转变。同时，加强企 业监管和企业服务体系建设，严格限制处于产业链低端环 节、特色优势不明显、高消耗、低效益的行业或企业扩张产 能，支持优势企业做大做强，提升核心竞争力。
图3
Welch 法周期图法得到的功率谱
Matlab 程序示例： %信号 t = 0:0.001:0.6； y = sin（2*pi*50*t）+sin（2*pi*120*t） ； x = y+randn（size（t） ） ； N=length（x） ； %观察数据点数 n=1:1:N； %参数设置 Fs=1000； %每百年采样 100 次 nfft=512； %傅里叶变换点数 %window=hamming（25） ； window=boxcar（100） ；
论
著
noverlap=20 range='half'； %计算序列的 psd [Pxx，f]=pwelch（x，window，noverlap，nfft，Fs，range） ； plot_Pxx=10*log10（Pxx） ； %功率谱 dB plot（f，plot_Pxx） xlabel（'Hz'） ； ylabel（'功率谱 dB'） ； title（'用 welch 法计算的功率谱'） ； grid on；
x ( n )e
n 0
N 1
2 jn
i 1 n 0
L
M
2 i xm ( n ) e j n
。
图 2 给出了平均周期图法得到的功率谱。
图 1 给出了基本周期图法得到的功率谱。
图2 图1 基本的周期图法得到的功率谱
平均周期图法得到的功率谱
Matlab 程序示例： t = 0:0.001:0.6； y = sin（2*pi*50*t）+sin（2*pi*120*t） ； x = y + randn（size（t） ） ； N=length（x） ； %观察数据点数
5
多窗口的周期图法
由于普通的周期图功率谱估计只利用单一窗口， 因此在 序列始端和末端均会丢失相关信息。 多窗口法利用多个正交 窗口获得各自独立的近似功率谱估计，然后综合这些估计， 最终得到平稳信号序列的功率谱估计。 图 4 给出了多窗口周期图法得到的功率谱。
4
Welch 法
P.O.韦尔奇提出一种把加窗处理与平均处理结合起来的 方法。先把分段的数据乘以窗函数（进行加窗处理） ，分别 计算其周期图，然后进行平均，得到的功率谱值为：
中国西部科技
2013 年 09 月第 12 卷第 09 期总第 290 期
1
周期图法的功率谱估计
黄 英
（苏州经贸职业技术学院，江苏 苏州 215009） 摘 要：功率谱估计是数字信号处理的重要研究内容之一。对目前用得较多且最具代表性的周期图法进行研究，并针对它谱 分辨率比较低等缺点，利用 Matlab 实现了几种改进的周期图法。 关键词：谱估计；周期图；Matlab DOI：10.3969/j.issn.1671-6396.2013.09.001
Matlab 程序示例： %信号 t = 0:0.001:0.6； y = sin（2*pi*50*t）+sin（2*pi*120*t） ； x = y+randn（size（t） ） ； N=length（x） ； %观察数据点数 n=1:1:N；
2 %参数设置 Fs=1000； %采样频率 nfft=512； %傅里叶变换点数 window=boxcar（30） ； noverlap=0； p=0.9； %计算序列的 psd [Pxx]=psd（x，nfft，Fs，window，noverlap，p） ； index=0：round（nfft/2-1） ； k=index*Fs/nfft； plot_Pxx=10*log10（Pxx（index+1） ） ； %功率谱 dB plot（k，plot_Pxx） xlabel（'Hz'） ； ylabel（'功率谱 dB'） ； title（'用平均周期图法计算的功率谱'） ； grid on；
[作者简介] 黄英（1965—），女，副教授，研究方向为嵌入 式系统应用开发。
Power spectrum Estimation Based on Periodogram HUANG Ying
(Suzhou Institute of Trade & Commerce,Suzhou,Jiangsu 215009,China)
Matlab 程序示例： %信号 t = 0:0.001:0.6； y = sin（2*pi*50*t）+sin（2*pi*120*t） ； x = y+randn（size（t） ） ； N=length（x） ； %观察数据点数 n=1:1:N； %参数设置 Fs=1000； %每百年采样 100 次 nfft=512； %傅里叶变换点数 %window=hamming（25） ； window=boxcar（100） ； noverlap=20 range='half'； %计算序列的 psd nw=3.5；p=0.99； [Pxx，Pxxc，f]=pmtm（x，nw，nfft，Fs，p） ； plot_Pxx=10*log10（Pxx） ； %功率谱 dB plot（f，plot_Pxx） xlabel（'Hz'） ； ylabel（'功率谱 dB'） ； title（'用多窗口法计算的功率谱'） ； grid on； （下转第 8 页）
Abstract:Spectral estimation is an important part of digital signal processing.This paper studied the basic period diagram method,and realized some improved methods on Matlab. Keywords:Spectral Estimation;Periodogram;Matlab
plot（k，plot_Pxx） xlabel（'Hz'） ； ylabel（'功率谱 dB'） ； ； title（'用基本周期图法（直接计算）计算功率谱'） grid on； 基本的周期图法可以提高计算效率， 不需要计算自相关 函数，但谱分辨率较低。
周期图的基本原理是对观测到的数据直接进行傅立叶 变换，然后取模的平方就是功率谱。取平稳随机信号 x（n） 的有限个观察点 x（0） 、x（1） 、…x（n-1） ，则傅立叶变换 为： X N (e j )
参考文献 [1] 李小建.经济地理学[M].北京:高等教育出版社,2005 [2] [法]弗朗索瓦·佩鲁.增长极概念[J].原载经济学译丛,1988,(9) [3] 鄂冰.中心城市产业结构优化与升级理论研究[J].城市发展研究,2012,(4) [4] 未江涛.城市中心区产业结构调整优化的基本思路与对策研究[J].天 津行政学院学报,2009,(29) [5] 隋广军 . 中心城市经济结构转型与区域产业结构优化升级 [J]. 改革 , 2008,(1) [6] 唐志鹏.能源约束视角下北京市产业结构的优化模拟与演进分析[J]. 资源科学,2012,(1) [7] 朱群.重庆市产业结构调整及其优化思路研究[J].中国集体经济,2011,(36) [8] 贺亚锋.大连市产业结构优化升级研究[J].云南地理环境研究,2012,(1) [9] 周振华.论城市能级水平与现代服务业[J].社会科学,2005,(9)
1
引言
功率谱估计就是用已观测到的一组数据去估计一个平 稳信号的功率谱， 目的是根据有限数据在频域内提取被淹没 在噪声中的有用信号， 它对于认识一个随机信号或其他应用 方面来讲都是重要的，是数字信号处理的重要研究内容之 一， 在雷达系统追踪目标速度信息、 声呐系统寻找目标信息、 语音识别中去除噪声及进一步声学处理、生物学、天文学等 众多领域有着重要的应用。 功率谱估计可以分为经典谱估计 （非参数估计）和现代谱估计（参数估计） 。经典谱估计中 的周期图法是用得较多且最具代表性的方法。 本文首先研究 了基本的周期图法，接着针对它谱分辨率比较低等缺点，利 用 Matlab 实现了几种改进的周期图法。
通过实验可以看出，基本的周期图法谱分辨率较低，会 有很大的失真。平均周期图法和 Welch 法可以减小随机起 伏，收敛性较好，曲线平滑，估计的结果方差较小，但是功 率谱主瓣较宽， 分辨率低。 但是， 如果信号序列不是足够长， 由于每段序列长度变短， 功率谱估值对不同频率成分的分辨 能力也随之下降。 多窗口的周期图法功率谱比平均周期图法 和 Welch 法的功率谱谱峰窄一些，谱峰高度也有所增大，因 此多窗口的周期图法的功率谱表现得更加准确。 普通周期图 法的功率谱估计利用单一窗口， 所以在序列的始端和末端均 会丢失相关信息， 而多窗口的周期图法以增加窗口来充分利 用这些相关信息。
[作者简介] 黄吉秀， 女， 副教授， 四川财经职业学院副院长， 研究方向为经济学。
（上接第 2 页）
6
结论
参考文献 [1] 丁玉美.数字信号处理——时域离散随机信号处理[M].西安:西安电 子科技大学出版社,2002 [2] 陈怀琛,吴大正.MATLAB 及在电子信息课程中的应用 （第二版） [M]. 北京:电子工业出版社,2004 [3] 李富强,万红.MATLAB6 的语谱图显示与分析[J].微计算机信息,2005, 10(3):1-5 [4] 魏鑫,张平.周期图法功率谱估计中的窗函数分析[J].现代电子技术, 2005,(3)