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

一、功率谱估计可以分为经典谱估计(非参数估计) 和现代谱估计(参数估计)。经典谱估计的方法主要方法有自相关估计法和周期图法以及对周期图的改进方法现代谱估计的内闭返容极其丰富，涉及的学科及应用领域也相当广泛，方法大致可分为参数模型谱估计和非参数模型谱估计，前者有AR 模型法(最大熵谱分析法)、MA模型，ARMA模型、Prony 指数模型等；后者有最小方差法，多分量的MUSIC方法等。其中周期图法和AR 模型法是用得较多且最具代表性的方法。从信号的来源分，又可分为一维谱估计、二维谱估计及多维谱估计。从使用的统计量来分，目前大部分工作是建立在二阶矩基础上的，但由于功率谱密度是频率的实函数，缺少相位信息，因此，建立在兄正高阶矩基础上的谱估计方法正引起人们的注意。从信号的特征来分，在这之前所说的方法都是对平稳随机信号而言，其谱分量不随时间变化，对非平稳随机信号，其谱是时变的，近20年来，以wigner分析为代表的时域分析引起了人们的广泛兴趣，形成了现代谱估计的一个新的研究领域。

二、MATLAB的简单介绍：

MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数羡态悔值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

好眼熟啊。。。通院备孙08级的随机信号实验题目？

这份是我在网备滚卖络上找的：

（1）.周期图法： 思想：周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度 N。由于序列 x(n)的离散傅里叶变换 X(k) 具有周期性，因而这种功率谱也具有周期性，常称为周期图。周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。

程序：

%首先，生成输入信号的程序为：

clear fs=20000 n=0:1/ fs: 0. 1N=lengt h( n)W=2000*pi%因方波频率 F=1000HZ 所以角频率 W=2000pi X1n=square( W*n) %方波信号 X2n=randn( 1, N)%白噪声信号 xn=X1n+0. 2*X2nsubplot (3, 1, 1) plot (n, xn) xlabel( ' n' ) ylabel('输入信仿逗号') %其次，开始用周期图法进行估计； clear all fs=20000n=0:1/fs:0.1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nsubplot(2,1,1) plot(n,xn)Nfft=256N=256%傅里叶变换的采样点数256 Pxx=abs(fft(xn,Nfft).^2)/Nf=(0:length(Pxx)-1)*fs/length(Pxx)subplot(2,1,2), plot(f,10*log10(Pxx)),%转成DB 单位

（2）.自相关函数法： 思想： 随机信号 x(n)的相关函数是在时间域内描述随机过程的重要特征。自相关函数是随机信号在不同时刻的值之间的依赖性的量度，是一个很有用的统计平均量，其定义为自相关函数 (1) 式中E【·】表示数学期望，*表示共轭值,m 为时间滞后数。在随机信号处理中，自相关函数可以用来检测淹没在随机噪声干扰中的信号，随机信号的自功率谱等于它的自相关函数的傅里叶变换。因此，通过自相关估计可求得信号的功率谱。利用计算机计算自相关估值有两种方法。一种是直接方法，先计算出随机信号和它的滞后序列的乘积，再取其平均值即得相关函数的估计值。另一种是间接方法，先用快速变换算法计算随机序列的功率谱密度，再作反变换计算出相关函数。

程序： n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(2) subplot(3,1,1) plot(n,xn)%输入信号 axis([0 0.01 -5 5]) m =-100:100 [r,lag]=xcorr(xn,100,'biased')%求XN 的自相关函数R，biased 为有偏估计,lag 为R 的序列号 subplot(3,1,2) hndl=stem(m,r)%绘制离散图，分布点从-100—+100 set(hndl,'Marker','.') set(hndl,'MarkerSize',2)ylabel('自相关函数R(m)') %利用间接法计算功率谱 k=0:1000%取1000 个点 w=(pi/500)*kM=k/500X=r*(exp(-li*pi/500).^(m'*k))%对R 求傅里叶变换,li 就是j magX=abs(X)subplot(3,1,3) plot(M,10*log10(magX))xlabel('功率谱的改进直接法估计')

（3）.自协方差法： 思想：在实际中,有时用随机信号在两个不同时刻t 1 ,t 2 的取值X(t 1 )和X(t 2 )之间的二阶混合中心矩来描述随机信号 X(t)在任意两个时刻取值起伏变化的相依程度, 即自协方差。自协方差函数与自相关函数描述的特性是一致的，所以其原理与相关函数法近似。

clear allfs=20000n=0:1/fs:0.1P=2000*piy=square(P*n)xn=y+0.2*randn(size(n))%绘制信号波形 figure(3) subplot(2,1,1) plot(n,xn) xlabel('时间（s）') ylabel('幅度') title('y+randn(size(n))') ymax_xn=max(xn)+0.2ymin_xn=min(xn)-0.2axis([0 0.3 ymin_xn ymax_xn]) %使用协方差法估计序列功率谱 p=floor(length(xn)/3)+1nfft=1024[xpsd,f]=pcov(xn,p,nfft,fs,'half')%绘制功率谱估计 pmax=max(xpsd)xpsd=xpsd/pmaxxpsd=10*log10(xpsd+0.000001)subplot(2,1,2) plot(f,xpsd) title('基于协方差的功率谱估计') ylabel('功率谱估计（db）') xlabel('频率') grid onymin=min(xpsd)-2ymax=max(xpsd)+2axis([0 fs/2 ymin ymax])

（4）.最大熵法 思想：上网查的最大熵法原理是取一组时间序列，使其自相关函数与一组已知数据的自相关函数相同，同时使已知自相关函数以外的部分的随机性最强，以所取时间序列的谱作为已知数据的谱估值。它等效于根据使随机过程的熵为最大的原则,利用 N 个已知的自相关函数值来外推其他未知的自相关函数值所得到的功率谱。最大熵法功率谱估值是一种可获得高分辨率的非线性谱估值方法，特别适用于数据长度较短的情况。最大熵法谱估值对未知数据的假定 ，一个平稳的随机序列，可以用周期图法对其功率谱进行估值。这种估值方法隐含着假定未知数据是已知数据的周期性重复。现有的线性谱估计方法是假定未知数据的自相关函数值为零，这种人为假定带来的误差较大。最大熵法是利用已知的自相关函数值来外推未知的自相关函数值，去除了对未知数据的人为假定，从而使谱估计的结果更为合理。熵在信息论中是信息的度量，事件越不确定，其信息量越大，熵也越大。对于上述问题来说，对随机过程的未知的自相关函数值，除了从已知的自相关函数值得到有关它的信息以外，没有其他的先验知识。因而，在外推时，不希望加以其他任何新的限制，亦即使之“最不确定”。换言之，就是使随机过程的熵最大。

程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(5) subplot(3,1,1) plot(n,xn)asis([0 0.01 -5 5]) Nfft=256%分段长度256 [Pxx,f]=pmem(xn,14,Nfft,fs)%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2), plot(f,10*log10(Pxx)), title(' 最大熵法，滤波器14'),xlabel('频率'),ylabel('功率谱db')

（5）.最大似然法： 思想：最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。可以证明,如果信号x 是由一个确定性信号S 加上一个高斯白噪声n 所组成，则上述滤波器的输出是信号S 的最大似然估值。如果n 不是高斯噪声，则上述滤波器的输出是信号S 的最小方差的线性的无偏估值，而且它能得到比使用固定的窗口函数的周期图法更高的分辨率。

程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(4) subplot(3,1,1) plot(n,xn)axis([0 0.01 -5 5]) %估计自相关函数 m=-500:500[r,lag]=xcorr(xn,500,'biased')R=[r(501) r(502) r(503) r(504)r(500) r(501) r(502) r(503)r(499) r(500) r(501) r(502)r(498) r(499) r(500) r(501)][V,D]=eig(R)V3=[V(1,3),V(2,3),V(3,3),V(4,3)].'V3=[V(1,4),V(2,4),V(3,4),V(4,4)].'p=0:3wm=[0:0.002*pi:2*pi]B=[(exp(-li)).^(wm'*p)]%li 就是虚数单位j A=B%最小方差功率谱估计 z=A*inv(R)*A'Z=diag(z')pmv=1./Zsubplot(2,1,2) plot(wm/pi,pmv)title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')

自己分下段就能运行了。

---