随机信号的功率谱密度估计方法研究
【摘要】对确定性信号进行傅里叶变换是在频率域分析研究的理论基础，但对于随机信号来说，其傅里叶变换并不存在，因此转型研究它的功率谱。这样就可以用功率谱来描述随机信号的频域特性。功率谱是功率谱密度的简称，是自相关函数的傅里叶变换。随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。对功率谱密度的估计又称功率谱估计。功率谱估计是数字信号处理中的一个重要组成部分，其理论研究已经非常成熟，但是传统的仿真实现由于计算量的庞大而往往令人望而却步。采用功能强大的MA TLAB软件对功率谱估计进行仿真，便于深刻理解各种方法的特点，从而在实际工作中做出合理的选择。本文简要介绍了MA TLAB仿真软件的基本功能、特点和优势，然后具体介绍了周期图法。
【关键词】自相关函数;功率谱;MA TLAB;随机信号; 周期图法.
Random signal power spectral density estimation
method research
Abstract:Fourier transform of the deterministic signals is the theoretical foundation of the frequency domain analysis, but for random signal, its Fourier transform does not exist, so we change to study the power spectrum. So it can be used to describe the random power spectrum of the signal frequency characteristics 。Power spectrum is the power spectral density abbreviation, is the autocorrelation function of Fourier transformation 。Random signal power spectral density of the energy used to describe the signal characteristics on frequency changes relationship。To power spectral density estimate and says power spectrum estimation 。Power spectrum estimation is an important part of the digital signal processing whose theory is well established, as the large of calculated quantity of traditional simulation often prevent people to walk more. we emulate several linear power spectrum by using MA TLAB software, we could get characteristics of those power spectrum estimation and discussions of these types of power spectrum estimation, we compare and analysis them to facilitate a deep understanding of the characteristics of various methods , and make a reasonable choice .in practical work .This paper introduces the basic functions, features and advantages of MATLAB simulation software, then describe the classics power spectrum estimation principle such as cycle graph。
Key words:Autocorrelation function power spectrum estimation MATLAB random signal cycle graph
目录
引言 (1)
1绪论 (1)
1.1功率谱研究的发展过程 (1)
1.2功率谱估计方法提出 (2)
1.3功率谱估计应用及用途 (2)
1.4本课题研究任务 (2)
2 随机信号简介 (2)
2.1随机变量 (2)
2.1.1 随机变量的分布函数 (3)
2.1.2 数学期望、方差、标准差 (3)
2.1.3 随机向量 (3)
2.1.4. 概率密度函数 (3)
2.1.5 相关函数 (3)
2.2随机信号的特征 (4)
2.3平稳随机信号 (4)
2.3.1 平稳随机信号的定义 (4)
2.3.2 平稳随机信号的自相关函数 (5)
2.3.3 平稳随机信号的功率谱 (6)
2.4估计质量的评价标准 (6)
2.5离散时间信号—序列 (7)
3随机信号的功率谱估计 (8)
3.1周期图法 (8)
3.1.1 周期图法的简介 (8)
3.1.2 周期图法的定义 (9)
3.1.3 周期图的性能 (9)
3.1.4 周期图法改进措施 (10)
3.2相关法谱估计（BT） (10)
3.3直接法和间接法的关系 (12)
3.4DFT应用中的问题与参数选择 (12)
3.4.1 混叠现象 (12)
3.4.2 栅栏效应 (12)
3.4.3 截断效应 (13)
4 MATLAB软件介绍 (13)
4.1MATLAB产生的历史背景 (13)
4.2基本功能 (14)
4.3应用 (14)
4.4特点优势 (15)
5估计方法编程与分析 (15)
5.1直接法 (15)
5.2确定信号的谱分析 (19)
5.3方法的总结 (20)
6 结论 (21)
致谢 (21)
参考文献 (21)
引言
对信号和系统进行分析研究、处理有两种方法:一类是在时域进行，一类是在频域进行。这两类方法都是信号处理的重要方法。功率谱是无限多个自相关函数的函数，单观测数据只有有限个，只能得到有限个自相关函数。要无限个观测数根据有限个样本数据，分析计算随机序列的真正功率谱，是求功率谱的中心问题，毫无疑问，这是一个功率谱的估计问题。
功率谱估计是随机信号处理的重要内容，功率谱的估计方法有很多，主要有经典估计，也称线性估计（有偏估计）。经典谱估计有可以分成两种，一种是BT法，也叫间接法或自相关法；另一种是周期图法。周期图法简单，不用估计自相关函数，且可以用FFT进行计算。因此，周期图法得到了更广泛的应用。经典谱估计致命的缺点是频率分辨率低，其原因是傅里叶变换域是无限大，而观察数据只能得到有限个，观察不到的数据被认为是0。如果只有N个观测数据，而对于N以外的数据，信号仍有较强的相关性，这样估计出功率谱会出现很大的偏差。对于有限个观测数据，相当于信号在时域中乘以矩形窗函数，在频域则使真正的功率谱卷积一个sinc函数，而sinc函数不同于 函数，它有主瓣和旁瓣，这样使卷积后的功率谱不同于真正的功率谱。sinc函数的主瓣引起功率谱向附近频域展宽，通常称这种展宽为泄露。显然泄露造成谱的模糊，降低谱的分辨率；而旁瓣则会引起谱间干扰，信号强的功率谱旁瓣影响信号弱的功率谱检测，严重时检测不出弱信号，或者把旁瓣误认为是信号，造成假信号。为了对经典谱估计进行改进，可以采用各种不同的窗函数，但都是以增加主瓣宽度来换取旁瓣的压低，因此谱分辨率低是经典谱估计的致命缺点。
1绪论
1.1 功率谱研究的发展过程
功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计的发展过程做简要回顾：英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅里叶提出了著名的傅里叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。
傅里叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅里叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”（periodogram）。这是经典谱估计的最早提法，这种提法至今仍然被沿用，只不过现在是用快速傅里叶变换（FFT）来计算离散傅立叶变换（DFT），用DFT的幅度平方作为信号中功率的度量。
周期图较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列。Y ule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。
Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule 和Walker都是开拓自回归模型的先锋。
1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即“功率谱密度是随机过程二阶统计量自相关函数的傅里叶变换”，这就是Wiener—Khintchine定理。该定理把功率谱密度定义为频率的连续函数，而不再像以前定义为离散的谐波频率的函数。
1949年，Tukey根据Wiener—Khintchine定理提出了对有限长数据进行谱估计的自相关法，即利用有限长数据估计自相关函数，再对该自相关函数求傅里叶变换，从而得到谱的估计。1958年，Blackman和Tukey在出版的有关经典谱估计的专著中讨论了自相关谱估计法，所以自相关法又叫BT法。周期图法和自相关法都可用快速傅里叶变换算法来实现，且物理概念明确，因而仍是目前较常用的谱估计方法。
1948年，Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下了良好的理论基础。1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。
1.2 功率谱估计方法提出
在通信系统中，往往需要研究具有某统种计特性的随机信号。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数却是确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。
1.3 功率谱估计应用及用途
功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)。例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。根据信号与噪声的功率谱(或()xx R m )才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器。
常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性)(ωH 时，可用一白噪声()n ω通过该系统。再从该系统的输出样本()y n 估计功率谱密度)(ωyy P 。
由于白噪声的功率谱密度(用)(ωωωP 表示)为一常数即2
)(ωωωσω=P ，于是有： ()()22
YY P w H w σ
= (1.1) 故通过估计输出信号的功率谱密度，可以估计出系统的频率特性)(ωH (模特性)。从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向。
功率谱估计就是通过信号的相关性估计出接收到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等（火控系统中的飞机航迹预判）。 1.4 本课题研究任务
本课题在于掌握周期图法的原理，在此基础之上编写MA TLAB 程序，产生一随机信号，用周期图法对随机信号进行谱估计，完成论文。
2 随机信号简介
2.1 随机变量
随机变量（random variable ）表示随机现象（在一定条件下，并不总是出现相同结果的现象称为随机现象）各种结果的变量（一切可能的样本点）。例如某一时间内公共汽车站等车乘客人数，电话交换台在一定时间内收到的呼叫次数等等，都是随机变量的实例。
随机变量在不同的条件下由于偶然因素影响，其可能取各种不同的值，具有不确定性和随机性，但这些取值落在某个范围的概率是一定的，此种变量称为随机变量。随机变量可以是离散型的，也可以是连续型的。如分析测试中的测定值就是一个以概率取值的随机变量，被测定量的取值可能在某一范围内随机变化，具体取什么值在测定之前是无法确定的，但测定的结果是确定的，多次重复测定所得到的测定值具有统计规律性。随机变量与模糊变量的不确定性的本质差别在于，后者的测定结果仍具有不确定性，即模糊性。
按照随机变量可能取得的值，可以把它们分为两种基本类型：①离散型随机变量，即在一定区间内变量取值为有限个，或数值可以一一列举出来。例如某地区某年人口的出生数、死亡数，某药治疗某病病人的有效数、无效数等。②连续型随机变量，即在一定区间内变量取值有无限个,或数值无法一一列举出来。例如某地区男性健康成人的身长值、体重值，一批传染性肝炎患者的血清转氨酶测定值等。
2.1.1 随机变量的分布函数
设X 是随机变量，对任意实数x ，事件{}X x <的概率{}P X x <称为随机变量X 的分布函数。记为()F x ，即(){}F x P X x =<，易知，对任意实数(),,a b a b <
{}{}{}()()P a X b P X b P X a F b F a <<=<-<=- （2.1）
分布函数的性质
1.单调不减性：若12X X <， 则()()12F X F X <；
2.归一性：对任意实数x ，()01F x <<，且
()()()()0,1lim lim x x F F x F F x →-∞
→+∞
-∞==+∞== （2.2）
3.左连续性：对任意实数x ，
()()0F x F x -= (2.3）
2.1.2 数学期望、方差、标准差
定义：{}()x u E x xp x dx +∞
-∞
==
?
，为X 的数学期望值，或简称为均值。
{}()2
2
2x
D E x
x p x dx +∞
-∞
==?
（2.4） {}()2
2
2x
x
x E x u x u p x dx σ+∞
-∞
=-=-?
（2.5）
以上分别称为X 的标准差和方差。若X 为离散型随机变量，则上述的求均值运算将有积分改为求和。例如{}x k k k
u E x x p ==∑
，式中的k P 是X 取值为k x 时的概率。
2.1.3 随机向量
在某些实际问题中，往往需要同时用两个或两个以上的随机变量来描述试验的结果。 设E 是一个随机试验, 样本空间是{}e Ω=，设()X X e =和()Y Y e =是定义在Ω上的随机变量，由它们构成的一个向量(),X Y 叫做二维随机向量或二维随机变量。
（注: 二维随机向量 (),X Y 的性质不仅与X 和Y 有关,而且还依赖于这两个随机变量的相互关系。）
2.1.4. 概率密度函数
概率密度函数是为了表示瞬时数据落在指定幅值范围的概率。
其定义为：()()0
01lim
lim
lim x x x T prob x x t x x T p x x
x T ?→?→→∞<≤+????
???
=
=??????
(2.6) 瞬时值()x t 小于或等于某值x 的概率定义为概率分布函数或累计概率分布函数
()()()x
p x prob x t x p d ξξ-∞
=≤=
?????
(2.7)
2.1.5 相关函数
表征了一个随机过程自身在不同时刻的状态间，或者两个随机过程在某个时刻状态间线性依从关系的数字特征。
相关函数是两随机变量之积的数学期望，称为相关性。统计学中用相关系数xy 来描述变量x ，y 之间的相关性，函数的相关系数，简称相关函数：
()()()(){
}
1
2
2
2x y xy
xy x y
x y E x u y u E x u E y u σρσσ??--??
==
????--??????
（2.8）
2.2 随机信号的特征
随机信号
[]
1----具有不重复性、不确定性，通常用概率与统计方法研究其中是否存在某些重复、
确定的成分。
随机过程在某一时刻1t 的均值（一阶矩）可将总体中各样本函数在1t 的瞬时值相加，然后除以
样本函数的个数而得到 ()()111
1lim N
x k N k u t x t N →∞==∑ (2.9)
自相关函数即为随机过程两不同时刻之值的相关性，又称二阶矩。用1t 和1t +两时刻瞬时值乘积的总
体平均值得到()()()11111
1,lim N
xx k k N k R t t x t x t N ττ→∞=+=+∑ （2.10）
自相关函数的性质：
(1) 自相关函数是τ的偶函数
()()x x R R ττ=-; （2.11）
(2) 当0τ=时，自相关函数具有最大值
()2
2
0x x x R u σ=+； （2.12）
(3) 周期信号的自相关函数仍然是同频率的周期信号，但不保留原信号的相位信息。
(4) 当随机信号中含有周期信号时，()x R τ中也必定有周期性分量，且周期相同。 (5) 对变化迅速的信号（宽带随机过程），相关的程度在τ很小时就完全丧失 。 2.3 平稳随机信号
2.3.1 平稳随机信号的定义
平稳信号分严平稳和宽平稳，严平稳的条件在信号处理中太严格，不实用，一般所说的平稳是指宽平稳，满足三个条件：
(1) 均值为与时间无关的常数； (2) 均方有界；
(3) 自相关函数与信号时间的起始点无关，只和时间差有关（宽平稳信号的方差和均方也是与时间无关的）。
1. 平稳随机过程[]
2的定义
如果对于任意n 和12,,...,n t t t 以及τ有
()()12121212,,...,,,...,,,...,,,...,n n n n n n f x x x t t t f x x x t t t τττ=+++ （2.13）
则称()X t 为严平稳随机过程，或称狭义平稳随机过程。
2. 平稳随机过程的数字特征
(1)()X
E X t m =????，平稳随机过程的数学期望与时间无关； (2)()2
x D X t σ=????，平稳随机过程的方差与时间无关；
(3)()()()121221212,,,X X R t t x x f x x dx dx R ττ+∞+∞
-∞
-∞
=
=??
其中：12t t τ-=； (2.14)
(4) ()()()2
12,X X X
X C t t R m C ττ=-=。 (2.15) 平稳随机过程的数学期望及方差与无关，它的自相关函数和协方差函数只与时间间隔
有关；随机过程的这种“平稳”数字特征，有时就直接用来判断随机过程是否平稳。即
若一个随机过程的数学期望及方差与时间无关，而其相关函数仅与有关，即我们就称这个随机过程是广义平稳的。
3. 宽平稳随机过程（广义平稳） 若()X t 的数学期望()X E X t a =????为常数，且自相关函数()()12,X X R t t R τ=只与12t t τ-=有关，则称()X t 为宽平稳随机过程，或称广义平稳随机过程。不难看出，严平稳过程一定是宽平稳过程，反之不一定。但对于正态随机过程两者是等价的。本论文若不加特别说明，平稳过程均指宽平稳过程。
2.3.2 平稳随机信号的自相关函数
实随机信号()X t 的自相关函数定义
()()()()1212122121212,,;,xx R t t E X t X t x x f x x t t dx dx +∞
+∞
-∞
-∞
==?????
?
(2.16)
由于平稳随机信号的统计特性与时间的起点无关，设21t t t =+，有
()()21212212,;,,;f x x t t f x x t = (2.17)
所以，平稳随机信号的自相关函数是时间间隔t 的函数，记为()xx R t 。
平稳随机信号自相关函数的性质：设()X t 为平稳随机过程，其自相关函数为()xx R t ，自协方
差函数()XX C t ，则它们有如下性质[]
2：
（1）0τ=时的自相关函数等于均方差，自协方差函数等于方差， 即
()()()2
2
00xx xx x
R E X t C σ
??=??
= (2.18)
（2）当平稳随机信号是实函数时，其相关函数是偶函数，即：
()()
()()
x x x x
x x x x R R C C
ττττ=-=- (2.19) （3）0τ=时的自相关函数、自协方差函数取最大值，即
()()()()
00xx xx xx xx R R C C ττ≥= (2.20)
（4）若()()X t X t T =+，则其自相关函数也是周期为T 的周期函数，即
()()()()
xx xx xx xx R R T C C T ττττ≥+=+ (2.21)
若均值0mx =，当τ→∞时，()X t 与()X t τ+相互独立，有
()0lim xx
R τ
τ→∞
=，即对于零均值的平稳随机信号，当时间间隔τ很大时，()X t 与()X t τ+相互独立，互不相关。
2.3.3 平稳随机信号的功率谱 1. 功率谱密度
定义： 设(){}
,X t t -∞<<+∞是均方连续的随机过程，称
()21
2lim T
T T P E X t dt T -→∞??=?
???
?为()X t 的平均功率。称()()21,2lim X X T S w E F w T T →∞??=??为()X t 的功率谱密度，简称谱密度。
2. 功率谱密度的性质
（1）若
()X R d ττ∞
-∞
<∞?
，则()X S w 是()X R τ的傅里叶变换；
()()()()12iw X
X iw X X S w R e d R S w e dw ττττ
τπ∞--∞∞
-∞?=???=?
?? (2.22) (2) ()X S w 是w 的非负实函数；
(3) 实平稳过程的谱密度是偶函数； 当()X S w 是w 的有理函数时，其形式必为
()2222220
222220
......n n n n X m
m m a w a w a S w w b w b ----+++=+++ (2.23) 其中2n a ，2m b 为常数，且20,n a m n >>，分母无实根。 3.随机序列的功率谱
随机序列X(n)，它的相关函数()x R m 满足其功率谱密度()x S m 具有如下式子：
()()()()12iwm
X X m iwm X X S w R m e R m S w e dw π+∞
-=-∞
∞
-∞?=???
?=??
∑? (2-24) 2.4 估计质量的评价标准
1. 无偏性
对于待估参数，不同的样本值就会得到不同的估计值，这样，要确定一个估计量的好坏，就不能仅仅依据某次抽样的结果来衡量，而必须由大量抽样的结果来衡量。对此，一个自然而基本的衡量标准是要求估计量无系统偏差。也就是说，尽管在一次抽样中得到的估计值不一定恰好等于待估参数的真值，但在大量重复抽样时，所得到的估计值平均起来应与待估参数的真值相同。换句话说，我们希望估计量的均值（数学期望）应等于未知参数的真值，这就是所谓无偏性(Unbiasedness)的要求。
定义：设来自总体X 的一个样本，是总体参数θ的一个估计量，若E θθ∧
=，则称θ∧
是θ的无偏估计量(Unbiased Estimator)。一个估计量如果不是无偏的就称它是有偏估计量。E θθ∧
-称为估计量θ∧
的偏差。无偏估计的实际意义就是无系统偏差，估计量是否无偏是评价估计量好坏的一个重要标准，若E θθ∧
≠，但有
lim n E θθ∧→∞
=，则称θ∧
是θ的渐近无偏估计。
2. 有效性
比较两个无偏估计量优劣的一个重要标准就是观察它们哪一个取值更集中于待估参数的真值附近，即哪一个估计量的方差更小，这就是下面给出的有效性(Effectiveness)概念。
定义：设与()2212,,...n X X X θθ∧∧
=都是总体参数θ的无偏估计，若12D D θθ∧∧
????≤ ? ?????，则称1θ∧比2θ∧
更有效。在θ的所有无偏估计量中，如果存在一个估计量0θ∧
，它的方差最小，则此估计量应当最好，并称此估计量0θ∧
为θ的最小方差无偏估计，也称其为最有效的。 3. 相合性
估计量θ∧
的无偏性和有效性都是在样本容量n 固定的情况下讨论的。由于估计量θ∧
和样本容量n 有关，我们自然希望当n 很大时，一次抽样得出的的θ值能以很大的概率充分接近被估参数θ，这
就提出了相合性(Consistency)（一致性）的要求。
定义：
设()12,,...n X X X θθ∧
∧
=是总体参数θ的估计量，如果对任意0ε>都有1lim n P θθ∧→∞??
-=????
， 则
称θ∧是θ的相合估计量（或一致估计量）。θ是θ∧的相合估计就意味着θ∧
依概率收敛于θ。根据大数
定律，无论总体X 服从什么分布，只要其k 阶原点矩()
k
k u E X =存在，则对任意0ε>都有
111lim n i n i P X EX n ε→∞=??
-<=????
∑ (2.25)
所以样本的k 阶原点矩1
1n k
k i i A X n ==∑始终是总体k 阶原点矩k u 的相合估计。 进一步地， 可以证明：
只要相应的总体矩存在，矩估计必定是相合估计。特别地，u X ∧
-
= 总是u EX =的相合估计， 样
本方差2S 和样本的二阶中心矩2
n S 都是总体方差2
σ的相合估计S 和s 又都是σ的相合估计。由相合
性定义可以看出，若θ∧
是θ的相合估计，当样本容量很大时，一次抽样得到的θ∧
值便可作为θ的较好近似值。
2.5 离散时间信号—序列
离散时间序列可由连续时间信号的采样得到[]
3。
采样过程是通过采样脉冲序列()p t 与连续时间信号()x t 相乘来完成的，采样过程如下图所示。
其采样脉冲序列
()()s
n p t t nT δ∞
=-∞
=
-∑ (2.26)
采样信号
()()()s x t x t p t = (2.27)
如果()()FT x t X w =????，()()FT p t P w =????那么，根据频域卷积定理，有
()()()
2s X w P w X w π
*=
采样脉冲序列()p t 的频谱是间隔为s w 的周期延拓，所以
()()s s
n X w X w nw ∞
=-∞
=
-∑此式表明，对连续信号进行等间隔采样形成采样信号，采样信号的频谱
是原连续信号的频谱以采样频率为周期进行周期性的延拓形成的。
3 随机信号的功率谱估计
3.1 周期图法
3.1.1 周期图法的简介
周期图法又称直接法。周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换
[]
3，
然后取其幅频特性的平方并除以序列长度N 。由于序列()x n 的离散傅里叶变换()X K 具有周期性，因而这种功率谱也具有周期性，常称为周期图。它是从随机信号()x n 中截取N 长的一段，把它视
为能量有限()x n 真实功率谱()
jw x S e 的估计()
jw
x S e 的抽样。
周期图这一概念早在1899年就提出了，但由于点数N 一般比较大，该方法的计算量过大而在
当时无法使用。只是1965年FFT []
3出现后，此法才变成谱估计的一个常用方法。周期图法包含了下列二条假设：
(a )
(b )
x a )(c )
(e )(d )
( f )
1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本()x n 中的一段()N x n 来估计该随机序列的功率谱。这当然必然带来误差。
2．由于对()N x n 采用DFT ，就默认()N x n 在时域是周期的，以及()N x k 在频域是周期的。这种方法把随机序列样本()x n 看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数()x R m 时将()N x n 以外是数据全都看成零，因此相关法认为除()N x n 外()x n 是全零序列，这种处理方法显然与周期图法不一样。
但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。我们发现：如果相关法中M N =，不加延迟窗，那么就和补充()1N -个零的周期图法一样了。简单地可以这样说：周期图法是M N =时相关法的特例。因此相关法和周期图法可结合使用。 3.1.2 周期图法的定义
周期图法，它是把随机序列()x n 的N 个观测数据视为一能量有限的序列，直接计算()x n 的离散傅里叶变换，得()X K ，然后再取其幅值的平方，并除以N ，作为序列()x n 真实功率谱的估计。
周期图谱估计定义为：
()()()
2
1
1
exp 2N PER n P f x n j fn N
π-∧
==
-∑ (3.1)
可以证明，周期图等于估计出的自相关序列的傅里叶变换，或
()()
()()1
1exp 2N PER xx k N P f r k j fk π-∧
∧
=--=
-∑
(3.2)
其中()xx r k ∧
是有偏自相关函数的估计值，定义为： ()()()10
1
N k
xx n r k x n x n k N
--∧
==
+∑ (3.3)
3.1.3 周期图的性能
周期图的期望值是：
()()(){}()()12xx -12
=P d PER B xx B E P f F w k r k W f ξξξ∧??
=-????
? (3.4)
式中()B W f 是Bartlett 窗（三角形窗）的傅里叶变换。由式（3-4）可知周期图的均值是真实PSD 和
Bartlett 窗傅里叶变换的卷积，在平均意义上得到真实功率谱密度（PSD ）的平滑形式。因此对有限记录数据，周期图一般有偏的，但是当N →∞时，它是无偏的。
()()xx N P lim PER E P f f ∧→∞
??=????这是由于()B W f 收敛到狄拉克δ函数。 对于高斯白噪声的特殊情况，()()()2
2
xx ,P xx x x r k k f σδσ==，结果为：
()()2
xx P PER x E P f f σ∧??==????
(3.5) 对于白噪声情况，即使有限记录数据，周期图也是无偏的。对于白高斯过程，
()()2
xx P PER x E P f f σ∧
??==????，方差()()2
2xx sin 2var P 1sin 2PER Nf P f f N f ππ∧
??????=+?? ??????????? (3.6) 对任何非白过程，只要记录数据足够长，
()()2
2
xx sin 2var P 1sin 2PER Nf P f f N f ππ∧
??????≈+?? ???????????
(3.7)
对于不靠近0或1
2
±
的频率，且N →∞时，上式近似地退化为： ()()2xx var P PER P f f ∧??≈????
(3.8)
可以看出，方差与数据长度N 无关，即方差不随N 的增大而减小至零。由此可得出一个重要的
看法：周期图估计器是不可靠的，因为标准差和均值一样大，因而周期图不是一致估计而其均值近似地等于要估计的量值。
上述论证表明，我们不能寄希望于直接用周期图方法获得良好的谱估计，必须采用适当的修正措施减小估计方差，才能使之成为一种实用的方法。 3.1.4 周期图法改进措施 1. 加窗周期图
周期图法只用了N 个样本，这可以看作是用一长度为N 的矩形窗函数与原来无限长的序列相乘的结果，我们知道，时域中两函数相乘对应于频域中它们的傅立叶变换的卷积。由此可以想到，当用周期图方法作谱估计时，它的谱分辨率约与N 成反比，且和信号本身的特征例如信噪比等无关。此外，如果序列是由多个正弦波信号组成的，而各分量强度不等，则弱信号分量可能淹没在强信号谱的旁瓣中而无法发现。这种所谓信号能量（向旁瓣）泄漏现象如果不设法消除，也将妨碍周期图谱估计法的应用。因此提出了周期图的改进方法：
周期图改进的方法[]
1之一是将长度为N 的序列()x n 乘以同一长度的数据窗()w n 。
数据加窗后的周期图谱估计值的数学期望值等于谱的真实值与窗谱函数的平方的卷积。显然它不等于功率谱的真实值，因而是有偏估计。
若序列为单频信号，则()xx P f 为δ函数，这样，数据加窗后的谱估计值的均值与窗谱函数的平方形状相同，因此选用低旁瓣的数据窗可使得杂散响应减少。但旁瓣的降低必然使主瓣加宽，而且降低了分辨率。
数据加窗后，周期图谱估计值的方差大于或近似等于谱估计值的平方，且不随数据长度的增大而减小到0。
从以上的分析可知，数据加窗用于周期图谱估计可以降低谱估计值的旁瓣，但要降低谱估计的分辩率，而用数据加窗的办法不能减小估计方差，因而无法降低分辨率。 3.2 相关法谱估计 （BT ）
这种方法以相关函数为媒介来计算功率谱，所以又叫间接法[]
1。自相关函数是描述随机信号()x t 在任意两个不同时刻12,t t 的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个判断指标，把两测点之间信号的互谱与各自的自谱联系了起来。它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生的误差非常有效。维纳-辛钦定理：随机信号的相关函数与其功率谱是一组傅里叶变换对,即相关函数的傅里叶变换为功率谱,而功率谱的逆傅里叶变换为相关函数。该定理揭示了从时间角度描述随机信号统计规律和从频率角度描述随机信号统计规律之间的联系。它是1958年由Blackman 和Tukey 提出。这种方法的具体步骤是：
第一步：从无限长随机序列()x n 中截取长度N 的有限长序列列()N x n 第二步：由N 长序列()N x n 求()21M -点的自相关函数()x R m ∧
序列。即
()()()1
1N x N N n R m x n x n m N -∧
==+∑ (3.9)
这里，()1,...1,0,1...1,m M M M N =----≤，()x R m 是双边序列，但是由自相关函数的偶
对称性式，只要求出0,...,1m M =-的傅里叶变换，另一半也就知道了。
第三步：由相关函数的傅式变换求功率谱。即
()()(
)
()
11M jw
jwm x x m M S e
R m e -∧
∧
-=--=∑
(3.10)
以上过程中经历了两次截断，一次是将()x n 截成N 长，称为加数据窗，一次是将()x n 截成()21M -长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M N <<，因为只有当M 较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT 问世之前，相关法是最常用的谱估计方法。
当FFT 问世后，情况有所变化。因为截断后的()N x n 可视作能量信号，由相关卷积定理可得
()()()1
x N N R m x m x m N
∧
=
*-???? (3.11) 这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上式两边取（2N-1）
点DFT ，则有
()()()()212212111
N x N N R m x K x K X K N N
-∧
--??=
*=?? (3.12) 于是将时域卷积变为频域乘积，用快速相关求)(m R x ∧
的完整方案如下： (1) 对N 长()N x n 的充()1N -个零，成为()21N -长的。 (2) 求()21N -点的FFT ，得()()21
2221210
N N mk
N N N X K x n W
-----==∑。
(3.13)
(3) 求
()2211
N X K N
-。由DFT 性质，()21N x n -是纯实的，()21N x k -满足共轭偶对称，而()2211
N X K N
-一定是实偶的，且以()21N -为周期。 (4) 求()21N -点的IFFT ：
()()()1
22121111
21N mk x N N K N R m X K W N N -∧
---=--=-∑ (3.14)
这里
()2211
N X K N
-是实偶的，()1...0...1m N N =---。本来IFFT 求和范围是0至22N -，由于()2211
N X K N
-的实偶性与周期性，求和范围改为()1N --至()1N -不影响计算结果。同理可将m 的范围改为()1N --）至()1N -。
上述的快速相关中，充零的目的是为了能用圆周卷积代替线性卷积，以便进一步采用快速卷积算法。快速相关输出是()1N --至()1N -的2N-1点，加)(m W M 窗后截取的是()1M --至()1M -的频段，最后作()21M -点FFT ，得()S k ∧
。我们注意到：如果数据点数与自相关序列点数相同即
M N =，则()21N -点的IFFT 后紧跟一个()21N -点的FFT ，利用()x R m ∧
的对称性，FT 运算
框的计算式变为
()()
()1
21211N mk X X N N m N S K R m W X -∧
---=--=
∑
(3.15)
由于M N =并假设窗形状是矩形的，第二次()m W M 的截断就不需要了。比较式（3.5）和式（3.6），
()()()()2211,x x X N S K FFT R m R m IFFT X K N ∧∧
-????==????????
,正反傅氏变换可以抵消，直接得
()()2211
X N S K X K N
-=
(3.16) 为了实行基2FFT ，也可将()21N -点换成2N 点，这样做不影响结果的正确性。 3.3直接法和间接法的关系
由以下式子可以得出： ()()()()()()()()()
2110
11
0011
1
1jw jw jw N N N jwk
jwn
k n N N jw k n k n I w X e X e X e N N
x k e x n e N
x k x n e
N
*----==----===
==
=∑∑∑∑ (3.17)
m k n =-，则k n m =+，得
()()()()()()1
11
10
11
N N m
N jwm
jwm N N m N n m N I w x n x n m e R m e N
------=--==--??=+=?
???
∑∑
∑ (3.18)
可见周期图法就是BT 法中1M N =-时的功率谱估计
3.4 DFT 应用中的问题与参数选择
在应用DFT 解决实际问题时常遇到下列几个问题[]
4：（1）混叠现象（2）栅栏效应（3）截断效应
这些问题与应用中的信号和DFT 的参数选择有关。 3.4.1 混叠现象 1. 混叠产生的机理
混叠主要由采样引起，可以指时域混叠，也可以指频域混叠。因为对()x t 带限信号，最高频率为c f ,以采样频率s f 进行采样，所得采样序列的频谱是以s f 为周期，进行周期延拓而成的。所以，处理连续时间信号前，常令其通过截止频率为h f 的前置低通滤波器，保证滤波器输出信号的双边谱带宽为2h f ，即得到最高频率h f ≤的模拟信号()a x t 。如果满足2s c f f ≥，基带谱与其他周期延拓形成的谱不重叠，但如果选择的采样频率太低，或者是信号最高截止频率过高，使2s c f f ≤，则以采样频率s f 周期延拓时，形成频谱混叠现象。 2. 避免混叠现象、保证处理要求的主要措施
理论上必须满足2s c f f ≥（c f 为连续信号的最高频率）。对s f 确定的情况下，一般在采样前进行预滤波，滤除高于折叠频率2s f 的频率成分，以免发生频谱混叠现象。 3.4.2 栅栏效应
1. 栅栏效应产生的机理
非周期信号的频谱是连续的。对非周期信号N 点采样进行DFT 得到的频谱只能是连续频谱上的有限离散频点采样。 2. 栅栏效应
N 点DFT 是在频率区间[]0,2π上对时域离散信号的频谱进行N 点等间隔采样，而采样点之间的频谱是看不到的。这就好像从N 个栅栏缝隙中观看信号的频谱情况，仅得到N 个缝隙中看到的频谱函数值。这种现象就是栅栏效应。不管是时域采样还是频域采样，都有相应的栅栏效应。只是当时域采样满足采样定理时，栅栏效应不会有什么影响。而频域采样的栅栏效应则影响很大，“挡住”或丢失的频率成分有可能是重要的或具有特征的成分，使信号处理失去意义。 3. 克服的措施
为了把被“栅栏”挡住的频率分量检测出来，对有限长序列，可在原序列尾部补零；对无限长序列，可以增大截取长度及DFT 变换区间长度，从而使频域采样间隔变小，增加频域采样点数和采样点位置，使原来漏掉得的某些频率分量被检测出来。对连续信号的谱分析，只要采样速率足够高，且采样点数满足频率分辨率要求，就可认为DFT 后所得离散谱的包络近似代表原信号的频谱。间隔小，频率分辨力高，被“挡住”或丢失的频率成分就会越少。但会增加采样点数，使计算工作量增加。解决此项矛盾可以采用如下方法：在满足采样定理的前提下，采用频率细化技术(ZOOM ），亦可用把时域序列变换成频谱序列的方法。 3.4.3 截断效应
实际中遇到的序列可能是无限长的，用DFT 对其进行谱分析时，必须将其截短，形成有限长序列。截断后序列的频谱与原序列的频谱必然有差别。这种差别对谱的影响主要表现为泄露和谱间干扰两个方面。
（1）泄露。信号经系统处理后，所得的频谱在原来没有频谱的区间出现了频谱。显然，泄露使频谱变模糊，使谱分辨率降低。具体的泄漏分布取决于所采用的窗的连续域傅里叶变换，对于没有使用窗的，相当于使用了矩形窗，而对于非矩形窗，窗本身就会产生一定的泄漏，频谱泄露程度与窗函数幅度谱的主瓣宽度直接相关。在所有窗函数中，矩形窗的主瓣是最窄的，但其旁瓣的幅度也最大。
（2）谱间干扰。在主谱线两边形成很多旁瓣，引起不同频率分量间干扰，强信号谱的旁瓣肯能湮没弱信号的主谱线，或是把强信号谱的旁瓣误认为是另一频率的信号的谱线，从而形成假信号。这样就会使谱分析产生较大偏差。
克服的措施：增加N 使主瓣变窄，减小泄露，提高频率分辨率，但旁瓣的相对幅度并不减小。为了减小谱间干扰，应使用形状合适的窗函数。总之，减小截断效应的最好方法是近代谱估计的方法。
4 MATLAB 软件介绍
4.1 MATLAB 产生的历史背景
在70年代中期，CleveMoler 博士和其同事在美国国家科学基金的资助下开发了调用EISPACK 和LINPACK 的FORTRAN 子程序库。EISPACK 是特征值求解的FOETRAN 程序库，LINPACK 是解线性方程的程序库。在当时这两个程序库代表矩阵运算的最高水平。
到70年代后期,身为美国NewMexico 大学计算机系系主任的CleveMoler,在给学生讲授线性代数课程时，想教学生使用EISPACK 和LINPACK 程序库，但他发现学生用FORTRAN 编写接口程序很费时间，于是他开始自己动手，利用业余时间为学生编写EISPACK 和LINPACK 的接口程序。CleveMoler 给这个接口程序取名为MATLAB ，该名为矩阵(matrix)和实验室(labotatory)两个英文单词的前三个字母的组合。在以后的数年里。MATLAB 在多所大学里作为教学辅助软件使用，并作为面向大众的免费软件广为流传。 1983年春天，Cleve Moler 到Standford 大学讲学，MATLAB 深深地吸引了工程师John Little.John Little 敏锐地觉察到MA TLAB 在工程领域的广阔前景。同年他和Cleve Moler SteveBangert 一起，用C 语言开发了第二代专业版。这一代的MA TLAB 语言同时具备了数值计算和数据图示化的功能。1984
年CleveMoler和JohnLittle成立了MathWorks公司，正式把MATLAB推向市场，并继续进行MA TLAB 的研究和开发。
在当今30多个数学类科技应用软件中，就软件数学处理的原始内核而言，可分为两大类：一类是数值计算型软件，如MATLAB、Xmath、Gauss等，这类软件长于数值计算，对处理大批数据效率高。另一类是数学分析型软件，Mathematica、Mapl等，这类软件以符号计算见长，能给出解析解和任意精确解，其缺点是处理大量数据时效率较低。Math Works公司顺应多功能需求之潮流，在其卓越数值计算和图示能力的基础上，又率先在专业水平上开拓了其符号计算。文字处理，可视化建模和实时控制能力，开发了适合多学科多部门要求的新一代科技应用软件MATLAB。经过多年的国际竞争，MATLAB以经占据了数值软件市场的主导地位。
在MA TLAB进入市场前，国际上的许多软件包都是直接以FORTRANC语言等编程语言开发的。这种软件的缺点是使用面窄，接口简陋，程序结构不开放以及没有标准的基库，很难适应各学科的最新发展，因而很难推广。MA TLAB的出现，为各国科学家开发学科软件提供了新的基础。在MATLAB问世不久的80年代中期，原先控制领域里的一些软件包纷纷被淘汰或在MA TLAB上重建。
MathWorks公司1993年推出了MATLAB4.0版，1995年推出4.2C版1997年推出5.0版。1999年推出5.3版。MA TLAB5 X版较MATLAB4，X版无论是界面还是内容都有长足的进展，其帮助信息采用超文本格式和PDF格式，在Netscape3.0或IE4.0及以上版本，AcrobatReader中可以方便地浏览。
时至今日，经过MathWorks公司的不断完善，MATLAB已经发展成为适合多学科，多种工作平台的功能强大的大型软件。在国外，MA TLAB已经经受了多年考验。在欧美等高校，MATLAB已经成为线性代数，自动控制理论，数理统计，数字信号处理，时间序列分析，动态系统仿真等高级课程的基本教学工具；成为攻读学位的大学生，硕士生，博士生必须掌握的基本技能。在设计研究单位和工业部门，MATLAB被广泛用于科学研究和解决各种具体问题。在国内，特别是工程界，MATLAB一定会盛行起来。可以说，无论你从事工程方面的哪个学科，都能在MA TLAB里找到合适的功能。
4.2 基本功能
MATLAB和Mathematica、Maple并称为三大数学软件[]5。它在数学类科技应用软件中在数值计算方面首屈一指。MA TLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MA TLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++ ，JA V A的支持。可以直接调用，用户也可以将自己编写的实用程序导入到MA TLAB 函数库中方便自己以后调用，此外许多的MA TLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
4.3 应用
MATLAB 产品族可以用来进行以下各种工作[]6：
（1）数值分析。
（2）数值和符号计算。
（3）工程与科学绘图。
（4）控制系统的设计与仿真。
（5）数字图像处理技术。
（6）数字信号处理技术。
（7）通讯系统设计与仿真。
（8）财务与金融工程。
在欧美大学里，诸如应用代数、数理统计、自动控制、数字信号处理、模拟与数字通信、时间
序列分析、动态系统仿真等课程的教科书都把MA TLAB 作为内容。这几乎成了九十年代教科书与旧版书籍的区别性标志。
在国际学术界，MATLAB 已经被确认为准确、可靠的科学计算标准软件。在许多国际一流学术刊物上，（尤其是信息科学刊物），都可以看到MA TLAB 的应用。如今MATLAB 已应用于下列领域：
数值和符号计算、工程与科学绘图、数值分析、建模和仿真、控制系统的设计与仿真、数字图像处理、数字信号处理、通讯系统设计与仿真、财务与金融工程。具体表现为：自动控制、航天工程、汽车工业、生物医学工程、语音处理、图像处理、雷达工程、信号分析、计算机技术等各行各业中。
4.4 特点优势
特点：
（1）此高级语言可用于技术计算。
（2）此开发环境可对代码、文件和数据进行管理。
（3）交互式工具可以按迭代的方式探查，设计及求解问题。
（4）数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数值积分等。 （5）二维和三维图形函数可用于可视化数据。 （6）各种工具可用于构建自定义的图形用户界面。
（7）各种函数可将基于MATLAB 的算法与外部应用程序和语言（如 C 、C++、Fortran 、Java 、COM 以及 Microsoft Excel ）集成。
（8）不支持大写输入，内核仅仅支持小写。 优势：
（1）友好的工作平台和编程环境。 （2）简单易用的程序语言。
（3）强大的科学计算机数据处理能力。 （4）出色的图形处理功能。
（5）应用广泛的模块集合工具箱。 （6）实用的程序接口和发布平台。 （7）应用软件开发（包括用户界面）。
5 估计方法编程与分析
5.1 直接法
直接法又称周期图法，它是把随机序列()x n 的N 个观测数据视为一能量有限的序列，直接计算()x n 的离散傅里叶变换，得()X k ，然后再取其幅值的平方，并除以N ，作为序列()x n 真实功率谱的估计。
(1) 周期图法
Matlab 代码示例： clear;
Fs=1000; %采样频率
n=0:1/Fs:1;%产生含有噪声的序列改变n 的取值范围观察图形的变换 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));
图5.1 周期图法仿真图
（2）加窗周期图
集平均周期图法与平滑周期图法的优点于一体，形成了平滑周期图平均法。方差大大减小，减少泄露，且功率谱是渐进无偏的。
以下是利用matlab对加窗周期图法[]8进行仿真，并得到仿真图。程序如下：
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(100);%矩形窗
window1=hamming(100);%汉明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
振动台在使用中经常运用的公式 1、 求推力（F ）的公式 F=（m 0+m 1+m 2+ ……）A …………………………公式（1） 式中：F —推力（激振力）（N ） m 0—振动台运动部分有效质量（kg ） m 1—辅助台面质量（kg ） m 2—试件（包括夹具、安装螺钉）质量（kg ） A — 试验加速度（m/s 2） 2、 加速度（A ）、速度（V ）、位移（D ）三个振动参数的互换运算公式 2.1 A=ωv ……………………………………………………公式（2） 式中：A —试验加速度（m/s 2） V —试验速度（m/s ） ω=2πf （角速度） 其中f 为试验频率（Hz ） 2.2 V=ωD ×10-3 ………………………………………………公式（3） 式中：V 和ω与“2.1”中同义 D —位移（mm 0-p ）单峰值 2.3 A=ω2 D ×10-3 ………………………………………………公式（4） 式中：A 、D 和ω与“2.1”，“2.2”中同义 公式（4）亦可简化为： A= D f ?250 2 式中：A 和D 与“2.3”中同义，但A 的单位为g 1g=9.8m/s 2 所以： A ≈D f ?25 2 ,这时A 的单位为m/s 2 定振级扫频试验平滑交越点频率的计算公式 3.1 加速度与速度平滑交越点频率的计算公式 f A-V = V A 28.6 ………………………………………公式（5） 式中：f A-V —加速度与速度平滑交越点频率（Hz ）（A 和V 与前面同义）。
3.2 速度与位移平滑交越点频率的计算公式 D V f D V 28.6103?=- …………………………………公式（6） 式中：D V f -—加速度与速度平滑交越点频率（Hz ）（V 和D 与前面同义）。 3.3 加速度与位移平滑交越点频率的计算公式 f A-D =D A ??2 3 )2(10π ……………………………………公式（7） 式中：f A-D — 加速度与位移平滑交越点频率（Hz ），（A 和D 与前面同义）。 根据“3.3”，公式（7）亦可简化为： f A-D ≈5× D A A 的单位是m/s 2 4、 扫描时间和扫描速率的计算公式 4.1 线性扫描比较简单： S 1= 1 1 V f f H - ……………………………………公式（8） 式中： S1—扫描时间（s 或min ） f H -f L —扫描宽带，其中f H 为上限频率，f L 为下限频率（Hz ） V 1—扫描速率（Hz/min 或Hz/s ） 4.2 对数扫频： 4.2.1 倍频程的计算公式 n=2Lg f f Lg L H ……………………………………公式（9） 式中：n —倍频程（oct ） f H —上限频率（Hz ） f L —下限频率（Hz ） 4.2.2 扫描速率计算公式 R= T Lg f f Lg L H 2/ ……………………………公式（10） 式中：R —扫描速率（oct/min 或）
功率谱定义 从确定性信号功率计算开始 ()()221 11lim lim 222T T T T T P x t dt X d T T ωωπ∞--∞→∞→∞==?? ()()21lim 2T T S X T ωω→∞= S(w)为功率谱密度，简称功率谱 则 ()12P S d ωωπ+∞-∞= ? 随机信号的功率谱密度 （1）样本功率谱与功率谱密度 ()()21,lim ,2X T T S X T ωξωξ→∞= 针对一个具体的样本而言，其是一个确定性的信号 (2) 随机信号的平均功率及平均功率谱密度 ()X X P E P ξ=???? 需要对具体的样本取概率均值才能计算出功率 ()()()21,lim ,2X X T T S E S E X T ωωξωξ→∞??==?????? 故功率谱密度是对所有概率取期望的反应。 （3）自相关函数与功率谱密度 ()()R S τω? （4）信号的自相关函数计算 分为确定信号和随机信号 确定信号 02002*0 1()lim ()()T T x T R x t x t dt T ττ-→∞=-? 周期信号 0202*0 1()()()T T x R x t x t dt T ττ-=-? 随机信号 *()[()()]x R E x t x t ττ=- 2 功率计算 （1）根据定义来计算
（2）周期信号如何计算 0cos()A t ω的计算 200()()1()[]2 A A s d T πσωωπσωωωω+∞-∞-++==?不好算因此放弃，但是应该可以类推得出结论 (3)自相关函数计算 0cos()A t ω的计算 /2 200/2 /222000/2201()cos()cos(())cos()cos(2)1[]2 cos()2 T T T T r A t t d T A A t d T A τωωτωωτωωτωωτ+-+-=-+-==?? 所以其功率谱为 200()2 A πσωωσωω（-）+（+） 0j t Ae ω的计算 0000/2()2/2 /22/2 21()1T j t j t T T j T j r A e e dt T A e dt T A e ωωτωτωτ τ+---+-===?? 总结：因此周期函数，首先转换成傅里叶级数，然后再通过自相关函数的定义计算自相关函数，得到其功率谱密度。
t=0:0.0001:0.1; %时间间隔为0.0001，说明采样频率为10000Hz x=square(2*pi*1000*t); %产生基频为1000Hz的方波信号 n=randn(size(t)); %白噪声 f=x+n; %在信号中加入白噪声 figure(1); subplot(2,1,1); plot(f); %画出原始信号的波形图 ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号'); y=fft(f,1000); %对原始信号进行离散傅里叶变换，参加DFT采样点的个数为1000 subplot(2,1,2); m=abs(y); f1=(0:length(y)/2-1)'*10000/length(y);%计算变换后不同点对应的幅值plot(f1,m(1:length(y)/2)); ylabel('幅值的模'); xlabel('时间(s)'); title('原始信号傅里叶变换'); %用周期图法估计功率谱密度 p=y.*conj(y)/1000; %计算功率谱密度 ff=10000*(0:499)/1000; %计算变换后不同点对应的频率值 figure(2); plot(ff,p(1:500)); ylabel('幅值'); xlabel('频率(Hz)'); title('功率谱密度（周期图法）'); 功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计
matlab实现功率谱密度分析psd及详细解说 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 直接法： 直接法又称周期图法，它是把随机序列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);
谱让人联想到的Fourier变换，是一个时间平均（time average）概念，对能量就是能量谱，对功率就是功率谱。 功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变化情况。保留了频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。 有两点需要注意： 1. 功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列） 2. 功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶矩是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。 频谱分析： 对动态信号在频率域内进行分析，分析的结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。 功率谱密度： 功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆(ohm)时的实际功率。
由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。维纳-辛钦定理（Wiener-Khinchin theorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。 信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。 随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。 功率谱与自相关函数是一个傅氏变换对。 功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。从名字分解来看就是说，观察对象是功率，观察域是谱域。 通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于一条直线。 一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。 1. 用相关函数的傅立叶变换来定义谱密度； 2. 用随机过程的有限时间傅立叶变换来定义谱密度； 3. 用平稳随机过程的谱分解来定义谱密度。 三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周
信号的功率谱分析 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 函数和周期信号的频谱——δ函数实行卷积，因此截断后的周期函数的频谱已不再是脉冲函数，原来为无限大的谱线高度变成有限长，谱
F F T在功率谱密度计算 中的应用 集团标准化工作小组 #Q8QGGQT-GX8G08Q8-GNQGJ8-MHHGN#
FFT在功率谱密度计算中的应用 一、FFT算法理论依据和编程思想 FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解） 先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组： 其中X 1（k）和X 2 （k）分别是x 1 （r）和x 2 （r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面 （1）式合成为一个N点DFT，注意到，X 1（k），X 2 （k）有N/2个点，即k=0，1，…， N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0， 1，…，N-1，要用X 1（k），X 2 （k）表示全部X（K）值，还必须应用系数w的周期性和 对称性。 (k)的(N/2)～N-1点表示： 由X(k)= X 1（k）+w k N X 2 （k）, k=0,1,2,…,N/2-1
随机信号的功率谱密度估计和相关函数
随机信号的功率谱密度估计和相关函数 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(…输入信号?) %绘输入信号图
功率谱密度的估计 原始波=余弦波+白噪声 这个实验采用了两个输入，一个是白噪声，一个是有用信号和噪声信号作为输入时，他们的功率谱密度的仿真图像，并将他们进行对比。 平稳随机信号的功率谱密度（PSD ）是相关序列的离散傅里叶变换： ()()jw m XX x P w r m e ∞ --∞=∑ 采用间接法计算噪声信号的功率谱。 间接法，又称自相关法或者BT 法，在1985年由布莱克曼与图基首先开拓。间接法的理论基础是维纳-辛钦定理。他是由N 个观察值x(0),x(1),……,x(N-1),估计出自相关函数R （m ）,然后再求R （m ）的傅里叶变换作为功率谱密度的估计。 ()(),||1M jw jw m N m M S e R m e M N -=-=<=-∑ clear all; randn('state',0) NFFT=1024; %采样点数 Fs=1000; %取样频率（单位为Hz ） t=0:1/Fs:.2;
y1=cos(t*20*pi); %余弦序列 figure(1) plot(t,y1); ylabel('余弦序列'); grid on; %余弦序列的图像： %白噪声 m=(0:NFFT-1)/Fs; y=0.1*randn(size(m)); %产生高斯白噪声。 figure(2); plot(m,y); title('白噪声波形'); grid on;
%白噪声的自相关函数 [cory,lags]=xcorr(y,200,'unbiased'); %计算白噪声的自相关函数 figure(3) plot(lags,cory); %自相关函数(无偏差的)，其中，cory为要求的自相关函数，lag为自相关函数的长度。 title('白噪声相关函数'); grid on;
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 ω为常数，Θ为在[]-,ππ上均匀分布的随机变量。
功率谱密度 不同形式的数字基带信号具有不同的频谱结构，分析数字基带信号的频谱特性，以便合理地设计数字基带信号，使得消息代码变换为适合于给定信道传输特性的结构，是数字基带传输必须考虑的问题。 在通信中，除特殊情况（如测试信号）外，数字基带信号通常都是随机脉冲序列。因为，如果在数字通信系统中所传输的数字序列是确知的，则消息就不携带任何信息，通信也就失去了意义。故我们面临的是一个随机序列的谱分析问题。 考察一个二进制随机脉冲序列。设脉冲、分别表示二进制码“0”和“1”, 为 码元的间隔，在任一码元时间内，和出现的概率分别为p和1-p。 则随机脉冲序列x(t)可表示成: 其中 研究由上面二式所确定的随机脉冲序列的功率谱密度，要用到概率论与随机过程的有关知识。可以证明，随机脉冲序列x(t)的双边功率谱公式（1）： 其中、分别为、的傅氏变换，。 可以得出如下结论： （1）随机脉冲序列功率谱包括两部分：连续谱（第一项）和离散谱（第二项）。对于连续谱而言，由于代表数字信息的及不能完全相同，故，因此，连 续谱总是存在；而对于离散谱而言，则在一些情况下不存在，如及是双极性的脉冲，且出现概率相同时。 （2）当、、p及给定后，随机脉冲序列功率谱就确定了。 上式的结果是非常有意义的，它一方面能使我们了解随机脉冲序列频谱的特点，以及如何去具体地计算它的功率谱密度；另一方面根据它的离散谱是否存在这一特点，将使我们明确能否从脉冲序列中直接提取离散分量，以及采取怎样的方法可以从基带脉冲序列中获得所需的离散分量。这一点，在研究位同步、载波同步等问题时，将是十分重要的；再一方面，根据它的连续谱可以确定序列的带宽（通常以谱的第一个零点作为序列的带宽）。 下面，以矩形脉冲构成的基带信号为例，通过几个有代表性的特例对功率谱密度公式的应用及意义做进一步的说明，其结果对后续问题的研究具有实用意义。
随机信号分析题目及答 案 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 )非正交 且线性相关。
功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。 功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。 谱是个很不严格的东西，常常指信号的Fourier变换，是一个时间平均（time average）概念功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。有两个重要区别：1。功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列）2。功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。热心网友回答提问者对于答案的评价：谢谢解答。 频谱分析（也称频率分析），是对动态信号在频率域内进行分析，分析的 结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变 量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密 度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。 功率谱是个什么概念？它有单位吗? 随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。 功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。 一是用相关函数的傅立叶变换来定义谱密度；二是用随机过程的有限时间傅立叶变换来定义谱密度；三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零，这样才能保证相关函数在时差趋向于无穷时衰减，所以lonelystar说的不全对，光靠相关函数解决不了许多问题，要求太严格了；对于第二种方式，虽然一个平稳随机过程在无限时间上不能进行傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换（FFT）估计谱密度的依据；第三种方式是根据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进行重构，在依靠正交性来建立的。 另外，对于非平稳随机过程，也有三种谱密度建立方法，由于字数限制，功率谱密度的单位
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章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。 其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在0dB附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB）。 4. Welch法估计及其MATLAB函数 Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计（CSD）。 MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计，而不需要自己编程。 （1）函数psd利用Welch法估计一个信号自功率谱密度，函数调用格式为： [Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])
关于matlab 中噪声功率谱密度与方差之间的关系的理解 1. 连续时间系统 高斯白噪声的定义为：如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。 故对于连续时间系统，理想的高斯白噪声的功率谱密度是一个常数，设为n0，而带宽是无限宽的，其功率为： 0*n ∞=∞ （1） 在n0不是为无穷小的情况下，理想的噪声功率Pn 是无限大的。 而实际当中，噪声带宽是有限宽的，只需要在我们所关心的频带范围内，噪声功率谱密度是个常数，则我们可认为其是高斯白噪声。设噪声单边功率谱密度为0n ，低通带宽为W ，则其噪声功率为： 0*2n n P W = （2） 如图1.1所示： o W -W 幅度 频率/HZ 0 2 n 图1.1 我们知道，高斯白噪声的分布为2 ~(0,)X N σ，则其功率为： 222()()()()n P E x D x E x D x σ==+== （3） 故对于低通系统有： 20/2 n W σ= （4） 而对于带通系统，如图1.2所示，有： 200*2*2n n P W n W σ=== （5）
W -W 幅度 频率/HZ 0 2 n 2. 离散时间系统 对于离散时间系统而言，带宽受到抽样速率fs 的限制。设WGN 一秒内抽取的一组数据样本为： 12[],,....fs x n x x x = 22([])0;([])([])E x n D x n E x n σ=== 2.1理论分析 由于时间为单个的离散点，故理想功率为0；但有下列定义：对于序列[]x n 的能量E 定义为序列各抽样值的平方和，则数据样本的能量为： 2221()*[()]*s f s s E x n f E x n f σ===∑ （6） 将功率定义为序列能量除以序列的时间，即 2*t s b E P f T σ==（单位：J/S ） （7） 式中，Tb 为序列时间，此处等于1S 。 如果功率单位采用W/symbol ，则有： 2/s t s P P f σ==（单位：J/symbol ） 2.2另一种理解 而实际当中，抽样点是一个时间段，认为1/s s T f =时间内的幅值就等于此抽样时刻的幅值，则单位抽样时间内的噪声能量为： 22***t s s s E E T f T σσ=== （6） 则噪声功率（单位：J/symbol ）为：
功率谱密度的三种matlab 实现方法 一：实验目的： (1)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。 二;实验内容： （1）简单说明三种方法的原理。 （2）用三种方法编写程序，在matlab 中实现。 （3）将计算结果表示成图形的形式，给出三种情况的功率谱图。 （4）比较三种方法的特性。 （5）写出自己的心得体会。 三：实验原理： 1.周期图法： 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列 )(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧ 序列。 )()(1)(1 m n x n x N m R N n N N x += ∑-=∧ (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 jwm M M m X jw x e m R e S ----=∧∧ ∑= )()(1) 1( 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M< 功率谱密度估计方法的MATLAB实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。 以下程序运行平台：Matlab R2015a（8.5.0.197613） 一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率100kHz N=1024; %数据长度N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title('信号') xlabel('时间');ylabel('幅度'); 实验报告 实验课程：数字信号处理实验开课时间：2020—2021 学年秋季学期 实验名称：随机信号功率谱分析实验时间： 2020年9月30日星期三 学院：物理与电子信息学院年级：大三班级：182 学号：1843202000234 姓名：武建璋 一、实验预习 实验目的要求深刻理解随机信号的特性，掌握随机信号功率谱估计的基本原理，灵活运用各种随机信号功率谱估计的基本方法。 实验 仪器 用具 装有Matlab的计算机一台 实验原理 功率谱估计是随机信号处理中的一个重要的研究和应用领域.功率谱估计基本上可以非参数估计的经典方法和参数估计的近代方法.典型功率谱估计是基于FFT 算法的非参数估计,对足够长的记录数据效果较好。 在工程实际中，经典功率谱估计法获得广泛应用的是修正期图发。该方法采取数据加窗处理再求平均的办法。通过求各段功率谱平均，最后得到功率谱计P（m），即： 式中：为窗口函数ω[k]的方差。K表示有重叠的分数段。 由于采用分段加窗求功率谱平均，有效地减少了方差和偏差，提高了估计质量，使修正周期图法在经典法中得到普遍应用。但在估计过程存在两个与实际不 符的假设，即 （1）利用有限的N个观察数据进行自相关估计，隐含着在已知N个数据之外的全部数据均为零的假设。 （2）假定数据是由N个观察数据以N为周期的周期性延拓。同时在计算过程中采用加窗处理，使得估计的方差和功率泄露较大，频率分辨率较低，不适用于短系列的谱分析和对微弱信号的检测。 近代谱估计是建立在随机信号参数模型的基础上，通过信号参数模型或预测误差滤波器（一步预测器）参数的估计，实现功率谱估计。由于既不需要加窗，又不需要对相关函数的估计进行如经典法那样的假设，从而减少公里泄露，提高了频谱分辨率。常用的参数模型有自回归（AR）模型、滑动平均（MA）模型、自回归滑动平均（ARMA）模型。其中AR模型是基本模型，求解AR模型的参数主要有L—D算法和Burg算法。 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 一、直接法： 直接法又称周期图法，它是把随机序列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太小，谱的分辨率又不好，因此需要改进。功率谱密度估计方法的MATLAB实现
(实验六 随机信号功率谱分析)
计算功率谱密度
相关主题
文本预览