当前位置:首页 >> 电力/水利 >>

随机信号的功率谱密度估计方法仿真


随机信号的功率谱密度估计方法仿真 目录
摘要 .......................................................................................................................................................... I ABSTRA

CT........................................................................................................................................... II 引言 ..........................................................................................................................................................1 1 绪论 ....................................................................................................................................................2 1.1 功率谱研究的发展过程 ..............................................................................................................2 1.2 功率谱估计方法提出 ..................................................................................................................3 1.3 功率谱估计应用及用途 ..............................................................................................................3 2 MATLAB 软件介绍 .............................................................................................................................5 2.1 MATLAB 产生的历史背景 ............................................................................................................5 2.2 基本功能 ......................................................................................................................................6 2.3 应用 ..............................................................................................................................................6 2.4 特点优势 ......................................................................................................................................7 3 随机信号的古典法谱估计 ................................................................................................................9 3.1 周期图法 ......................................................................................................................................9 3.2 相关法谱估计 (BT) ...............................................................................................................9 3.3 巴特利特(BARTLETT)平均周期图的方法 ................................................................................. 11 3.4 4 WELCH 法 .................................................................................................................................14

各种估计方法编程与分析 ............................................................................................................16 4.1 直接法 ........................................................................................................................................16 4.2 间接法 ........................................................................................................................................18 4.3 改进的直接法: ........................................................................................................................19 4.3.1 4.3.2 Bartlett 法 ...........................................................................................................................19 Welch 法 .............................................................................................................................20

4.4 各种方法的比较总结 ................................................................................................................23 5 结论 ................................................................................................................................................24

参考文献 ................................................................................................................................................25

摘要
对确定性信号进行傅里叶变换是在频率域分析研究的理论基础,但对于随机信 号来说,其傅里叶变换并不存在,因此转型研究它的功率谱。功率谱估计是数字信 号处理中的一个重要组成部分,其理论研究已经非常成熟,但是传统的仿真实现由 于计算量的庞大而往往令人望而却步。采用功能强大的 MATLAB 软件对常用的几 种线性法功率谱估计进行仿真,通过仿真结果讨论了这几种功率谱估计法的特点, 并对其比较分析,便于深刻理解各种方法的特点,从而在实际工作中做出合理的选 择。本文简要介绍了 MATLAB 仿真软件的基本功能、特点和优势,然后具体介绍 了周期图法、平均周期图法、窗函数法等几种经典的功率谱估计方法的原理,并对 各种方法进行仿真,通过结果分析其优缺点,对这几种方法进行质量评价,并提出 了改进方向。 关键词:自相关函数; 功率谱; MATLAB; 随机信号

I

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. 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 MATLAB 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, the average cycle diagram, the window function and so on, and simulate various methods we could get its advantages and disadvantages of several methods by the results, assess the quality of these methods, and rise the improved direction. Keyword:matlab; power spectrum estimation; sigle process

II

随机信号的功率谱密度估计方法仿真

引言
对信号和系统进行分析研究、处理有两种方法:一类是在时域进行,我们学习的 维娜一卡尔曼滤波和自适应滤波都属于这种方法;本论文是在频率与急性研究的另 一类方法。这两类方法都是信号处理的重要方法。功率谱是无限多个自相关函数的 函数,单观测数据只有有限个,只能得到有限个自相关函数。要无限个观测数根据 有限个样本数据,分析计算随机序列的真正功率谱,是求功率谱的中心问题,毫无 疑问,这是一个功率谱的估计问题。 功率谱估计是随机信号处理的重要内容,功率谱的估计方法有很多,主要有经 典估计,也称线性估计。经典谱估计有可以分成两种,一种是 BT 法,也叫间接法; 另一种是周期图法。1958 年,布莱克曼(R,Blackman)和杜基(J.Tukey)首先提 出 BT 法,并命名为布莱克曼-杜基谱估计器(简称 BT 普估计器) 。这种方法是先 按照有限个观测数据估计自相关函数,再计算功率谱。在 1965 年 FFT 未出现以前, BT 法一直是最常用的方法。1898 年,舒斯特(Schuster)利用傅里叶级数去拟合待 分析信号, 试图寻找信号中隐藏的周期性, 由此提出了周期图的术语, 但一直到 FFT 出现以后, 周期图法才受到人们的重视。 这种方法是直接对观测数据进行 DFT(FFT), 取模的平方,再除以 N 得到功率谱。比较两种方法,周期图法简单,不用估计自相 关函数,且可以用 FFT 进行计算。因此,周期图法得到了更广泛的应用。经典谱估 计致命的缺点是频率分辨率低,其原因是傅里叶变换域是无限大,而观察数据只能 得到有限个,观察不到的数据被认为是 0。如果只有 N 个观测数据,而对于 N 以外 的数据,信号仍有较强的相关性,这样估计出功率谱会出现很大的偏差。对于有限 个观测数据,相当于信号在时域中乘以矩形窗函数,在频域则使真正的功率谱卷积 一个 sinc 函数,而 sinc 函数不同于 ? 函数,它有主瓣和旁瓣,这样使卷积后的功率 谱不同于真正的功率谱。sinc 函数的主瓣引起功率谱向附近频域扩展,造成谱的模 糊,降低谱的分辨率;而旁瓣则会引起谱间干扰,信号强的功率谱旁瓣影响信号弱 的功率谱检测,严重时检测不出弱信号,或者把旁瓣误认为是信号,造成假信号。 为了对经典谱估计进行改进,可以采用各种不同的窗函数,但都是一增加主瓣宽度 来换取旁瓣的压低,因此谱分辨率低是经典谱估计的致命缺点。

1

随机信号的功率谱密度估计方法仿真

1

绪论

1.1 功率谱研究的发展过程
功率谱估计是数字信号处理的主要内容之一,主要研究信号在频域中的各种特 征,目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计 的发展过程做简要回顾: 英国科学家牛顿最早给出了“谱”的概念。后来,1822 年,法国工程师傅立叶提 出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理 论基础。 傅立叶级数提出后,首先在人们观测自然界中的周期现象时得到应用。19 世纪 末,Schuster 提出用傅立叶级数的幅度平方作为函数中功率的度量,并将其命名为 “周期图”(periodogram) 。这是经典谱估计的最早提法,这种提法至今仍然被沿用, 只不过现在是用快速傅立叶变换( FFT)来计算离散傅立叶变换(DFT) ,用 DFT 的幅度平方作为信号中功率的度量。 周期图较差的方差性能促使人们研究另外的分析方法。 1927 年, Yule 提出用线 性回归方程来模拟一个时间序列。Yule 的工作实际上成了现代谱估计中最重要的方 法——参数模型法谱估计的基础。 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 的快速计算方法。这些工作为现代谱估计的发展打下
2

随机信号的功率谱密度估计方法仿真

了良好的理论基础。1965 年,Cooley 和 Tukey 提出的 FFT 算法,也促进了谱估计 的迅速发展。

1.2 功率谱估计方法提出

在通信系统中,往往需要研究具有目中统计特性的随机信号。由于随机信号是 一类持续时间无限长,具有无限大能量的功率信号,它不满足傅里叶变换条件,而 且也不存在解析表达式,因此就不能够应用确定信号的频谱计算方法去分析随机信 号的频谱。然而,虽然随机信号的频谱不存在,但其相关函数确实确定的。如果随 机信号是平稳的,那么其相关函数的福利叶变换就是它的功率谱密度函数,简称功 率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱 密度函数做出估计。

1.3 功率谱估计应用及用途
功率谱估计有着极其广泛的应用,不仅在认识一个随机信号时,需要估计它的 功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所,要求预先 知道信号的功率谱密度(或自相关函数)。例如,在最佳线性过滤问题中,要设计一 个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。 根据信号与噪声的功率谱(或 ? xx (m) )才能设计出能够尽量不失真的重现信号,而把 噪声最大限度抑制的维纳滤波器。 常常利用功率谱估计来得到线性系统的参数估计。例如,当我们要了解某一系 统的幅频特性 H (? ) 时,可用一白色噪声 ? ?n ? 通过该系统。再从该系统的输出样本 y(n) 估计功率谱密度 Pyy (? ) 。由于白色噪声的 PSD( 用 P?? (?) 表示 ) 为一常数即
2 ,于是有: P?? (?) ? ? ?
2 Pyy (? ) ? ? ? H (? ) 2

(1-1)

故通过估计输出信号的 PSD,可以估计出系统的频率特性 H (?) (模特性)。从 宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这 要求功率谱估计有足够好的频率的分辨率,否则就不一定能够清楚地检测出来。所 谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频 率间隙(距)。提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向 功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关 系,实际用途有滤波,信号识别(分析出信号的频率) ,信号分离,系统辨识等。 谱估计技术是现代信号处理的一个重要部分,还包括空间谱估计,高阶谱估计等。 维纳滤波、卡尔曼滤波,可用于自适应滤波,信号波形预测等(火控系统中的飞机
3

随机信号的功率谱密度估计方法仿真

航迹预判) 。

4

随机信号的功率谱密度估计方法仿真


2.1

matlab 软件介绍
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 敏锐地觉察到 MATLAB 在工程领域的广阔前景。同年他和 Cleve Moler SteveBangert 一起,用 C 语言开发了第二代专业版.这一代的 MATLAB 语 言同时具备了数值计算和数据图示化的功能。 1984 年 CleveMoler 和 JohnLittle 成立 了 MathWorks 公司,正式把 MATLAB 推向市场,并继续进行 MATLAB 的研究和 开发。 在当今 30 多个数学类科技应用软件中,就软件数学处理的原始内核而言,可 分为两大类:一类是数值计算型软件,如 MATLAB、Xmath、Gauss 等,这类软件长 于数值计算,对处理大批数据效率高。另一类是数学分析型软件,Mathematica、 Mapl 等,这类软件以符号计算见长,能给出解析解和任意精确解,其缺点是处理大 量数据时效率较低。Math Works 公司顺应多功能需求之潮流,在其卓越数值计算和 图示能力的基础上,又率先在专业水平上开拓了其符号计算。文字处理,可视化建 模和实时控制能力, 开发了适合多学科多部门要求的新一代科技应用软件 MATLAB. 经过多年的国际竞争,MATLAB 以经占据了数值软件市场的主导地位。 在 MATLAB 进入市场前,国际上的许多软件包都是直接以 FORTRANC 语言 等编程语言开发的。这种软件的缺点是使用面窄,接口简陋,程序结构不开放以及 没有标准的基库,很难适应各学科的最新发展,因而很难推广。MATLA 的出现, 为各国科学家开发学科软件提供了新的基础。 在 MATLAB 问世不久的 80 年代中期, 原先控制领域里的一些软件包纷纷被淘汰或在 MATLAB 上重建。 MathWorks 公司 1993 年推出了 MATLAB4.0 版,1995 年推出 4.2C 版 1997 年
5

随机信号的功率谱密度估计方法仿真

推出 5.0 版。1999 年推出 5.3 版。MATLAB5 X 版较 MATLAB4,X 版无论是界面 还是内容都有长足的进展, 其帮助信息采用超文本格式和 PDF 格式, 在 Netscape3.0 或 IE4.0 及以上版本,AcrobatReader 中可以方便地浏览。 时至今日,经过 MathWorks 公司的不断完善,MATLAB 已经发展成为适合多 学科,多种工作平台的功能强大大大型软件。在国外,MATLAB 已经经受了多年考 验。在欧美等高校,MATLAB 已经成为线性代数,自动控制理论,数理统计,数字 信号处理,时间序列分析,动态系统仿真等高级课程的基本教学工具;成为攻读学 位的大学生,硕士生,博士生必须掌握的基本技能。在设计研究单位和工业部门, MATLAB 被广泛用于科学研究和解决各种具体问题。在国内,特别是工程界, MATLAB 一定会盛行起来。可以说,无论你从事工程方面的哪个学科,都能在 MATLAB 里找到合适的功能。

2.2 基本功能
MATLAB 和 Mathematica、Maple 并称为三大数学软件。它在数学类科技应用 软件中在数值计算方面首屈一指。MATLAB 可以进行矩阵运算、绘制函数和数据、 实现算法、创建用户界面、连 接其他编程语言的程序等,主要应用于工程计算、 控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB 的基本数据单位是矩阵, 它的指令表达式与数学、 工程中常用的形式 十分相似,故用 MATLAB 来解算问题要比用 C,FORTRAN 等语言完成相同的事 情简捷得多,并且 mathwork 也吸收了像 Maple 等软件的优点,使 MATLAB 成为一 个强大的数学软件。 在新的版本中也加入了对 C, FORTRAN, C++ , JAVA 的支持。 可以直接调用,用户也可以将自己编写的实用程序导入到 MATLAB 函数库中方便自 己以后调用,此外许多的 MATLAB 爱好者都编写了一些经典的程序,用户可以直 接进行下载就可以用。

2.3 应用
MATLAB 产品族可以用来进行以下各种工作: ● 数值分析。 ● 数值和符号计算。 ● 工程与科学绘图。 ● 控制系统的设计与仿真。 ● 数字图像处理 技术。 ● 数字信号处理 技术。
6

随机信号的功率谱密度估计方法仿真

● 通讯系统设计与仿真。 ● 财务与金融工程。 在欧美大学里,诸如应用代数、数理统计、自动控制、数字信号处理、模拟与 数字通信、时间序列分析、动态系统仿真等课程的教科书都把 MATLAB 作为内容。 这几乎成了九十年代教科书与旧版书籍的区别性标志。 在国际学术界,MATLAB 已经被确认为准确、可靠的科学计算标准软件。在许 多国际一流学术刊物上, (尤其是信息科学刊物) ,都可以看到 MATLAB 的应用。 如今 MATLAB 已应用于下列领域: 数值和符号计算、工程与科学绘图、数值分析、建模和仿真、控制系统的设计 与仿真、数字图像处理、数字信号处理、通讯系统设计与仿真、财务与金融工程。 具体表现为:自动控制、航天工程、汽车工业、生物医学工程、语音处理、图像处 理、雷达工程、信号分析、计算机技术等各行各业中。

2.4 特点优势
特点: ●此高级语言可用于技术计算。 ●此开发环境可对代码、文件和数据进行管理。 ●交互式工具可以按迭代的方式探查,设计及求解问题。 ●数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数值积分等。 ●二维和三维图形函数可用于可视化数据。 ●各种工具可用于构建自定义的图形用户界面。 ●各种函数可将基于 MATLAB 的算法与外部应用程序和语言(如 C、C++、 Fortran、Java、COM 以及 Microsoft Excel)集成。 ●不支持大写输入,内核仅仅支持小写。 优势: (1)友好的工作平台和编程环境。 (2)简单易用的程序语言。 (3)强大的科学计算机数据处理能力。 (4)出色的图形处理功能。
7

随机信号的功率谱密度估计方法仿真

(5)应用广泛的模块集合工具箱。 (6)实用的程序接口和发布平台。 (7)应用软件开发(包括用户界面) 。

8

随机信号的功率谱密度估计方法仿真

3

随机信号的古典法谱估计

3.1 周期图法
周期图法又称直接法。它是从随机信号 x(n)中截取 N 长的一段,把它视为能量 有限 x(n)真实功率谱 Sx (e jw ) 的估计 Sx (e jw ) 的抽样. 周期图这一概念早在 1899 年就提出了,但由于点数N一般比较大,该方法的 计算量过大而在当时无法使用。只是 1965 年 FFT 出现后,此法才变成谱估计的一 个常用方法。周期图法包含了下列二条假设: 1.认为随机序列是广义平稳且各态遍历的,可以用其一个样本 x(n)中的一段

xN (n) 来估计该随机序列的功率谱。这当然必然带来误差。
2.由于对 xN (n) 采用 DFT,就默认 xN (n) 在时域是周期的,以及 xN (k ) 在频域 是周期的。这种方法把随机序列样本 x(n)看成是截得一段 xN (n) 的周期延拓,这也 就是周期图法这个名字的来历。与相关法相比,相关法在求相关函数 Rx (m) 时将

xN (n) 以外是数据全都看成零,因此相关法认为除 xN (n) 外 x(n)是全零序列,这种处
理方法显然与周期图法不一样。 但是,当相关法被引入基于 FFT 的快速相关后,相关法和周期图法开始融合。 比我们发现:如果相关法中 M=N,不加延迟窗,那么就和充(N-1)个零的周期图 法一样了。简单地可以这样说:周期图法是 M=N 时相关法的特例。因此相关法和 周期图法可结合使用。

3.2 相关法谱估计 (BT)
这种方法以相关函数为媒介来计算功率谱,所以又叫间接法。它是 1958 年由 Blackman 和 Tukey 提出。这种方法的具体步骤是: 第一步:从无限长随机序列 x(n)中截取长度 N 的有限长序列列 xN (n) 第二步:由 N 长序列 xN (n) 求(2M-1)点的自相关函数 R x (m) 序列。即
R x (m) ?
?
?

1 N ?1 ? xN (n) xN (n ? m) N n ?0

(3-1)

这里,m=-(M-1)…,-1,0,1…,M-1,M N, R x (m) 是双边序列,但是由自相关函 数的偶对称性式,只要求出 m=0,。 。 。,M-1 的傅里叶变换,另一半也就知道了。 第三步:由相关函数的傅式变换求功率谱。即

S x (e jw ) ?

?

M ?1 ? X m ? ? ( M ?1)

?R

(m)e ? jwm
9

(3-2)

随机信号的功率谱密度估计方法仿真

以上过程中经历了两次截断,一次是将 x(n)截成 N 长,称为加数据窗,一次是 将想 x(n)截成(2M-1)长,称为加延迟窗。因此所得的功率谱仅是近似值,也叫谱 估计,式中的 Sx (e jw ) 代表估值。一般取 M<<N,因为只有当 M 较小时,序列傅式 变换的点数才较小,功率谱的计算量才不至于大到难以实现,而且谱估计质量也较 好。因此,在 FFT 问世之前,相关法是最常用的谱估计方法。 当 FFT 问世后,情况有所变化。因为截断后的 xN (n) 可视作能量信号,由相关 卷积定理可得
R x ( m) ?
?

1 [ x N (m) ? x N (?m)] N

(3-3)

这就将相关化为线性卷积,而线性卷积又可以用快速卷积来实现。我们可对上 式两边取(2N-1)点 DFT,则有
R x ( m) ?
?

1 1 2 [ x2 N ?1 (k ) ? x2 N ?1 ( K ) ? X 2 N ?1 ( K ) N N
?

(3-4)

于是将时域卷积变为频域乘积,用快速相关求 R x (m) 的完整方案如下:
1. 2.

对 N 长 xN (n) 的充(N-1)个零,成为(2N-1)长的。 求(2N-1)点的 FFT,得 X 2 N ?1 ( K ) ? 求
2 N ?2 N ?0

?x

2 N ?1

(n)W2?Nmk ?1 。

3.

1 2 X 2 N ?1 ( K ) 。由 DFT 性质,x2 N ?1 (n) 是纯实的,x2 N ?1 (k ) 满足共轭偶对称, N 1 2 X 2 N ?1 ( K ) 一定是实偶的,且以(2N-1)为周期。 而 N

4. 求(2N-1)点的 IFFT:

R x (m) ?
这里

?

N ?1 1 1 2 X 2 N ?1 ( K ) W2?Nmk ? ?1 2 N ? 1 k ??( N ?1) N

(3-5)

1 2 X 2 N ?1 ( K ) 是实偶的, m=-(N-1)...0...N-1。 本来 IFFT 求和范围是 0 至 2N-2, N 1 2 X 2 N ?1 ( K ) 的实偶性与周期性,求和范围改为-(N-1)至(N-1)不影响计 由于 N

算结果。同理可将 m 的范围改为-(N-1)至(N-1)。 上述的快速相关中,充零的目的是为了能用圆周卷积代替线性卷积,以便进一 步采用快速卷积算法。快速相关输出是-(N-1)至(N-1)的 2N-1 点,加 WM (m) 窗 后截取的是-(M-1)至(M-1)的频段,最后作(2M-1)点 FFT,得 S x (k ) 。我们 注意到:如果数据点数与自相关序列点数相同即 M=N,则(2N-1)点的 IFFT 后紧
10
?

随机信号的功率谱密度估计方法仿真

跟一个(2N-1)点的 FFT,利用 R x (m) 的对称性,FT 运算框的计算式变为

?

S X (K ) ?

m ? ? ( N ?1)

?R

N ?1 ?

X

(m) W2?Nmk ?1 X 2 N ?1

(3-6)

由于 N=M 并假设窗形状是矩形的,第二次 WM ?m? 的截断就不需要了。比较式
? FFT [R x (m)] , R x (m) ? IFFT [ (3-5)和式(3-6), S( x k)
?
?

1 2 X 2 N ?1 ( K ) ] 正反傅氏 N

变换可以抵消,直接得

? S( x k)

1 2 X 2 N ?1 ( K ) N

(3-7)

为了实行基 2FFT, 也可将 (2N-1) 点换成 2N 点, 这样做不影响结果的正确性。

3.3 巴特利特(Bartlett)平均周期图的方法
首先让我们来看一下为什么周期图经过某种平均 ( 或平滑 ) 后会使它的方差当 N ? ? 时趋于零,达到一致估计的目的。如果 x1 , x2 , ?, x L 是不相关的随机变量, 每一个具有期望值 ? , 方差 ? 2 , 则可以证明它们的数学平均 x ? ( x1 ? x2 ? ...? xl ) / L 的 期望值等于 ? ,数学平均的方差等于 ? 2 / L ,即:
E?x ? ? 1 1 E?x1 ? x2 ? ? ? x L ? ? ? L? ? ? L L

Var?x ? ? E ( x ? E( x ))2 ? E[ x 2 ] ? ( E?x ?) 2
? 1 E ( x1 ? x 2 ? ? ? x L ) 2 ? ? 2 2 L

?

?

?

?

1 ? 2 L

? L L ? 2 2 2 E x ? x ? ? ? x ? E xi x j ? 1 ?? 2 L j ? 1 i ? 1 ? i? j ?

?

?

?

?

? ? 2 ??? ? ?

?? E xi x j ? ? E x j ? ? E?xi ?
j ?1 i ?1 i? j j ?1 i ?1 i? j

L

?

?

L

? ?

L

? L? ( L ?1)? ? L2 ? 2 ? L? 2

所以

Var?x ? ?

1 ?L 1 ?L ? 2 2 2 2? 2 E x ? L ? ? L ? ? ? ? E xi2 ? L? 2 ? i ? 2 ?? 2 ?? L ? i ?1 L ? i ?1 ? ?
11

? ?

? ?

随机信号的功率谱密度估计方法仿真

?

1 2 2 E x12 ? ( E[ x1 ]) 2 ? E x2 ? ( E[ x2 ]) 2 ? ? ? E x L ? ( E[ x L ]) 2 L2

?? ? ?

? ?? ?

?

?? ?

??
(3-8)

?

1 ?2 2 L ? ? L L2

由 (3-8) 可 见 , L 个 平 均 的 方 差 比 每 个 随 机 变 量 的 单 独 方 差 小 L 倍 。 当

L ? ? Var?x ? ? 0 ,可达到一致谱估计的目的。因而降低估计量的方差的一种有效
方法是将若干个独立估计值进行平均。把这种方法应用于谱估计通常归功于 Bartlett。 Bartlett 平均周期图的方法是将序列 x(n) (0 ? n ? N ? 1) 分段求周期图再平均。 设将 x(n) 分成 L 段,每段有 M 个样本,因而 N ? LM ,第 i 段样本序列可写成

x i (n) ? x(n ? iM ? M )
第 i 段的周期图为

0 ? n ? M ? 1, 1 ? i ? L

1 I (? ) ? M
i M

M ?1 n ?0

?x

2 j

( n )e

? j?n

i 如果 m ? M , ? xx (m) 很小, 则可假定各段的周期图 I M 对功率 (? ) 是互相独立的。

谱密度的概念的讨论,谱估计可定义为 L 段周期图的平均,即
L ? (? ) ? 1 ? I i (? ) P xx M L i ?1

(3-9)

于是它的期望值为
L ? (? ) ? 1 ? E I i (? ) ? E I i (? ) EP xx M M L i ?1 ? (? ) ? 1 ? P (? )W (e j (? ?? ) )d? EP xx xx B 2? ???

? ?

? ?

?

? ?

?

1 ? 2?M

? sin ??? ? ? ? M 2 ? ? d? ??? Pxx (? )? ? sin ??? ? ? ? / 2 ?
?
2

(3-10)

这里 M ? N / L ,因此 Bartlett 估计的期望值是真实谱 Pxx (? ) 与三角窗函数的卷 积。由于三角窗函数不等于 ? 函数,所以 Bartlett 估计也是有偏估计即 Bias ? 0 ,但 当 N ? ? 时, Bias ? 0 。 由于我们假定各段周期图是相互独立的,所以可按式(3-8)得到下式:
? (? ) ? 1 Var ?I (? )? Var P xx M L
12

?

?

随机信号的功率谱密度估计方法仿真
2 ? ?s i n 1 2 ?(M ) ? ? ? Pxx (? ) ?1 ? ? ? ? L ?? ? ? ? Msin ? ?

(3-11)

? (?) ? 0 。因 ? (?) 是下降的,当 L ? ? 时,Var P 由此可见,随着 L 的增加 Var P xx xx

?

?

?

?

此 Bartlett 估计是一致估计。
? (?) 与式(3-1)的 E?I (?)? 可见在二种情况的估计量的期望值 比较式(3-10)的 E P N xx

?

?

都是真值 Pxx (? ) 与窗口函数 WB (e j? ) 的卷积形式,但后者将前者 WB 中 N 改为 M,
M ? N/L? N。 因而使 WB (e j? ) 主瓣的宽度增大。 由于主瓣的宽度愈窄愈接近 ? 函数,

偏 倚 愈 小 。 今 式 (3-10) 中 WB (e j? ) 的 主 瓣 宽 度 大 于 后 式 中 的 主 瓣 宽 度 , 因 而
? (?) ? Bias?I (?)?,而主瓣愈宽显然使分辨率愈差。因此 Bias 可用来说明谱的 Bias P xx N

?

?

分辨率,Bias 愈大说明谱分辨率愈差。一个固定的记录长度 N,周期图分段的数目 L 愈大将使方差愈小,但 M 也愈小,因而使 Bias 愈大,谱分辨率变得愈差。因此 Bartlett 方法中 Bias 或谱分辨率和估计量的方差间是有互换关系的。M 和 N 的选择 一般是由对所研究的信号的预先了解来指导的。 例如, 如果我们知道谱有一个窄峰, 同时如果分辨出这个峰是重要的,那么我们必须选择 M 足够大。又从方差的表达 式我们可以确定谱估计的可接受的方差所要求的记录长度 N=(LM)。 由此可见 Bartlett 法使谱估计的方差减小是用增加 Bias 以及降低谱分辨率的代 价换来的。实际上,当 N 一定时,Bias 与 Var 的互换性是谱估计的一个固有特性。 [例如] 为了说明经平均后的周期图作为功率谱估计的实际效果, 设有一零均值 高斯分布的随机过程,其功率谱密度为
Pxx ( z ) ? 1 (1 ? az )(1 ? az)
?1

a ? 1 ? 0.1

这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其
H ( z) ? 1/(1 ? az?1 ) 的滤波器形成的。为了简单设选用 N ? 8 的矩形窗函数。其周期图的

期望值(用 E[ I 8 ( w)]表示)与真值 Pxx (w) 均表示在图 3.4 中。它说明 N ? 8 的周期图可以 得到的 Bias 的情况。图 3-2 表示 M=8 分 4 段与 16 段二种情况平均后的周期图。显 然 L=4 的方差比 L=16 的大。 将 L=16 的曲线与图 5.4 的曲线比较可见在这种估计中 误差的大部分起因于 Bias 而不是方差(因为二种情况均为 M
? 8 ,Bias

相仿,误差也

相仿)。M=16,L=2 及 8 的周期图表示在图 3-3。图 3-3 中 L 不同造成的影响也是明 显的。但是这二个曲线的起伏都很大,因此有理由认为误差主要起因于方差。比较
M ? 8与

M=16 的周期图可见, M 愈大(即§3.3 中的 N 愈大),将使周期图的起伏愈

增快的结论,在这里也同样成立。 比较图 3-2 与图 3-3 发现在这个例子中最好的选择是应用 L=16, M=8 的估计而
13

随机信号的功率谱密度估计方法仿真

不是 L=8、M=16 的估计,即宁可减小方差,牺牲 Bias。在实际中,当然功率谱密 度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预 先知道的。通过改变 M 和 L 以及利用预先知道的情况,通常可以很好地进行选择。 平均周期图的方法特别适合于应用 FFT 算法。因此在 FFT 出现以后这个方法 比下面将要讨论的利用窗口函数的处理法用得更多。 而在 FFT 出现以前主要是用窗 口函数处理法平滑周期图。

图3

图 3-1

Px (?) 与 E[ I 8 (? )]的特性

图 3-2

平滑后的周期图 (每段取 8 个数据)

图 3-3

平均后的周期图(每段取 16 个数据)

3.4

Welch 法
Welch 提出了对 Bartlett 法的修正使更适合于用 FFT 进行计算。他主要提出二

方面的修正,其一是选择适当的窗函数 ? (n) ,并在周期图计算前直接加进法,这样 得到的每一段的周期图为

1 I (? ) ? MU
(i ) M

? x (n)? (n)e
i n ?0

N ?1

2 ? j?n

(3-12)

这里 U ?

1 M

?W
n ?0

N ?1

2

(n) 为归一化因子,而 Bartlett 法每段的周期图为
14

随机信号的功率谱密度估计方法仿真

1 I (? ) ? M
(i ) M

? x ( n) e
i n ?0

N ?1

2 ? j?n

(3-13)

这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段 时,可使各段之间有重迭,这样将会使方差减小(当 N 与 M 一定时)。他建议可以重 迭 50%。

15

随机信号的功率谱密度估计方法仿真



各种估计方法编程与分析

4.1 直接法
直接法又称周期图法,它是把随机序列 x(n)的 N 个观测数据视为一能量有限的 序列,直接计算 x(n)的离散傅立叶变换,得 X(k),然后再取其幅值的平方,并除以 N,作为序列 x(n)真实功率谱的估计。 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));

图 4-1

周期图法中 N=1000 功率谱

16

随机信号的功率谱密度估计方法仿真

图 4-2

周期图中 N=100000 功率谱

图 4-3

周期图当 N=100 功率谱

随着采样点数的增加,该股集是渐进无骗的。从图中可以看出,采用周期突发 估计得出的功率谱很不平滑,相应的估计协方差比较大。而且采用增加采样点的办 法也不能吃周期图变得更加平滑,这是周期图法的缺点。周期图法得出的估计谱方 差特性不好: 当数据长度 N 太大时, 扑线的起伏加剧; N 太小时谱的分辨率又不好。 对其改进的主要方法有二种,即平均和平滑,平均就是将截取的数据段 x N (n) 再分 成 L 个小段,分别计算功率谱后取功率谱的平均,这种方法使估计的方差减少,但 偏差加大,分辨率下降。平滑是用一个适当的窗函数 W (e jw ) 与算出的功率谱
S X (e jw ) 进行卷积,使谱线平滑。这种方法得出的谱估计是无偏的,方差也小,但
?

分辨率下降。

17

随机信号的功率谱密度估计方法仿真

4.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);

图 4-4

直接发功率谱图

当 M=N-1 时,BT 法与周期图法估计出的功率谱是一样的;当 M < N-1 时,BT 法
18

随机信号的功率谱密度估计方法仿真

的偏差大于周期图法,在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的 方差;分辨率比周期图法低,与窗函数的选择有关。
?

BT 法的缺点在于当 M →N 时, R ( m ) 的方差很大,使谱估计质量下降;由 R ( m ) 得到的 S ( w) 不一定为正值,从而可能失去功率谱的物理意义。
?

?

4.3 改进的直接法:
对于直接法的功率谱估计, 当数据长度 N 太大时, 谱曲线起伏加剧, 若 N 太小, 谱的分辨率又不好,因此需要改进。 4.3.1 Bartlett 法

对周期图法改进的思想是将信号分段进行估计,然后再将这些估计结果进行平 均,从而减小估计的协方差,是功率谱图变得比直接法更平滑。增加分段数可以进 一步减低估计的协方差,然而若每段中的数据点数太少,就会使估计的频率分辨率 下降很多。从样本信号序列总点数一定的条件下,可以采用使分段相互重叠的方法 来增加分段数,从而保持每段信号点数不变,这样就在保证频率分辨率的前提下进 一步降低估计协方差。而一般的 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;
19

随机信号的功率谱密度估计方法仿真

figure(2) plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

图 4-5 4.3.2 Welch 法

bartlett 法估计功率谱图

现在比较常用的改进方法是 Welch 法,又叫加权交叠平均法,简记为 WOSA 法,这种方法以加窗(加权)求取平滑,以分段重叠求得平均,因此集平均与平滑 的优点于一体,同时也不可避免带有两者的缺点,因此归根到底是一种折中。其主 要步骤是: (1)将 N 长的数据段分成 L 个小段,每小段 M 点,相邻小段间交叠 M/2 点 (即 2:1 分段) 。因为 L(M/2)+M/2=N,所以段数
L? N ?M /2 M /2
(4-1)

(2)对各小段加同样的平滑窗
M ?1 n ?0

后起傅氏变换

X i (e ) ? ? xi (n)ax(n)e ? jwn , i ? 1,...L,
jw

(4-2)

(3)用下式求各小段功率谱的平均
S i (e jw ) ?
? 2 1 l 1 X i (e jw ) ? L i ?1 MU

(4-3)

这里,

代表窗函数平均功率,MU 是 M 长窗函数的能量。
20

随机信号的功率谱密度估计方法仿真

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(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); pause; figure(2) plot(f,plot_Pxx1); pause; figure(3) plot(f,plot_Pxx2);

21

随机信号的功率谱密度估计方法仿真

图 4-6

矩形窗

图 4-7

海明窗

图 4-8

blackman 窗

Welch 法对 Bartlett 法进行了两方面的修正,一是选择适当的窗函数 w(n),并 再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非 负。二是在分段时,可使各段之间有重叠,这样会使方差减小。 不同窗函数的 Welch 谱估计在选择窗函数时,一般有如下要求: (1)窗口宽度 M 要远小于样本序列长度 N ,以排除不可靠的自相关值; (2)当平稳信号为实过程时,为保证平滑周期图和真是功率谱也是实偶函数,平 滑窗函数必须是实偶对称的; (3)平滑窗函数应当在 m=0 出游峰值,并且 m 随绝对值增加而单调下降,使可 靠的自相关值有较大的权值;
22

随机信号的功率谱密度估计方法仿真

(4)功率谱是频率的非负函数且周期图是非负的,因而要求窗函数的 Fourier 变 换是非负的。 仍以上述平稳随机信号为例,采样频率、采样点数、FFT 点数、窗长度及重叠数 据不变,窗函数采用矩形窗、Blackman 窗、 Hamming 窗,仿真结果如图所示。 由仿真结果可以看出:使用不同的窗函数谱估计的质量是不一样的,矩形窗的 主瓣较窄,分辨率较好,但方差较大,噪声水平较高;而 Blackman 窗和 Hamming 窗的主瓣较宽,分辨率较低,但方差较小,噪声水平较低。因此,在进行谱分析时 选择何种窗函数,要视具体情况而定。如果强调高分辨率,能精确读出主瓣频率, 而不关心幅度的精度,例如测量震动物体的自震频率,可以选用主瓣宽度比较窄的 矩形窗;对受到强干扰的窄带信号若干扰靠近信号,则可选用旁瓣幅度较小的窗函 数,若离开通带较远,则可选用渐近线衰减速度比较快的窗函数。总之,要针对不 同的信号和不同的处理目的来选择合适的窗函数,这样才能得到良好的效果。

4.4 各种方法的比较总结
1) 周期图法和 BT 法的特点是离散性大,曲线粗糙,方差较大,但是分辨率较 高; (2) Bartlett 法和 Welch 法的收敛性较好,曲线平滑,方差较小,但是功率谱主 瓣较宽,分辨率低,这是由于对随机序列加窗截断所引起的 Gibbs 效应造成的; (3) 与 Bartlett 法相比,Welch 法的估计曲线比较粗糙,但是分辨率较好,原因 是 Welch 法中对数据进行截断时加的是 Hanning 窗,而在 Bartlett 法中使用的是矩 形窗,相对于矩形窗,Hanning 窗的主瓣包含更多的能量,因而使功率谱的主瓣较 窄,分辨率较高。 由上述理论分析及仿真实验可知,Welch 法采用加窗交叠求功率谱,可以有效 减小方差和偏差,一般情况下能接近一致估计的要求,因而得到广泛应用。同时还 可以发现,对信号加不同的窗函数,谱估计的质量是不同的。

23

随机信号的功率谱密度估计方法仿真



结论
在经典谱估计中, 无论是周期图法还是其改进的方法, 都存在着频率分辨率低、

方差性能不好的问题,原因是谱估计时需要对数据加窗截断,用有限个数据或其自 相关函数来估计无限个数据的功率谱,这其实是假定了窗以外的数据或自相关函数 全为零,这种假定是不符合实际的,正是由于这些不符合实际的假设造成了经典谱 估计分辨率较差。另外,经典谱估计的功率谱定义中既无求均值运算又无求极限运 算,因而使得谱估计的方差性能较差,当数据很短时,这个问题更为突出。如何选 取最佳窗函数、提高频谱分辨率,如何在短数据情况下提高信号谱估计质量,还需 要进一步研究。

参考文献
(1)美 John G.Proakis Dimitris G.Manolakis.数字信号处理 原理 算法与应用(第三版) (2)胡广书. 数字信号处理:理论、算法与实现 (第二版):出版日期: 2003-08-12 (3)英 艾费科,.Digital Signal Processing:A Practical Approach,Second Edition 第二版 (4) 苏金明.MATLAB 7.0 实用指南 [J].2004-11-1 (5) 魏巍.MATLAB 信息工程工具箱技术手册.2004-01-1. (6) 飞思科技产品研发中心 . Matlab7 辅助信号处理技术与应用 [M] . 北京 : 电子工业出版 社 ,2005. (7) John Wiley & Sons .Introduction to Discrete-Time Signal Processing, 1976. (8)J. Makhoul. Linear Prediction: A Tutorial Review Proc. IEEE, Vol. 63, No. 4, April 1975, pp. 561-580 (9)丁玉美.数字信号处理 时域离散随机信号处理 [M].西安电子科技大学出版社,2002. (10)姚武川,姚天任.1 经典谱估计方法的 Matlab 分析 [J].华中理工大学学报,2000 (4).

24


相关文章:
实验2 离散随机信号的计算机仿真
实验2 离散随机信号的计算机仿真(验证性实验) 一、实验目的(1)掌握指定分布随机信号的仿真方法; (2)掌握随机信号自相关函数和功率谱密度的分析方法; 二、实验步骤...
实验四 窄带随机信号的仿真与分析
实验四 窄带随机信号的仿真与分析_信息与通信_工程科技_专业资料。窄带随机信号的...信号的均值,均方值、方差、概率密度、频谱及功率谱密度、相关 函数,用图示法...
基于Burg算法的AR模型功率谱估计简介
然后根据 2 S x (e jw ) ??w H (e jw ) 2 即可求得该随机信号的功率谱密度。 3 实验仿真分析假定信号序列为两个正弦信号与白噪声的叠加,其功率谱估计...
随机信号的谱估计方法及实现+孟波+20091342083
实验三、随机信号的功率... 18页 免费 随机信号的功率谱密度估... 26页 1...采 用功能强大的 MATLAB 软件对常用的几种线性法功率谱估计进行仿真,通过对 于...
功率谱估计性能分析及Matlab仿真
利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计谱估计...4096 2 图 2 周期图法 N ? 128 说明: (1) 本报告仿真中所采用的用于...
二进制基带信号功率谱密度的研究
? 4. 数字基带信号功率谱密度仿真与分析 4. 1 随机二进制序列的产生 编写 bianry 函数,用以产生指定长度的随机二进制序列: function[out]= binary(nsize) r=...
随机信号分析仿真
(t ) 为一个具有零均值的平稳随机过程,其功率谱密度均匀分布在()整个频率区间...输出信号包络检测 本文相关仿真所使用的软件为 matlab,以滤波器 1 为例,仿真...
实验三 随机信号的功率谱估计
随机信号的功率谱估计 一、 实验目的(1)了解估计功率谱密度的几种方法; (2)...但对于用 matlab 进行仿真估计,也十分不熟悉。 于是,我们小组成员找出课本复习了...
随机信号处理基础 matlab仿真
随机信号处理基础 matlab仿真_数学_自然科学_专业资料。南京理工大学随机信号处理基础...(t ) 为确知信号, n(t ) 为均值为零的平稳白噪声,其功率谱密度为 No ...
随机信号及其自相关函数和功率谱密度的MATLAB实现
实现 摘要: 摘要: 学习用 rand 和 randn 函数产生白噪声序列; 学习用 MATLAB 语 言产生随机信号;学习用 MATLAB 语言估计随机信号的自相关函数 和功率谱密度。利...
更多相关标签:
随机信号功率谱密度 | 随机信号的功率谱密度 | 随机信号的仿真与分析 | 功率谱密度估计 | 随机振动功率谱密度 | 随机振动的功率谱密度 | 随机过程的功率谱密度 | 随机序列的功率谱密度 |