周期图法的功率谱估计_黄英

[作者简介] 黄英（1965—），女，副教授，研究方向为嵌入 式系统应用开发。
Power spectrum Estimation Based on Periodogram HUANG Ying
(Suzhou Institute of Trade & Commerce,Suzhou,Jiangsu 215009,China)
[作者简介] 黄吉秀， 女， 副教授， 四川财经职业学院副院长， 研究方向为经济学。
（上接第 2 页）
6
结论
参考文献 [1] 丁玉美.数字信号处理——时域离散随机信号处理[M].西安:西安电 子科技大学出版社,2002 [2] 陈怀琛,吴大正.MATLAB 及在电子信息课程中的应用 （第二版） [M]. 北京:电子工业出版社,2004 [3] 李富强,万红.MATLAB6 的语谱图显示与分析[J].微计算机信息,2005, 10(3):1-5 [4] 魏鑫,张平.周期图法功率谱估计中的窗函数分析[J].现代电子技术, 2005,(3)
参考文献 [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)
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；
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
基本的周期图法
中国西部科技
2013 年 09 月第 12 卷第 09 期总第 290 期
1
周期图法的功率谱估计
黄 英
（苏州经贸职业技术学院，江苏 苏州 215009） 摘 要：功率谱估计是数字信号处理的重要研究内容之一。对目前用得较多且最具代表性的周期图法进行研究，并针对它谱 分辨率比较低等缺点，利用 Matlab 实现了几种改进的周期图法。 关键词：谱估计；周期图；Matlab DOI：10.3969/j.issn.1671-6396.2013.09.001
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
图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） ；
论
著
与创新研发的投入，提高科技与创新投入的效益，优先保证 关键技术研发创新和重大产业化项目建设等的要素需求， 促 进经济发展从外延扩张型向内涵提高型转变。同时，加强企 业监管和企业服务体系建设，严格限制处于产业链低端环 节、特色优势不明显、高消耗、低效益的行业或企业扩张产 能，支持优势企业做大做强，提升核心竞争力。
5
多窗口的周期图法
由于普通的周期图功率谱估计只利用单一窗口， 因此在 序列始端和末端均会丢失相关信息。 多窗口法利用多个正交 窗口获得各自独立的近似功率谱估计，然后综合这些估计， 最终得到平稳信号序列的功率谱估计。 图 4 给出了多窗口周期图法得到的功率谱。
4
Welch 法
P.O.韦尔奇提出一种把加窗处理与平均处理结合起来的 方法。先把分段的数据乘以窗函数（进行加窗处理） ，分别 计算其周期图，然后进行平均，得到的功率谱值为：
论
著
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；
3
平均周期图法
x ( n )e
n 0
N 1
j n
为了减小随机起伏，M.S.巴特利特提出平均周期图法， 即先把 N 点信号序列 x（n）分为 L 段，每段长度为 M。对 每段分别计算其周期图， 然后取各个周期图的平均作为功率
ˆ 1 谱的估值，可得功率谱值为： P ML
，进行谱估计：
2 ˆ ( ) 1 X (e j ) 1 P N N N
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 页）
1
引言
功率谱估计就是用已观测到的一组数据去估计一个平 稳信号的功率谱， 目的是根据有限数据在频域内提取被淹没 在噪声中的有用信号， 它对于认识一个随机信号或其他应用 方面来讲都是重要的，是数字信号处理的重要研究内容之 一， 在雷达系统追踪目标速度信息、 声呐系统寻找目标信息、 语音识别中去除噪声及进一步声学处理、生物学、天文学等 众多领域有着重要的应用。 功率谱估计可以分为经典谱估计 （非参数估计）和现代谱估计（参数估计） 。经典谱估计中 的周期图法是用得较多且最具代表性的方法。 本文首先研究 了基本的周期图法，接着针对它谱分辨率比较低等缺点，利 用 Matlab 实现了几种改进的周期图法。
ˆ P
1 MUL

i 1 n 0
L
M
2 i xm (n)e jn
，其中 U
(n) ，ω
n
（n）为窗函数。 为了得到较好的功率谱估值， 加窗和平均处理均应兼顾 减小随机起伏和保证有足够的谱分辨率两个方面。 图 3 给出 了 Welch 法周期图法得到的功率谱。
图4
多窗口周期图法得到的功率谱
8 而提升城市竞争力的重要产业支撑[9]。西部大开发以来的成 都实践表明，金融业、商务服务等现代服务业发展在解决就 业、推动产业结构优化升级方面发挥了重要作用。成都应深 入推进服务业结构战略性优化调整，加快发展现代物流、商 务服务、会展服务、文化创意、信息与科技服务等先导服务 业，不断提升商贸、金融、旅游等支柱服务业，大力拓展电 子商务、健康产业、环保服务等新兴服务，构建现代服务业 体系。立足“全域成都”资源条件和产业基础，优化现代服务 空间布局， 加快推进现代服务业市级战略性功能区与现代服 务业重点集聚区建设， 形成分工明确、 优势互补的布局导向， 构建“双核引领、三城辐射、两带带动”的战略性空间格局。 推进服务综合改革，创新服务业发展方式，完善服务业要素 保障体系，引导服务业向信息化、网络化发展，提升现代服 务在经济联系中的连通性和协同性功能。 3.3 发挥雄厚的产业积淀优势，做优做强特色优势产业 从“三线”建设到西部大开发，成都凭借优越的区位条件 和富庶的自然资源，积淀形成实力雄厚的产业基础，成为全 市经济社会发展的支撑。当前及未来很长一段时间内，传统 的特色优势产业仍将主导成都产业的发展。 这就需要进一步 明确电子信息产品制造、汽车、机械、石化、食品、冶金、 建材、家具以及制鞋等轻化工的特色优势产业地位，加大本 土企业培育力度，塑造本地企业竞争优势。积极支持生产要 素向高端技术、高端产品领域集聚，鼓励企业加大知识技术
通过实验可以看出，基本的周期图法谱分辨率较低，会 有很大的失真。平均周期图法和 Welch 法可以减小随机起 伏，收敛性较好，曲线平滑，估计的结果方差较小，但是功 率谱主瓣较宽， 分辨率低。 但是， 如果信号序列不是足够长， 由于每段序列长度变短， 功率谱估值对不同频率成分的分辨 能力也随之下降。 多窗口的周期图法功率谱比平均周期图法 和 Welch 法的功率谱谱峰窄一些，谱峰高度也有所增大，因 此多窗口的周期图法的功率谱表现得更加准确。 普通周期图 法的功率谱估计利用单一窗口， 所以在序列的始端和末端均 会丢失相关信息， 而多窗口的周期图法以增加窗口来充分利 用这些相关信息。
plot（k，plot_Pxx） xlabel（'Hz'） ； ylabel（'功率谱 dB'） ； ； title（'用基本周期图法（直接计算）计算功率谱'） grid on； 基本的周期图法可以提高计算效率， 不需要计算自相关 函数，但谱分辨率较低。
周期图的基本原理是对观测到的数据直接进行傅立叶 变换，然后取模的平方就是功率谱。取平稳随机信号 x（n） 的有限个观察点 x（0） 、x（1） 、…x（n-1） ，则傅立叶变换 为： X N (e j )
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） ； %观察数据点数