周期图法与ar模型.doc

周期图法：周期图法是直接将信号的采样数据x(n)进行傅立叶变换来求取功率谱密度估计。假设有限长随机信号序列x(n),它的傅立叶变换和功率谱密度估计S存在下列关系： Sn=1Nx(n)2式中，N为随机信号x(n)的长度，在离散的频率点处可得： Sk=1Nx(k)2= 1NDFTxn2, k=0,1,N-1由于DFTx(n)的周期为N，求得的功率谱估计也是以N为周期，故成为周期图法。AR模型：信号的AR模型定义为： xn=-k=1pakxn-k+u(n) （1） Hz=1A(z)=11+k=1pakz-k （2） Px(ej)=21+k=1pake-jk2 （3）式中x(n)为信号X(k)的下一点的估计，因此，可以用AR模型来进行预报。式中中P为AR模型阶次，一般可以用模型拟合残差曲线来确定。当模型阶次增加，残差不在显著下降时，该阶次就是正确的模型阶次。 另外，AR模型是信号的参数化模型，可以从AR模型系数计算出信号的频谱，称为最大墒谱，公式如下： 实际工作中，往往只是知道N点数据x(n); 因此，可以按如下步骤估计x(n)的功率谱：（1）、由x(n)估计自相关函数，得r(m);（2）、用r(m)根据Yule-Wslker方程求出AR模型真实参数的估计值，即a(1), a(2), , ;（3）、将这些参数带入上式（3）中，得到x(n)的功率谱估计程序如下：1， 周期图法load y.txt;load x.txt;data=y;M=100;PS=zeros(M,1);vecter=myditfft(data);u=(1/M)*sum(data);for n=1:M PS(n)=(1/M)*(vecter(n)*conj(vecter(n)+u;endplot(xnian,PS); function y=myditfft(x)%FFT变换，实现DIT-FFT基2算法，%x为给定序列，y为x的离散傅立叶变换%点数取大于等于x的长度的2的幂次m=nextpow2(x);N=2m; if length(x)N x=x,zeros(1,N-length(x); end nxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1;y=x(nxd); for mm=1:m Nmr=2mm;u=1; WN=exp(-i*2*pi/Nmr);)for j=1:Nmr/2 for k=j:Nmr:N kp=k+Nmr/2; t=y(kp)*u; y(kp)=y(k)-t; y(k)=y(k)+t; endu=u*WN;endend2、AR模型load x,txt;load y.txt;data_burg=y; N=length(y);k=input(k=); A=zeros(N,1);PSD=zeros(1,N);r,a,sigma=burg_unknown(data_burg,k);s=sigma(k+1);aa=1 a(k,1:k) zeros(1,N-k-2);aa_fft=myditfft(aa);for l=1:N A(l)=aa_fft(l)*conj(aa_fft(l); PSD(l)=s/A(l);endplot(xnian,PSD);function r,a,sigma=burg_unknown(x,K)N=max(size(x);r=;ef(1,:)=x;eb(1,:)=x;sigma(1)=sum(x.2)/N;a=zeros(K+1,K+1);a(:,1)=1;for p=2:K+1; r(p-1)=2*(ef(p-1,p:N)*eb(p-1,p-1:N-1)/(ef(p-1,p:N)*ef(p-1,p:N)+eb(p-1,p-1:N-1)*eb(p-1,p-1:N-1); sigma(p)=(1-r(p-1).2)*sigma(p-1); for i=2:p-1 a(p,i)=a(p-1,i)-r(p-1)*a(p-1,p-i+1); end a(p,p)=-r(p-1); for n=2:N ef(p,n)=ef(p-1,