[转载]matlab的功率谱计算正确表达
已有 7044 次阅读
|个人分类:|系统分类:
功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。
但是周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。matlab中，welch法用函数psd实现。
调用格式如下：
[Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)X：输入样本数据NFFT：FFT点数Fs:采样率WINDOW：窗类型NOVERLAP，重叠长度现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多，这里不再详述，可以参考有关资料。matlab中，现代谱估计的很多方法都可以实现。music方法用pmusic命令实现；pburg函数利用burg法实现功率谱估计；pyulear函数利用yule-walker算法实现功率谱估计等等。另外，sptool工具箱也具有功率谱估计的功能。窗口化的操作界面很方便，而且有多种方法可以选择。t=linspace(0,1,1000);signal=4*sin(2*pi*50*t)+5*sin(2*pi*200*t);noise=randn(size(t));symbol=signal+Y=fft(symbol,128);f=)/128;P1=Y.*conj(Y)/128; %直接法subplot(1,3,1);plot(f,10*log10(P1(1:65)));xlabel('frequency')ylabel('power')title('直接法') 直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。MFs=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));
%间接法cxn=xcorr(symbol,'unbiased'); %计算序列的自相关函数P2=fft(cxn,128);subplot(1,3,2);plot(f,10*log10(P2(1:65)));xlabel('frequency')ylabel('power')title('间接法')
%加窗法window2=blackman(100); %blackman窗noverlap=20; %数据无重叠range='onesided'; %频率间隔为[0 1000/2]，只计算一半的频率[P3(1:65),f]=pwelch(symbol,window2,noverlap,128,1000,range);plot_P3=10*log10(P3(1:65));subplot(1,3,3);plot(f,plot_P3(1:65));xlabel('frequency')ylabel('power')title('加窗法')
间接法：间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。Matlab代码示例：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/plot_Pxx=10*log10(Pxx(index+1));plot(k,plot_Pxx);
改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。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/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_Pxxc plot_Pxx+plot_Pxxc]);
2. Welch法Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。Matlab代码示例：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);
figure(2)plot(f,plot_Pxx1);
figure(3)plot(f,plot_Pxx2);
本文引用地址：&此文来自科学网彭鹏博客，转载请注明出处。
上一篇：下一篇：
当前推荐数：0
评论 ( 个评论)
作者的其他最新博文
热门博文导读
Powered by
Copyright &AR模型功率谱估计的典型算法比较及MATLAB实现_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
AR模型功率谱估计的典型算法比较及MATLAB实现
阅读已结束，如果下载本文需要使用
想免费下载本文？
下载文档到电脑，查找使用更方便
还剩3页未读，继续阅读
你可能喜欢matlab里功率谱估计的幅度问题如x1=4*cos(2*pi*200*n)+4*cos(2*pi*320*n)+ 4*cos(2*pi*420*n); 三个不同频率同幅度的正选信号为什么作图 幅度不一样？_百度作业帮
matlab里功率谱估计的幅度问题如x1=4*cos(2*pi*200*n)+4*cos(2*pi*320*n)+ 4*cos(2*pi*420*n); 三个不同频率同幅度的正选信号为什么作图 幅度不一样？
matlab里功率谱估计的幅度问题如x1=4*cos(2*pi*200*n)+4*cos(2*pi*320*n)+&4*cos(2*pi*420*n);&三个不同频率同幅度的正选信号为什么作图&幅度不一样？
功率谱的话幅度应该是一样的，可能是你的数据处理有误。以下是一段Matlab验证程序：n=(1:;x1=4*cos(2*pi*200*n)+4*cos(2*pi*320*n)+ 4*cos(2*pi*420*n);plot(abs(fft(x1))*2/1000);你自己试一下吧。查看: 19249|回复: 73|关注: 0
Matlab功率谱估计的详尽分析——绝对原创
功率谱估计是信息学科中的研究热点,在过去的30多年里取得了飞速的发展。现代谱估计主要是针对经典谱估计(周期图和自相关法)的分辨率低和方差性能不好的问题而提出的。其内容极其丰富,涉及的学科和领域也相当广泛,按是否有参数大致可分为参数模型估计和非参数模型估计,前者有AR模型、MA模型、ARMA模型、PRONY指数模型等;后者有最小方差方法、多分量的MUSIC方法等。
ARMA谱估计叫做自回归移动平均谱估计,它是一种模型化方法。由于具有广泛的代表性和实用性,ARMA谱估计在近十几年是现代谱估计中最活跃和最重要的研究方向之一。
二： AR参数估计及其SVD—TLS算法。
谱分析方法要求ARMA模型的阶数和参数以及噪声的方差已知.然而这类要求在实际中是不可能提供的,即除了一组样本值x(1),x(2),…,x(T)以供利用(有时会有一定的先验知识)外,再没有其它可用的数据.因此必须估计有关的阶数和参数,以便获得谱密度的估计.在ARMA定阶和参数之估计中,近年来提出了一些新算法，如本文介绍的SVD—TLS算法便是其中之一。
三：实验结果分析和展望
1，样本数多少对估计误差的影响。（A=[1,0.8,-0.68,-0.46]）
& & 图1 上部分为N=1000；下部分为取相同数据的前N=50个数据产生的结果。
图1 N数不同：子图一N=1000,子图二N=200，子图三 N=50
由图可知，样本数在的多少，在对功率谱估计的效果上有巨大的作用，特别在功率谱密度函数变化剧烈的地方，必须有足够多的数据才能完整的还原原始功率谱密度函数。
2，阶数大小对估计误差的影响。
A=[1,-0.9,0.76]
A=[1,-0.9,0.76,-0.776]
图二 阶数为二阶和三阶功率密度函数图
A=[1,-0.9,0.86,-0.96,0.7]
A=[1,-0.9,0.86,-0.96,0.7,-0.74]
图三 阶数为三阶和四阶功率密度函数图
如图所示，阶数相差不是很大时，并不能对结果产生较大的影响。但是阶数太低，如图二中二阶反而不能很好的估计出原始值。
3，样本点分布对估计误差
对于相同的A=[1,-0.9,0.86,-0.96,0.7]；样本的不同，在估计时的误差是不可避免的。因此，我们在取得样本时，应该尽可能的减少不必要的误差。
图四：不同的样本得到不同的估计值
4，奇异值的阈值判定范围不同对结果的影响。
上图是取奇异值的阈值大于等于0.02，而下图是取阈值大于等于0.06，显然在同种数据下，阈值的选取和最终结果有密切关系。由于系数矩阵和其真实值的逼近的精确度取决于被置零的那些奇异值的平方和。所以选取太小，导致阶数增大，选取太大会淘汰掉真实的系数。根据经验值，一般取0.05左右为最佳。
14:33 上传
点击文件名下载附件
1.03 KB, 下载次数: 1057
好帖，必须支持LZ
支持原创，谢谢楼主分享，可是没有图~~~~(&_&)~~~~
楼主是实验报告吗，很想看一下，可以以附件形式发上来吗，或者发到我邮箱可以吗？
我的邮箱是
O(∩_∩)O谢谢
回复 1# logic115634 的帖子
好贴必须顶之
好东西啊，支持
正在学习，支持，下来学习学习
是要顶的啦，，那么好的贴，，可惜没有看见图，，，
收下了,慢慢学习,谢谢楼主
找的好辛苦啊，终于找到了！
谢谢LZ的分享！！！
好帖，顶上去，
可是，还得好好研究一下才行
站长推荐 /2
Powered by苹果/安卓/wp
苹果/安卓/wp
积分 66, 距离下一级还需 19 积分
道具: 彩虹炫, 雷达卡, 热点灯, 雷鸣之声, 涂鸦板, 金钱卡, 显身卡下一级可获得
权限: 自定义头衔
购买后可立即获得
权限: 隐身
道具: 金钱卡, 雷鸣之声, 彩虹炫, 雷达卡, 涂鸦板, 热点灯
利用一组数据估计其分布函数参数时出错，请高手指教估计的正确方法和结果，谢谢！
本帖是的后续工作。
原始数据：
<font color="#.008
<font color="#.008348
<font color="#.0092
<font color="#.0099
<font color="#.012
<font color="#.0123
<font color="#.0124
<font color="#.0145
<font color="#.0146
<font color="#.0148
<font color="#.0187
<font color="#.0207
<font color="#.0217
<font color="#.0228
<font color="#.023
<font color="#.0239
<font color="#.0244
<font color="#.0246
<font color="#.0253
<font color="#.0268
<font color="#.0273
<font color="#.0293
<font color="#.03
<font color="#.031
<font color="#.0319
<font color="#.0337
<font color="#.0349
<font color="#.0355
<font color="#.0358
<font color="#.037
<font color="#.0388
<font color="#.0388
<font color="#.0419
<font color="#.042
<font color="#.0422
<font color="#.044
<font color="#.0446
<font color="#.045
<font color="#.045
<font color="#.0492
<font color="#.0502
<font color="#.0506
<font color="#.0515
<font color="#.0516
<font color="#.0529
<font color="#.055
<font color="#.0564
<font color="#.0584
<font color="#.0591
<font color="#.0594
<font color="#.0595
<font color="#.0597
<font color="#.0637
<font color="#.0646
<font color="#.0656
<font color="#.0685
<font color="#.0686
<font color="#.0686
<font color="#.0737
<font color="#.0756
<font color="#.0759
<font color="#.0778
<font color="#.0798
<font color="#.0821
<font color="#.085
<font color="#.0859
<font color="#.0879
<font color="#.0885
<font color="#.0906
<font color="#.0919
<font color="#.0938
<font color="#.0984
<font color="#.099
<font color="#.0994
<font color="#.1008
<font color="#.1022
<font color="#.103
<font color="#.1062
<font color="#.1065
<font color="#.1126
<font color="#.1134
<font color="#.1175
<font color="#.1232
<font color="#.1346
<font color="#.1362
<font color="#.1373
<font color="#.1729
<font color="#.1736
<font color="#.1779
<font color="#.1782
<font color="#.2139
<font color="#.2168
<font color="#.2333
<font color="#.3126
<font color="#.3148
<font color="#.3178
<font color="#.3178
运行过程：
&& mixedpdf=@(x,mu1,mu2,s1,s2,rho)(rho*normpdf(x,mu1,s1)+(1-rho)*normpdf(x,mu2,s2))
fnegpdf=@(x,k0,k1,lamda,mu1,mu2,s1,s2,rho)((1-k0*exp(k1*x))* mixedpdf(x,mu1,mu2,s1,s2,rho))
fpospdf=@(x,k0,k1,lamda,mu1,mu2,s1,s2,rho)(lamda*exp(-lamda*x)*(rho*k0* normcdf((-mu1-k1*s1^2)/s1)*exp(mu1*k1+s1^2*k1^2/2)+(1-rho)*k0*normcdf((-mu2-k1*s2^2)/s2)*exp(mu2*k1+s2^2*k1^2/2))+ mixedpdf(x,mu1,mu2,s1,s2,rho))
roepdf=@(x)fnegpdf.*(x&0)+fpospdf.*(x&=0)
mixedpdf =
@(x,mu1,mu2,s1,s2,rho)(rho*normpdf(x,mu1,s1)+(1-rho)*normpdf(x,mu2,s2))
@(x,k0,k1,lamda,mu1,mu2,s1,s2,rho)((1-k0*exp(k1*x))* mixedpdf(x,mu1,mu2,s1,s2,rho))
@(x,k0,k1,lamda,mu1,mu2,s1,s2,rho)(lamda*exp(-lamda*x)*(rho*k0* normcdf((-mu1-k1*s1^2)/s1)*exp(mu1*k1+s1^2*k1^2/2)+(1-rho)*k0*normcdf((-mu2-k1*s2^2)/s2)*exp(mu2*k1+s2^2*k1^2/2))+ mixedpdf(x,mu1,mu2,s1,s2,rho))
@(x)fnegpdf.*(x&0)+fpospdf.*(x&=0)
&& [phat1,pci1]=mle(gap_roe,'pdf',roepdf,'start',[.9,25,120,-0.0,0.9,0.3012])
??? Error using ==& stats\private\mlecustom&checkFunErrs
The following error occurred while trying to evaluate
the user-supplied pdf function '@(x)fnegpdf.*(x&0)+fpospdf.*(x&=0)':
Error using ==& @(x)fnegpdf.*(x&0)+fpospdf.*(x&=0)
Too many input arguments.
Error in ==& stats\private\mlecustom at 159
checkFunErrs('pdf',pdfFun,start,uncensData,[],[],pdfAddArgs);
Error in ==& mle at 219
[phat, pci] = mlecustom(data,varargin{:});
附件1为主要的推导过程（其中α在实际估计时取为0），最终要估计（5）式中的参数。
附件2为参考文献：王亚平，吴联生，白云霞，2005，中国上市公司盈余管理的频率与幅度，经济研究，2005年 12期
所有公式推导均引自参考文献第三部分。
(34.54 KB)
17:09:14 上传
17:09:14 上传
载入中......
自己先顶一下，差点看不见了
悄悄地用力再顶一下
再用一下力
本帖最后由 yuyi3860 于
11:49 编辑
错误的关键点可能在于对概率密度函数的定义，这里定义了一个分段概率密度函数
roepdf=@(x)fnegpdf.*(x&0)+fpospdf.*(x&=0)
问题在于，在matlab中应如何定义是分段函数的概率密度函数？
论坛好贴推荐
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
为做大做强论坛，本站接受风险投资商咨询，请联系（010-）
邮箱：service@pinggu.org
合作咨询电话：(010)
广告合作电话：（刘老师）
投诉电话：(010)
不良信息处理电话：(010)
京ICP证090565号
京公网安备号
论坛法律顾问：王进律师