关于求最大熵功率谱

弯弓射大雕

求出预测误差滤波系数rm,k和预测误差功率后，再怎样求最大熵功率谱，我看了计算公式为见附件；
我参考《振动数字信号处理程序库》第281页，这里提供的程序是通过直接求和（功率谱计算公式分母部分）来计算功率谱
我又看了胡广书提供的利用burg法求最大熵功率谱，他是给一个2^n谱线数，即在频率轴上显示的频率个数，然后用FFT法求分母部分
不知那个对呀，并且不知小弟理解的对不对，请教版主和各位高手
附件里是计算公式

zhlong

分母的那个求和项就是离散傅立叶变换的形式啊。
另外最大熵谱方法和AR法是等价的。
matlab里面，
Pxx=pyulear(x,order);
Pxx=pburg(x,order);
[Pxx,f]=pcov(x,20,[],Fs);
这些函数得到的Pxx都是你公式里的s(f)，你看一下这些函数的程序就可以知道怎么编写上述计算公式了。

弯弓射大雕

不知楼上说的看看这些函数的程序是啥意思呀，从matlab里是看不到的吧

[ 本帖最后由 zhlong 于 2007-10-22 13:48 编辑 ]

jingrenzhi

你好,《振动数字信号处理程序库》这本书是谁编的?哪里出版的?谢谢

zhlong

怎么看不到。你试过吗？

弯弓射大雕

中国电子仪器仪表学会信号处理学会编的，科学出版社1988年编的，最大熵谱的程序在281页，按照那里的程序，最大熵谱的分辨率是1/2/N/dt，N为采样点个数

弯弓射大雕

对，我用命令edit+pburg，看到一些程序，但没有看到求出预测误差滤波器系数和预测误差功率后，谱是怎么求得

[ 本帖最后由 弯弓射大雕 于 2007-10-22 15:08 编辑 ]

jingrenzhi

谢谢!这个资料有点老了!不妨查点新的!

zhlong

edit pburg
可以看到[Pxx,freq,msg,units,Sxx,options] = arspectra(method,x,p,varargin{:});
arspectra在pburg文件所在的文件夹里的一个名为private的文件夹内，在arspectra里面你又可以看到：
[Pxx,w,units] = computepsd(Sxx,w,options.range,options.nfft,options.Fs,'psd');
再去找computepsd，看这个函数名应该就是计算谱的了。
就是麻烦了一点，学东西嘛！

弯弓射大雕

感谢zhlong的热心解答，:victory:

songzy41

我在这里给出一个函数，便是按楼主笫一层给的公式来计算功率谱的：
function pf=lpcar2pf(ar,np)
%LPCAR2PF Convert AR coefs to power spectrum PF=(AR,NP)
% the power spectrum will contain np+2 points ranging from Dc to nyquist frequency
% The routine is faster if NP+1 is a power of 2.
[nf,p1]=size(ar);
if nargin<2 np=p1-1; end
pf=abs(rfft(ar.',2*np+2).').^(-2);

其中用到rfft(实数序列的FFT)，在附件中。RFFT

弯弓射大雕

谢谢宋老师，还有一个问题，就是Matlab中的改进协方差方法求最大熵谱的方法即[pxx,w]=Pmcov(x,order,nfft,fs),采用的算法是不是就是Marple1980年提出的那种方法呀。因为要往论文里写，所以要知道它的算法
谢谢了！

zhangnan3509

你可以好好看看这个函数的程序啊
比较一下算法就可以得出结论

sarain

function pf=lpcar2pf(ar,np)这个里面的ar是什么啊？是原序列还是什么系数？？

cjy1987

sarain 发表于 2009-6-8 17:29

                               
登录/注册后可看大图

function pf=lpcar2pf(ar,np)这个里面的ar是什么啊？是原序列还是什么系数？？

是Ar模型的系数