AR模型

转自:http://blog.csdn.net/xiahouzuoxin/article/details/9904147


本文目标:分析AR模型并求解AR模型的输出x(n)的功率谱。


1. AR模型概念观

数字信号处理功率谱估计方法分经典功率谱估计和现代功率谱估计,现代功率谱估计以参数模型功率谱估计为代表,参数功率谱模型如下:             

u(n) ——>  H(z)   ——> x(n)

参数模型的基本思路是:

—— 参数模型假设研究过程是由一个输入序列u(n)激励一个线性系统H(z)的输出。

—— 由假设参数模型的输出x(n)或其自相关函数来估计H(z)的参数

—— 由H(z)的参数估计x(n)的功率谱

因此,参数模型功率谱的求解有两步:

(1)H(z)模型参数估计

(2)依据模型参数求功率谱

AR模型(自回归模型,Auto Regression Model)是典型的现代参数功模型。其定义为

其中,输入设定为方差为的白噪声序列,ak是模型的参数,p是模型的阶数,Px为x(n)功率谱,也即本文要求解的目标。

AR模型是一个全极点模型,“自回归”的含义是:现在的输出是现在的输入和过去p个输出的加权和。

现在我们希望建立AR参数模型和x(n)的自相关函数的关系,也即AR模型的正则方程:


上面的正则方程也称Yule-Walker方程,其中的rx为自相关函数。由方程可以看出,一个p阶的AR模型有p+1个参数()。

通过推导可以发现,AR模型与线性预测器是等价的,AR模型是在最小平方意义上对数据的拟合。

2. AR模型参数求解——Levinson-Durbin Algorithm

定义为p阶AR模型在m阶次时的第k个系数,k=1,2,...,m。定义为m阶系统时的,这也是线性预测器中前向预测的最小误差功率。此时, 一阶AR模型时有



我们定义初始时,则


由PART1中矩阵的对称性质,将上面的公式推广到高阶AR模型,可以推导出Levinson-Durbin递推算法:




Levinson-Durbin递推算法从低阶开始递推,,给出了每一阶次时所有参数,。这一特点有利于我们选择合适的AR模型阶次。

因为必须大于0,由式,如果,递推应该停止。

到此,选择最佳阶次的参数代入到中,求得功率谱。


3. matlab实现

matlab工具箱中提供了现成的函数实现AR模型功率谱计算。参考 [2] ,我们将内容摘录如下:


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)

说明: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语言实现

  1. /* 
  2.  * ar_model.h 
  3.  * 
  4.  *  Created on: 2013-8-11 
  5.  *      Author: monkeyzx 
  6.  */  
  7.   
  8. #ifndef AR_MODEL_H_  
  9. #define AR_MODEL_H_  
  10.   
  11. typedef struct {  
  12.     float real;  
  13.     float imag;  
  14. } complex;  
  15.   
  16. extern void maryuwa(complex x[],complex a[],complex r[],int n,int ip,  
  17.         float *ep,int *ierror);  
  18. extern void mpsplot(float psdr[],float psdi[],int mfre,float ts);  
  19.   
  20. extern void zx_ar_model(void);  
  21.   
  22. #endif /* AR_MODEL_H_ */  

  1. /* 
  2.  * ar_model.c 
  3.  * 
  4.  *  Created on: 2013-8-11 
  5.  *      Author: monkeyzx 
  6.  */  
  7.   
  8. #include <stdio.h>  
  9. #include <stdlib.h>  
  10. #include <math.h>  
  11. #include <stdlib.h>  
  12. //#include "msp.h"  
  13. #include "ar_model.h"  
  14. #include "time.h"  
  15.   
  16. float mabs(complex a)  
  17. {  
  18.      float m;  
  19.   
  20.      m=a.real*a.real+a.imag*a.imag;  
  21.      m=sqrt(m);  
  22.   
  23.      return(m);  
  24. }  
  25.   
  26. /*--------------------------------------------------------------------- 
  27.   Routine MCORRE1:To estimate the biased cross-correlation function 
  28.   of complex arrays x and y. If y=x,then it is auto-correlation. 
  29.   input parameters: 
  30.      x  :n dimensioned complex array. 
  31.      y  :n dimensioned complex array. 
  32.      n  :the dimension of x and y. 
  33.      lag:point numbers of correlation. 
  34.   output parameters: 
  35.      r  :lag dimensioned complex array, the correlation function is 
  36.          stored in r(0) to r(lag-1). 
  37.                                       in Chapter 1 and 11 
  38. ---------------------------------------------------------------------*/  
  39. void mcorre1(complex x[],complex y[],complex r[],int n,int lag)  
  40. {  
  41.     int m,j,k;  
  42.   
  43.     for(k=0;k<lag;k++) {  
  44.         m=n-1-k;  
  45.         r[k].real=0.0f;  
  46.         r[k].imag=0.0f;  
  47.         for(j=0;j<=m;j++) {  
  48.             r[k].real+=y[j+k].real*x[j].real+y[j+k].imag*x[j].imag;  
  49.             r[k].imag+=y[j+k].imag*x[j].real-y[j+k].real*x[j].imag;  
  50.         }  
  51.         r[k].real=r[k].real/n;  
  52.         r[k].imag=r[k].imag/n;  
  53.     }  
  54.     return;  
  55. }  
  56.   
  57. /*--------------------------------------------------------------------- 
  58.   Routine maryuwa: To determine the autoregressive coefficients by 
  59.           solving Yule-Walker equation with Levinson algorithm. 
  60.   Input Parameters: 
  61.      n     : Number of data samples (integer) 
  62.      ip    : Order of autoregressive model 
  63.      x     : Array of complex data values, x(0) to x(n-1) 
  64.   Output Parameters: 
  65.      ep    : Driving noise variance (real) 
  66.      a     : Array of complex autoregressive coefficients, a(0) to 
  67.              a(ip) 
  68.   ierror=0 : No error 
  69.         =1 : ep<=0 . 
  70.  
  71.         r  : complex work array, auto-correlation 
  72.                                        in chapter 12 
  73. --------------------------------------------------------------------*/  
  74. void maryuwa(complex x[],complex a[],complex r[],int n,int ip,  
  75. float *ep,int *ierror)  
  76. {  
  77.     complex sum;  
  78.     int i,k;  
  79.     float r0;  
  80.   
  81.     *ierror=1;  
  82.     mcorre1(x,x,r,n,ip+1);  
  83.     a[0].real=1.0;  
  84.     a[0].imag=0.0;  
  85.     r0=r[0].real;  
  86.     a[1].real=-r[1].real/r0;  
  87.     a[1].imag=-r[1].imag/r0;  
  88.     *ep=r0*(1.0f-pow(mabs(a[1]),2));  
  89.     for(k=2;k<=ip;k++) {  
  90.         sum.real=0.;  
  91.         sum.imag=0.;  
  92.         for(i=1;i<k;i++) {  
  93.             sum.real+=r[k-i].real*a[i].real-r[k-i].imag*a[i].imag;  
  94.             sum.imag+=r[k-i].real*a[i].imag+r[k-i].imag*a[i].real;  
  95.         }  
  96.         sum.real+=r[k].real;  
  97.         sum.imag+=r[k].imag;  
  98.         a[k].real=-sum.real/(*ep);  
  99.         a[k].imag=-sum.imag/(*ep);  
  100.         (*ep)*=1.-pow(mabs(a[k]),2);  
  101.         if(*ep<=0.0)  
  102.             return;  
  103.         for(i=1;i<k;i++) {  
  104.             x[i].real=a[i].real+a[k-i].real*a[k].real+  
  105.                     a[k-i].imag*a[k].imag;  
  106.             x[i].imag=a[i].imag+a[k-i].real*a[k].imag-  
  107.                     a[k-i].imag*a[k].real;  
  108.         }  
  109.         for(i=1;i<k;i++) {  
  110.             a[i].real=x[i].real;  
  111.             a[i].imag=x[i].imag;  
  112.         }  
  113.     }  
  114.     *ierror=0;  
  115. }  
  116.   
  117. /*---------------------------------------------------------------------- 
  118.   routinue mrelfft:To perform  split-radix DIF fft algorithm. 
  119.  
  120.   input parameters: 
  121.    xr,xi:real and image part of complex data for DFT/IDFT,n=0,...,N-1 
  122.    N    :Data point number of DFT compute . 
  123.    isign:Transform direction disignator , 
  124.                isign=-1: For Forward Transform. 
  125.                isign=+1: For Inverse Transform. 
  126.  
  127.   output parameters: 
  128.    xr,xi:real and image part of complex result of DFT/IDFT,n=0,...,N-1 
  129.  
  130.   Note: N  must be a power of 2 . 
  131.                                        in chapter 5 
  132. ---------------------------------------------------------------------*/  
  133. void mrelfft(float xr[],float xi[],int n,int isign)  
  134. {  
  135.     float e,es,cc1,ss1,cc3,ss3,r1,s1,r2,s2,s3,xtr,xti,a,a3;  
  136.     int m,n2,n4,j,k,is,id,i0,i1,i2,i3,n1,i,nn;  
  137.   
  138.     for(m=1;m<=16;m++) {  
  139.         nn=pow(2,m);  
  140.         if(n==nn)break;  
  141.     }  
  142.     if(m>16) {  
  143. #ifdef _DEBUG  
  144.         printf(" N is not a power of 2 ! \n");  
  145. #endif  
  146.         return;  
  147.     }  
  148.     n2=n*2;  
  149.     es=-isign*atan(1.0)*8.0;  
  150.     for(k=1;k<m;k++) {  
  151.         n2=n2/2;  
  152.         n4=n2/4;  
  153.         e=es/n2;  
  154.         a=0.0;  
  155.         for(j=0;j<n4;j++) {  
  156.             a3=3*a;  
  157.             cc1=cos(a);  
  158.             ss1=sin(a);  
  159.             cc3=cos(a3);  
  160.             ss3=sin(a3);  
  161.             a=(j+1)*e;  
  162.             is=j;  
  163.             id=2*n2;  
  164.             do {  
  165.                 for(i0=is;i0<n;i0+=id) {  
  166.                     i1=i0+n4;  
  167.                     i2=i1+n4;  
  168.                     i3=i2+n4;  
  169.                     r1=xr[i0]-xr[i2];  
  170.                     s1=xi[i0]-xi[i2];  
  171.                     r2=xr[i1]-xr[i3];  
  172.                     s2=xi[i1]-xi[i3];  
  173.                     xr[i0]+=xr[i2];  
  174.                     xi[i0]+=xi[i2];  
  175.                     xr[i1]+=xr[i3];  
  176.                     xi[i1]+=xi[i3];  
  177.                     if(isign!=1) {  
  178.                         s3=r1-s2;  
  179.                         r1=r1+s2;  
  180.                         s2=r2-s1;  
  181.                         r2=r2+s1;  
  182.                     } else {  
  183.                             s3=r1+s2;  
  184.                             r1=r1-s2;  
  185.                             s2=-r2-s1;  
  186.                             r2=-r2+s1;  
  187.                     }  
  188.                     xr[i2]=r1*cc1-s2*ss1;  
  189.                     xi[i2]=-s2*cc1-r1*ss1;  
  190.                     xr[i3]=s3*cc3+r2*ss3;  
  191.                     xi[i3]=r2*cc3-s3*ss3;  
  192.                 }  
  193.                 is=2*id-n2+j;  
  194.                 id=4*id;  
  195.             }while(is<n-1);  
  196.         }  
  197.     }  
  198. /*   ------------ special last stage -------------------------*/  
  199.     is=0;  
  200.     id=4;  
  201.     do {  
  202.         for(i0=is;i0<n;i0+=id) {  
  203.             i1=i0+1;  
  204.             xtr=xr[i0];  
  205.             xti=xi[i0];  
  206.             xr[i0]=xtr+xr[i1];  
  207.             xi[i0]=xti+xi[i1];  
  208.             xr[i1]=xtr-xr[i1];  
  209.             xi[i1]=xti-xi[i1];  
  210.         }  
  211.         is=2*id-2;  
  212.         id=4*id;  
  213.     } while(is<n-1);  
  214.     j=1;  
  215.     n1=n-1;  
  216.     for(i=1;i<=n1;i++) {  
  217.         if(i<j) {  
  218.             xtr=xr[j-1];  
  219.             xti=xi[j-1];  
  220.             xr[j-1]=xr[i-1];  
  221.             xi[j-1]=xi[i-1];  
  222.             xr[i-1]=xtr;  
  223.             xi[i-1]=xti;  
  224.         }  
  225.         k=n/2;  
  226.         while(1) {  
  227.             if(k>=j)break;  
  228.             j=j-k;  
  229.             k=k/2;  
  230.         }  
  231.         j=j+k;  
  232.     }  
  233.     if(isign==-1) return;  
  234.     for(i=0;i<n;i++) {  
  235.         xr[i]/=n;  
  236.         xi[i]/=n;  
  237.     }  
  238. }  
  239.   
  240. /*--------------------------------------------------------------------- 
  241.    Routine mpsplot: To plot the normalized power spectum curve on the 
  242.    normalized frequency axis from -.5 to  +.5 . 
  243.         mfre : Points in frequency axis and must be the power of 2. 
  244.         ts   : Sample interval in seconds (real). 
  245.         psdr : Real array of power spectral density values. 
  246.         psdi : Real work array. 
  247.                                        in chapter 11,12 
  248. --------------------------------------------------------------------*/  
  249. void mpsplot(float psdr[],float psdi[],int mfre,float ts)  
  250. {  
  251.     FILE *fp;  
  252.     char filename[30];  
  253.     int k,m2;  
  254.     float pmax,fs,faxis;  
  255.   
  256.     m2=mfre/2;  
  257.     for(k=0;k<m2;k++){  
  258.         psdi[k]=psdr[k];  
  259.         psdr[k]=psdr[k+m2];  
  260.         psdr[k+m2]=psdi[k];  
  261.     }  
  262.     pmax=psdr[0];  
  263.     for(k=1;k<mfre;k++)  
  264.         if(psdr[k]>pmax)  
  265.             pmax=psdr[k];  
  266.         for(k=0;k<mfre;k++) {  
  267.             psdr[k]=psdr[k]/pmax;  
  268.         if(psdr[k]<=0.0)  
  269.             psdr[k]=.000001;  
  270.     }  
  271.     fs=1./ts;  
  272.     fs=fs/(float)(mfre);  
  273.     printf("Please input filename:\n");  
  274.     scanf("%s",filename);  
  275.     if((fp=fopen(filename,"w"))==NULL) {  
  276.         printf("cannot open file\n");  
  277.         exit(0);  
  278.     }  
  279.     for(k=0;k<mfre;k++) {  
  280.         faxis=fs*(k-m2);  
  281.         fprintf(fp,"%f,%f\n",faxis,10.*log10(psdr[k]));  
  282.     }  
  283.     fclose(fp);  
  284.     return;  
  285. }  
  286.   
  287. /*---------------------------------------------------------------------- 
  288.    Routine mar1psd: To compute the power spectum by AR-model parameters. 
  289.    Input parameters: 
  290.           ip : AR model order (integer) 
  291.           ep   : White noise variance of model input (real) 
  292.           ts   : Sample interval in seconds (real) 
  293.           a    : Complex array of AR  parameters a(0) to a(ip) 
  294.    Output parameters: 
  295.           psdr : Real array of power spectral density values 
  296.           psdi : Real work array 
  297.                                         in chapter 12 
  298. ---------------------------------------------------------------------*/  
  299. void mar1psd(complex a[],int ip,int mfre,float *ep,float ts)  
  300. {  
  301.     static float psdr[4096];  
  302.     static float psdi[4096];  
  303.     int k;  
  304.     float p;  
  305.   
  306.     for(k=0;k<=ip;k++) {  
  307.         psdr[k]=a[k].real;  
  308.         psdi[k]=a[k].imag;  
  309.     }  
  310.     for(k=ip+1;k<mfre;k++) {  
  311.         psdr[k]=0.;  
  312.         psdi[k]=0.;  
  313.     }  
  314.     mrelfft(psdr,psdi,mfre,-1);  
  315.     for(k=0;k<mfre;k++) {  
  316.         p=pow(psdr[k],2)+pow(psdi[k],2);  
  317.         psdr[k]=(*ep)*ts/p;  
  318.     }  
  319.   
  320.     mpsplot(psdr,psdi,mfre,ts);  
  321.   
  322.     return;  
  323. }  
  324.   
  325.   
  326. /* 
  327.  * Below are examples for using @maryuwa and @mar1psd 
  328.  */  
  329. #define PI            (3.1415926)  
  330. #define N             (1024)  
  331. #define AN            (10)  
  332. complex x[N];  
  333. complex r[N];  
  334. complex a[AN];  
  335.   
  336. /* 
  337.  * generate random number which satify guass distribution 
  338.  */  
  339. double guass_rand(void)  
  340. {  
  341.     static double V1, V2, S;  
  342.     static int phase = 0;  
  343.     double X;  
  344.   
  345.     if ( phase == 0 ) {  
  346.         do {  
  347.             double U1 = (double)rand() / RAND_MAX;  
  348.             double U2 = (double)rand() / RAND_MAX;  
  349.   
  350.             V1 = 2 * U1 - 1;  
  351.             V2 = 2 * U2 - 1;  
  352.             S = V1 * V1 + V2 * V2;  
  353.         } while(S >= 1 || S == 0);  
  354.   
  355.         X = V1 * sqrt(-2 * log(S) / S);  
  356.     } else {  
  357.         X = V2 * sqrt(-2 * log(S) / S);  
  358.     }  
  359.   
  360.     phase = 1 - phase;  
  361.   
  362.     return X;  
  363. }  
  364.   
  365. void zx_ar_model(void)  
  366. {  
  367.     int i=0;  
  368.     float ep = 0;  
  369.     int ierror = 0;  
  370.   
  371.     /* 
  372.      * generate x[N] 
  373.      */  
  374.     srand(time(NULL));  
  375.     for (i=0; i<N; i++) {  
  376.         x[i].real = sin(2*PI*i/N) + guass_rand();  
  377.         x[i].imag = 0;  
  378.     }  
  379.   
  380.     /* Find parameters for AR model */  
  381.     maryuwa(x, a, r, N, AN, &ep, &ierror);  
  382.   
  383.     /* Calculate power spectum using parameters of AR model */  
  384.     mar1psd(a, AN, N, &ep, 1);  
  385. }  

  1. /* 
  2.  * main.c 
  3.  * 
  4.  *  Created on: 2013-8-11 
  5.  *      Author: monkeyzx 
  6.  */  
  7. #include "ar_model.h"  
  8.   
  9. int main(void)  
  10. {  
  11.     zx_ar_model();  
  12.   
  13.     return 0;  
  14. }  

上面的实例中给定输入信号为余弦信号,采样点数为1024个点,通过计算后的功率谱通过mpsplot函数保存到文本文件output.txt中,保存格式如下:

-0.500000,-15.334630
-0.499023,-15.334833
-0.498047,-15.335444
-0.497070,-15.336456
-0.496094,-15.337864
-0.495117,-15.339655
-0.494141,-15.341816
-0.493164,-15.344331
-0.492188,-15.347179
-0.491211,-15.350342
-0.490234,-15.353794
-0.489258,-15.357505
-0.488281,-15.361453
-0.487305,-15.365603
-0.486328,-15.369924
-0.485352,-15.374381
......

最后借助matlab读取该文件,绘制出功率谱的图形

[plain]  view plain copy print ? 在CODE上查看代码片 派生到我的代码片
  1. data = load('output.txt');  
  2. plot(data(:,1),data(:,2));  


关于上面的C程序,这里只提与主题无关的,double guass_rand(void)是C语言中典型的生成高斯分布随机数的发生器,这里用于在余弦函数上加上一个高斯的噪声。关于更多的随机数生成器可参考关于怎样产生随机数的彻底研究 [自行理解],我将该博文转载过来,感谢作者。


Refrences:

[1] 胡广书《数字信号处理——理论、算法与实现 第二版》

[2] AR模型matlab相关函数描述http://blog.sina.com.cn/s/blog_62f573ad0100sfh1.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值