AR模型更新

look_w

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)
说明：Pxx =Pyulear(x,ORDER,NFFT)中，采用Yule--Walker方法估计序列x的功率谱，参数ORDER用来指定AR模型的阶数，NFFT为FFT算法的长度，默认值为256，若NFFT为偶数，则Pxx为（NFFT/2 + 1）维的列矢量，若NFFT为奇数，则Pxx为（NFFT + 1）/2维的列矢量；当x为复数时，Pxx长度为NFFT。
[Pxx,W]=Pyulear(x,ORDER,NFFT)中，返回一个频率向量W.
[Pxx,W]=Pyulear(x,ORDER,NFFT,Fs)中，可以在F向量得到功率谱估计的频率点，Fs指定采样频率。
Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)中，直接画出功率谱估计的曲线图。
2． Pburg函数：
功能：利用Burg方法进行功率谱估计。
格式：Pxx=Pburg(x,ORDER,NFFT)
[Pxx,W]=Pburg(x,ORDER,NFFT)
[Pxx,W]=Pburg(x,ORDER,NFFT,Fs)
Pburg(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
说明：Pburg函数与Pyulear函数格式相同，只是计算AR模型时所采用的方法不同，因此格式可以参照Pyulear函数。
3． Pcov函数：
功能：利用协方差方法进行功率谱估计。
格式：Pxx=Pcov(x,ORDER,NFFT)
[Pxx,W]=Pcov(x,ORDER,NFFT)
[Pxx,W]=Pcov(x,ORDER,NFFT,Fs)
Pcov(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
说明：Pcov函数采用协方差法估计AR模型的参数，然后计算序列x的功率谱。协方差法与改进的协方差法相比，前者仅令前向预测误差为最小，其他步骤是一样的。：Pcov函数与Pyulear函数格式相同，只是计算AR模型时所采用的方法不同，因此格式可以参照Pyulear函数.
4.Pmcov:
功能：利用改进的协方差方法进行功率谱估计。
格式：Pxx=Pmcov(x,ORDER,NFFT)
[Pxx,W]=Pmcov(x,ORDER,NFFT)
[Pxx,W]=Pmcov(x,ORDER,NFFT,Fs)
Pmcov(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
例如：输入下面语句：
figure 8.10--8.11
Fs=1000; %采样频率
n=0:1/Fs:3;
xn=cos(2*pi*n*200)+randn(size(n));
%设置参数
order=20;
nfft=1024;
%Yule-Walker方法
figure(1)
pyulear(xn,order,nfft,Fs);
%Burg方法
figure(2)
pburg(xn,order,nfft,Fs);
%协方差法
figure(3)
pcov(xn,order,nfft,Fs);
%改进协方差方法
figure(4)
pmcov(xn,order,nfft,Fs);


AR谱的分辨率：
经典谱估计的分辨率反比与信号的有效长度,但是现代谱估计的分辨率可以不受此限制. 这是因为对于给定的N点有限长序列x(n)，虽然其估计出的相关函数也是有限长的，但是现代谱估计的一些方法隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，因而AR谱的分辨率较高。
例如：序列x(n)由两个正铉信号组成，其频率分别为f1=20Hz和f2=21Hz,并含有一定的噪声量。试分别用周期图法，Burg方法与改进的协方差法估计信号的功率谱，且AR模型的阶数取30和50两种情况讨论。
上面的例子可以通过下面程序实现：
Fs=200;
n=0:1/Fs:1;
xn=sin(2*pi*20*n)+sin(2*pi*21*n)+0.1*randn(size(n));
window=boxcar(length(xn));
nfft=512;
[Pxx,f]=periodogram(xn,window,nfft,Fs);
figure(1)
plot(f,10*log10(Pxx)),grid
xlabel('Frequency(Hz)')
ylabel('Power Spectral Density(dB/Hz)')
title('Periodogram PSD Estimate')
order1=30;
order2=50;
figure(2)
pburg(xn,order1,nfft,Fs)
figure(3)
pburg(xn,order2,nfft,Fs)
figure(4)
pmcov(xn,order1,nfft,Fs)
figure(5)
pmcov(xn,order1,nfft)

4. C语言实现

• /*
• * ar_model.h
• *
• *  Created on: 2013-8-11
• *      Author: monkeyzx
• */
• #ifndef AR_MODEL_H_
• #define AR_MODEL_H_
• typedef
  struct {
•     float real;
•     float imag;
• } complex;
• extern
  void maryuwa(complex x[],complex a[],complex r[],int n,int ip,
•         float *ep,int *ierror);
• extern
  void mpsplot(float psdr[],float psdi[],int mfre,float ts);
• extern
  void zx_ar_model(void);
• #endif /* AR_MODEL_H_ */