当前位置：文档库 › (完整版)功率谱估计性能分析及Matlab仿真

(完整版)功率谱估计性能分析及Matlab仿真

功率谱估计性能分析及Matlab 仿真

1 引言

随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。

信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的

N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。

2 经典功率谱估计

经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。

2.1 周期图法( Periodogram )

Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。

取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换

()()N j j n N n X e

x n e ω

ω---==∑

然后进行谱估计

1()()j N S X e N

ωω-=

周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT 快速算法来计算。但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据。同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低[1]。

该方法基于Matlab 实现的程序： clear all; load test x; N=4096;

Fn=-0.5:1/N:0.5-1/N; px=fft(x,N);

pmax=max(px);%归一化 px=px/pmax;

px=10*log10(px+0.000001); plot(Fn,fftshift(px));grid on;

图1 周期图法 4096N =

图2 周期图法 128N =

说明：

(1) 本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献[3]的test.dat 。该数据为128点复序列（图3），由复数噪声加上四个复正弦组成。其归一化频率分别是：12340.15,0.16,0.252,0.16f f f f ====-。

图3 复序列 (),[0,127]x n n ∈

(2) 从仿真图可以清晰看到，1f 和2f 不能完全分开，仅在波形的顶部能看出

是两个频率分量；此外，当数据长度N 太大时（图1），谱曲线呈现较大的起伏；当数据长度N 太小时（图2），谱的分辨率又不好。据此，周期图法不满足一致性估计条件。

2.2 自相关法( BT 法)

自相关法的理论基础是维纳—辛钦定理。1958年Blackman 和Tukey 给出了这一方法的具体实现。

对于平稳随机信号来说，其自相关函数是确定性函数，故其功率谱也是确定的。这样可由平稳随机离散信号的有限个离散值(0),(1),...,(1)x x x n -求出自相关函数

()()()N m x n R m x n x n m N

--==

+∑

然后在(,)M M -内对()x R m 做傅里叶变换，得到功率谱

()()M

j n x m M

S R m e ωω-=-=

∑

该方法基于Matlab 实现的程序： clear all; load test x; N=4096;

Fn=-0.5:1/N:0.5-1/N; Mlag=64;

rx=xcorr(x,Mlag,'unbiased'); px=fft(rx,N);

pmax=max(px);%归一化 px=px/pmax;

px=10*log10(px+0.000001); plot(Fn,fftshift(px)); grid on;

图4 自相关法不加窗64

图5 自相关法不加窗32

图6 自相关法使用汉明窗( Hamming )

说明：

(1) 该方法先由序列()x n 估计出自相关函数()x R m ，然后对()x R m 进行傅里叶变换，便得到()x n 的功率谱估计。当延迟与数据长度之比很小时，可以有良好的估计精度。

(2) 图4是用自相关法（BT 法）求出的功率谱，64M =没有加窗；图5也是用自相关法（BT 法）求出的功率谱，32M =，没有加窗；图6同样是采用自相关法求出的功率谱，32M =，使用了汉明窗。显然，自相关函数的延迟M 越小，谱变得越平滑。

2.3 Welch 法

该方法的基本原理是在对随机序列分段时，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均。这样可得功率谱

110

1()()L M i

j m i n S x n e MUL ωω--===∑∑ 其中

0()M n U n ω-==∑

为窗函数。

这里()n

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

Fn=-0.5:1/N:0.5-1/N;

xpsd=pwelch(x,hamming(33),16,N,'whole');

mmax=max(xpsd);%归一化

xpsd=xpsd/mmax;

xpsd=10*log10(xpsd+0.000001);

plot(Fn,fftshift(xpsd));

grid on;

图7 Welch法不叠合使用汉明窗( Hamming )

图8 Welch法叠合16点使用汉明窗( Hamming )

图9 Welch法叠合16点使用矩形窗( Boxcar )

图10 Welch法叠合16点使用布莱克曼窗( Blackman ) 说明：

(1) 因为Welch法允许各段数据交叠，所以数据段数L会增加，使方差得到更大的改善，但是数据的交叠减小了每一段数据的不相关性，使方差的减小不会达到理论程度。另外，采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率。

(2) 图7是利用Welch法求出的周期图，共分四段，每段32点，没有叠合，使用了汉明窗；图8也是利用Welch法求出的周期图，共分四段，每段32点，，使用了汉明窗；图9是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，且使用了矩形窗；图10是利用Welch法求出的周期图，共分四段，每段32点，交叠数为16，使用了布莱克曼窗。从图中可以看出，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好，但是其起伏性较大，所以其方差特性最差。由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。因此，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同。2.4 经典功率谱估计的性能比较

由以上的Matlab仿真图形和相关结果分析，我们得到了经典谱估计算法性能的直观比较：

(1) 周期图法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰。

(2) 自相关法（BT法）由于使用了平滑窗对周期图法估计的功率谱进行了平滑，因此方差性能较好，功率谱比周期图法估计的要平滑，但其分辨率比周期图法低。

(3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。

综合上述讨论，我们可以对经典谱估计的算法作大致的总结[3]：

(1) 功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而仍是目前较常用的谱估计方法。

π，N是所使用的数据长度。

(2) 谱的分辨率较低，它正比于2/N

(3) 方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。

(4) 周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率且增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。

3 现代功率谱估计

由前一章的讨论我们可知，经典功率谱估计方法的方差性能较差，分辨率较低。而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。参数模型法是现代谱估计的主要内容，参数模型主要分为AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能，因此是被研究最多并获得广泛应用的一种模型。本报告中现代功率谱估计的仿真基于的是AR模型。

3.1 自相关法

假定观察到得数据为(0),(1),...,(1)

x x x n-，而对于无法观察到得区间（即和），()

<>-

n n N

x n的样本假定为0，观测数据区间之外的数据为0，在均方

R是Toeplitz矩阵，误差意义下使得数据的预测误差功率最小。由于自相关矩阵?

而且又为正定的，故可利用Levinson-Durbin递归算法高效求解，得到AR模型

参数。

该方法基于Matlab实现的程序：

clear all;

load test x;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pyulear(x,20,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));

grid on;

图11 自相关法10

图12 自相关法 20p =

图13 自相关法 30p =

说明：

(1) 图11、12和13是用自相关法求出的AR 谱曲线，阶次p 分别等于10，20和30。可以看出，在阶次较低时（图11），分辨率和检测能力均不好。当30p =时，1f 和2f 处的两个正弦刚刚可以分开，在3f 和4f 处的两个正弦也可以检出。

因此必须通过提高阶次来达到分辨出间隔较小的频率点的效果。

(2) AR 模型的自相关法等效于对前向预测的误差序列前后加窗，加窗的结果是使得自相关法的分辨率降低。数据越短，分辨率越不好。

3.2 协方差法

协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取得不同。由于协方差法对于观察区间[]0,1N -外()x n 样本并未假定为0，故预测误差功率表达式中的()x n k -总是落在观察区间[]0,1N -中，为此预测误差功率的求

和上下限必须在[],1p N -之间。但由此得到的自相关矩阵?p R 是对此的半正定矩阵，且不具有Toeplitz 性质，故不能采用Levinson-Durbin 递归算法求解，因此得到的AR 模型可能不稳定。

该方法基于Matlab 实现的程序： clear all; load test x;

N=4096;tn=-0.5:1/N:0.5-1/N; xpsd=pcov(x,10,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); plot(tn,fftshift(xpsd)); grid on;

图14 协方差法 10p =

说明：

(1) 可以看到谱图在信号源频率处：12340.15,0.16,0.252,0.16f f f f ====-谱线狭窄突出，其他处谱线起伏较为剧烈。

(2) 采用协方差法对信号进行建模，能够较好地反映出信号真正的模型。

3.3修正的协方差法

AR 谱估计的协方差算法基于的是最小化前向预测误差。而修正的协方差算法基于的是最小化前向和后向预测误差。这样使得它的误差功率的计算是在相对于协方差法多一倍的数据点上进行，这在观察数据长度很短的情况下，是非常有利的，但这要求信号在正反两个方向上呈现相同的特性。此外，由此得到的自相

关矩阵?p R 不具有Toeplitz 性质，故其正则方程不能采用Levinson-Durbin 递归算法求解。

该方法基于Matlab 实现的程序： clear all; load test x; N=4096;

fn=-0.5:1/N:0.5-1/N; xpsd=pmcov(x,10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));grid on;

图15 修正的协方差法10

图16 修正的协方差法15

说明：

(1) 修正的协方差法较协方差法而言，谱估计图大致相同，但前者在信号源频率处的谱峰更加突出、尖锐，易于辨别。

(2) 通过图15和16的对比，我们可发现，在阶次较高的情况下能得到非常满意的结果。

3.4 Burg法

Burg算法是较早提出的建立在数据基础上的AR系数求解的有效算法。它基于最小化前向后向预测误差的同时满足Levinson-Durbin递归。对比其它的AR 估计方法，Burg法避免了对自相关函数的计算，改而直接估计反射系数。

对于短数据的估计，Burg法求出的AR功率谱密度估计非常逼近于真值。另外，它能确保产生一个稳定的AR模型，并且能高效计算。

Burg法由于具有上述优点，所以分辨率比自相关法高，但对于混有白噪声的正弦信号，有时可能会出现谱线分裂现象。

该方法基于Matlab实现的程序：

clear all;

N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pburg(x,10,N);

pmax=max(xpsd);

xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001);

plot(fn,fftshift(xpsd));

grid on;

图17 Burg法10

图18 Burg法15

说明：

(1) 图17和18是采用Burg算法对数据文件所做的功率谱估计，在阶次p=的情况下，得到了非常满意的效果。

(2) 通过图13和17的对比，我们可明显看出Burg法的分辨率比自相关法高，在阶次较低的情况下（图17）也能较好的分辨出间隔小的频率点。

3.5 MUSIC 法

刚才所讨论的自相关法、协方差法和Burg 法都是基于参数建模的功率谱估计，而基于非参数建模的功率谱估计也是现代功率谱估计的重要内容。该方法是基于自相关矩阵的特征分析或者特征值分解的功率谱估计，它将相关矩阵的特征向量空间分解为信号子空间和噪声子空间，由此衍生出EV( Eigenvector )算法与MUSIC( Multiple Signal Classification )算法的信号功率谱估计[4]。其中EV 谱估计与MUSIC 算法谱估计都是基于噪声子空间的功率谱估计。这类方法对线谱（正弦信号的谱）最合适，对检测混有白噪声的正弦信号很有效，特别是低信噪比的情况。

MUSIC 估计由下面方程给出

()()()

()

MUSIC p p H H H k k k k M k M P V e e V V e ωωωω++=+=+=

???

∑

此处()e ω是复正弦信号向量。该方法基于Matlab 实现的程序： clear all; load test x; N=4096;

fn=-0.5:1/N:0.5-1/N; xpsd=pmusic(x',10,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); for i=1:N

xxpsd(i)=xpsd(N+1-i); end

plot(fn,fftshift(xxpsd)); grid on;

图19 MUSIC法10

说明：

(1) 从谱图上看，MUSIC算法的谱图较为圆滑，在信号源频率处，谱线峰值明显，但峰带频率跨度较大。

(2) 通过13和19的对比，我们可明显看出MUSIC法的分辨率要好于AR 模型的自相关法，在阶次较低的情况下也能较好的分辨出间隔小的频率点。

4 总结

本报告首先讨论了经典功率谱估计中三种主要方法的定义、算法及估计性能，其次讨论了现代功率谱估计中的主要内容——参数模型法谱估计，最后简要介绍了基于相关阵特征分解的谱估计算法。对所有的谱估计算法，报告均给出了Matlab程序、仿真结果及相应的分析。

参考文献

[1] 王春兴．基于MA TLAB实现经典功率谱估计[J] ．曲阜师范大学学报，2011,37(2) ．

[2] 宋宁，关华．经典功率谱估计及其仿真[J]．现代电子技术，2008,31(11)．

[3] 胡广书．数字信号处理[M]．北京：清华大学出版社，2003．

[4] 王福杰，潘宏霞．MA TLAB中几种功率谱估计函数的分析比较与选择[J]．计算机科学

与技术，2009,27(6)．

matlab功率谱估计

功率谱估计及其MATLAB仿真 1经典功率谱估计经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1(Btfangfa.M)： Fs=500;%采样频率 n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));%产生含有噪声的序列 nfft=512; cxn=xcorr(xn,'unbiased');%计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); %Round towards nearest integer. k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); figure(1); plot(k,plot_Pxx); 结果如下： 1.2周期图法(periodogram) 周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例2(PEriod.M)： Fs=600; n=0:1/Fs:1;

matlab实现功率谱密度分析psd

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

(完整版)功率谱估计性能分析及Matlab仿真

功率谱估计性能分析及Matlab 仿真 1 引言随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。 2 经典功率谱估计经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。周期图法( Periodogram ) Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换 1 ()()N j j n N n X e x n e ω ω---==∑ 然后进行谱估计

功率谱估计介绍(介绍了matlab函数)

功率谱估计介绍谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。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工具箱也具有功率谱估计的功能。窗口化的操作界面很方便，而且有多种方法可以选择在海杂波抑制的研究中，对海杂波谱分析一定要用到谱估计理论，一定得花时间学好！

功率谱密度估计方法的MATLAB实现

功率谱密度估计方法的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('幅度');

AR功率谱估计MatlAB

AR模型的谱估计是现代谱估计的主要内容 AR模型的谱估计是现代谱估计的主要内容。 1．AR 模型的Yule—Walker方程和Levinson-Durbin递推算法：在MATLAB中，函数levinson和aryule都采用 Levinson-Durbin递推算法来求解AR模型的参数a1,a2,……,ap及白噪声序列的方差，只是两者的输入参数不同，它们的格式为： A=LEVINSON(R,ORDER) A=ARYULE(x,ORDER) 两函数均为定阶ORDER的求解，但是函数levinson的输入参数要求是序列的自相关函数，而函数aryule的输入参数为采样序列。下面语句说明函数levinson和函数aryule的功能是相同的：例子： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; r=xcorr(x,'biased'); r(1:length(x)-1)=[]; A=levinson(r,5) B=aryule(x,5) 2．Burg算法：格式为：A=ARBURG(x,ORDER); 其中x为有限长序列，参数ORDER用于指定AR 模型的阶数。以上面的例子为例： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=arburg(x,5)

3.改进的协方差法：格式为：A=ARMCOV(x,ORDER); 该函数用来计算有限长序列x(n)的ORDER阶AR 模型的参数。例如：输入下面语句： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=armcov(x,5) AR模型阶数P的选择： AR 模型阶数P一般事先是不知道的，需要事先选定一个较大的值，在递推的过程中确定。在使用Levinson—Durbin递推方法时，可以给出由低阶到高阶的每一组参数，且模型的最小预测误差功率Pmin（相当于白噪声序列的方差）是递减的。直观上讲，当预测误差功率P达到指定的希望值时，或是不再发生变化时，这时的阶数即是应选的正确阶数。因为预测误差功率P是单调下降的，因此，该值降到多少才合适，往往不好选择。比较常见的准则是：最终预测误差准则：FPE(r)=Pr{[N+(r+1)]/ [N-(r+1)]} 信息论准则：AIC(r)=N*log(Pr)+2*r 上面的N为有限长序列x(n)的长度，当阶数r由1增加时，FPE(r) 和AIC(r)都将在某一r处取得极小值。将此时的r定为最合适的阶数p。 MATLAB中AR模型的谱估计的函数说明： 1． Pyulear函数：功能：利用Yule--Walker方法进行功率谱估计. 格式： Pxx=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT,Fs) Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)

功率谱估计的MATLAB实现

实验功率谱估计实验目的: 1、掌握最大熵谱估计的基本原理。 2、了解最终预测误差（FPE）准则。 3、掌握周期图谱估计的基本原理。 4、掌握传统谱估计中直接法与间接法之间的关系。 5、复习快速傅里叶变换与离散傅里叶变换之间关系。实验内容: 1、设两正弦信号的归一化频率分别为0.175和0.20，用最大熵法编程计算信噪比S/N=30dB、N=32点时该信号的最大熵谱估计结果。 2、用周期图法编程计算上述信号的谱估计结果。程序示例: 1、最大熵谱估计 clc; N=32; SNR=30; fs=1; t=1:N; t=t/fs; y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t); x = awgn(y,SNR); M=1; P(M)=0; Rx(M)=0; for n=1:N P(M)=P(M)+(abs(x(n)))^2; ef(1,n)=x(n); eb(1,n)=x(n); end P(M)=P(M)/N; Rx(M)=P(M); M=2;

A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); TH=FPE(M-1); for n=M:N ef(M,n)=ef(M-1,n)+xishu*eb(M-1,n-1); eb(M,n)=eb(M-1,n-1)+xishu*ef(M-1,n); end M=M+1; A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); for m=1:M-2 a(M-1,m)=a(M-2,m)+xishu*a(M-2,M-1-m); end while FPE(M-1)

matlab求功率谱

matlab实现经典功率谱估计 fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数 matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch法估计，即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。 1、直接法：直接法又称周期图法，它是把随机序列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)); 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;

(完整word版)自己编写算法的功率谱密度的三种matlab实现方法

功率谱密度的三种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实现经典功率谱估计]matlab功率谱估计

[matlab实现经典功率谱估计]matlab功率谱估计 1、直接法：直接法又称周期图法，它是把随机序列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)); 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); 3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。3.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; %数据无重叠

MATLAB处理信号得到频谱、相谱、功率谱全解

第一：频谱一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意：（1）函数FFT返回值的数据结构具有对称性。例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) → Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。（2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。二.FFT应用举例例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列

subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果： fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：（1）数据个数N=32，FFT所用的采样点数NFFT=32；（2）N=32，NFFT=128；（3）N=136，NFFT=128；（4）N=136，NFFT=512。

功率谱估计MATLAB实现

功率谱估计性能分析及其MATLAB实现一、经典功率谱估计分类简介 1．间接法根据维纳-辛钦定理，1958年Blackman和Turkey给出了这一方法的具体实现，即先由N个观察值，估计出自相关函数，求自相关函数傅里叶变换，以此变换结果作为对功率谱的估计。 2．直接法直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N 个观察值直接进行傅里叶变换，得到，然后取其幅值的平方，再除以N，作为对功率谱的估计。 3．改进的周期图法将N点的观察值分成L个数据段，每段的数据为M，然后计算L个数据段的周期图的平均，作为功率谱的估计，以此来改善用N点观察数据直接计算的周期图的方差特性。根据分段方法的不同，又可以分为Welch法和Bartlett法。 Welch法所分的数据段可以互相重叠，选用的数据窗可以是任意窗。 Bartlett法所分的数据段互不重叠，选用的数据窗是矩形窗。

二、经典功率谱估计的性能比较 1．仿真结果为了比较经典功率谱估计的性能，本文采用的信号是高斯白噪声加两个正弦信号，采样率Fs=1000Hz，两个正弦信号的频率分别为f1=200Hz，f2=210Hz。所用数据长度N=400. 仿真结果如下： Figure1(a)示出了待估计信号的时域波形；

Figure2(b)示出了用该数据段直接求出的周期图，所用的数据窗为矩形窗； Figure2(c)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为矩形窗，长度M=128，数据没有加窗； Figure2(d)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为Hamming 窗，长度M=64，数据没有加窗； Figure2(e)是用Welch平均法求出的功率谱曲线，每段数据的长度为64点，重叠32点，使用的Hamming窗； Figure2(f)是用Welch平均法求出的功率谱曲线，每段数据的长度为100点，重叠48点，使用的Hamming窗； 2．性能比较 1)直接法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰； 2)间接法由于使用了平滑窗对直接法估计的功率谱进行了平滑，因此方差性能比直接法好，功率谱比直接法估计的要平滑，但其分辨率比直接法低。 3)Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。 3．关于经典功率谱估计的总结 1)功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而仍是目前较常用的谱估计方法。 2)谱的分辨率较低，它正比于2π/N，N是所使用的数据长度。 3)方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。 4)周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图的方差性能，但往往又减小了分辨率和增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。三、AR模型功率谱估计 1．A R模型功率谱估计简介 AR模型功率谱估计是现代谱估计中最常用的一种方法，这是因为AR模型参数的精确估计可以用解一组线性方程（Yule-Walker方程）的方法求得。其核心思想是：将信号看成是一个p 阶AR过程，通过建立Yule-Walker方程求解AR模型的参数，从而得到功率谱的估计。由于已知的仅仅是长度有限的观测数据，因此AR模型参数的求得，通常是首先通过某种算法求得自相关函数的估计值，进而求得AR模型参数的估计值。常用的几种AR模型参数提取方法有： 1)自相关法假定观测数据区间之外的数据为0，在均方误差意义下使得数据的前向预测误差最小。

基于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太小，谱的分辨率又不好，因此需要改进。 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));

MATLAB中AR模型功率谱估计中AR阶次估计的实现

MATLAB中AR模型功率谱估计中AR阶次估计的实现（最近看了几个关于功率谱的问题，有关AR模型的谱估计，在此分享一下，希望大家不吝指正）（声明：本文内容摘自我的毕业论文——心率变异信号的预处理及功率谱估计）（按：AR模型功率谱估计是对非平稳随机信号功率谱估计的常用方法，但是其模型阶次的估计，除了HOSA工具箱里的arorder函数外，没有现成的函数可用，arorder函数是基于矩阵SVD分解的阶次估计方法，为了比较各种阶次估计方法的区别，下面的函数使用了'FPE', 'AIC', 'MDL', 'CAT'集中准则一并估计，并采用试验方法确定那一个阶次更好。） ………………………………以上省略…………………………………………………………………… 假设原始数据序列为x，那么n阶参数使用最小二乘估计在MATLAB中实现如下：复制内容到剪贴板代码: Y = x; Y(1:n) = []; m = N-n; X = [];% 构造系数矩阵 for i = 1:m for j = 1:n X(i,j) = xt(n+i-j); end end beta = inv(X'*X)*X'*Y'; beta即为用最小二乘法估计出的模型参数。此外，还有估计AR模型参数的Yule-Walker方程法、基于线性预测理论的Burg算法和修正的协方差算法等[26]。相应的参数估计方法在MATLAB中都有现成的函数，比如aryule、arburg以及arcov等。 4.3.3 AR模型阶次的选择及实验设计

文献[26]中介绍了五种不同的AR模型定阶准则，分别为矩阵奇异值分解（Singular Value Decomposition, SVD）定阶法、最小预测定误差阶准则（Final Prediction Error Criterion, FPE）、AIC定阶准则（Akaika’s Information theoretic Criterion, AIC）、MDL定阶准则以及CAT定阶准则。文献[28]中还介绍了一种BIC定阶准则。SVD方法是对Yule-Walker方程中的自相关矩阵进行SVD分解来实现的，在MATLAB工具箱中arorder函数就是使用的该算法。其他五种算法的基本思想都是建立目标函数，阶次估计的标准是使目标函数最小化。以上定阶准则在MATLAB中也可以方便的实现，下面是本文实现FPE、AIC、MDL、CAT定阶准则的程序（部分）：复制内容到剪贴板代码: for m = 1:N-1 …… % 判断是否达到所选定阶准则的要求 if strcmp(criterion,'FPE') objectfun(m+1) = (N+(m+1))/(N-(m+1))*E(m+1); elseif strcmp(criterion,'AIC') objectfun(m+1) = N*log(E(m+1))+2*(m+1); elseif strcmp(criterion,'MDL') objectfun(m+1) = N*log(E(m+1))+(m+1)*log(N); elseif strcmp(criterion,'CAT') for index = 1:m+1 temp = temp+(N-index)/(N*E(index)); end objectfun(m+1) = 1/N*temp-(N-(m+1))/(N*E(m+1)); end if objectfun(m+1) >= objectfun(m) orderpredict = m; break; end end orderpredict变量即为使用相应准则预测的AR模型阶次。

用matlab做经典功率谱估计

用matlab做经典功率谱估计经典功率谱估计 1、直接法：直接法又称周期图法，它是把随机序列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)); 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);

3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N 太小，谱的分辨率又不好，因此需要改进。 3.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/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot_Pxxc=10*log10(Pxxc(index+1)); figure(1) plot(k,plot_Pxx); pause; figure(2) plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]); 3.2、Welch法 Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数 w(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));

Matlab功率谱计算

【转】matlab的功率谱计算功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合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模型等；非参数模型谱估计有最小方差法和

MATLAB仿真实现功率谱估计(采用自相关法)

数字信号处理课程实验报告实验指导教师：黄启宏实验名称 MATLAB 仿真实现功率谱估计（采用自相关法）专业、班级电子与通信工程姓名张帅实验地点仿古楼301 实验日期 2013.11.17 一、实验内容采用自相关法（参数模型）实现功率谱估计。二、实验目的（1）掌握自相关法（参数模型）估计出功率谱的步骤和方法；（2）在实验的过程中找到影响功率谱估计的因素；（3）比较自相关法（参数模型）和周期图法（直接法）估计功率谱的优点。三、实验原理（1）假定研究的()x n 过程是由一个输入序列()u n 激励一个线性系统()H z 的输出；（2）由已知的 ()x n ，或者其自相关函数()x r m 来估计()H z 的参数；（3）由()H z 的参数来估计()x n 的功率谱。四、涉及实验的相关情况介绍（包含使用软件或实验设备等情况）一台安装MATLAB 软件的电脑

五、实验记录程序、相关的图形、相关数据记录及分析）（ %自相关法（参数模型）实现功率谱估计； fs = 1000; n = 0:1/fs:0.3; xn = cos(200*pi*n)+0.1*randn(size(n)); nfft = 1024; window=boxcar(length(xn)); %输出随机信号xn； subplot(211); plot(n,xn);xlabel('时间');ylabel('幅度');title('输入信号x（n）'); axis([0 0.3 -2 2]); grid on; %自相关法功率谱估计 p1=floor(length(xn)/3)+1;%阶数 [px1,f1] = pyulear(xn,100,nfft,fs); px1max = max(px1); px1 = px1/px1max; px1=10*log10(px1+0.000001); [Pxx,f] = periodogram(xn,window,nfft,fs);%直接法功率谱估计PxxMAX =max(Pxx); Pxx = Pxx/PxxMAX;