AR模型更新(3)

look_w

•    Routine mpsplot: To plot the normalized power spectum curve on the
•    normalized frequency axis from -.5 to  +.5 .
•         mfre : Points in frequency axis and must be the power of 2.
•         ts   : Sample interval in seconds (real).
•         psdr : Real array of power spectral density values.
•         psdi : Real work array.
•                                        in chapter 11,12
• --------------------------------------------------------------------*/
  
• void mpsplot(float psdr[],float psdi[],int mfre,float ts)  
• {  
•     FILE *fp;  
•     char filename[30];  
•     int k,m2;  
•     float pmax,fs,faxis;  
• 
  
•     m2=mfre/2;  
•     for(k=0;k<m2;k++){  
•         psdi[k]=psdr[k];  
•         psdr[k]=psdr[k+m2];  
•         psdr[k+m2]=psdi[k];  
•     }  
•     pmax=psdr[0];  
•     for(k=1;k<mfre;k++)  
•         if(psdr[k]>pmax)  
•             pmax=psdr[k];  
•         for(k=0;k<mfre;k++) {  
•             psdr[k]=psdr[k]/pmax;  
•         if(psdr[k]<=0.0)  
•             psdr[k]=.000001;  
•     }  
•     fs=1./ts;  
•     fs=fs/(float)(mfre);  
•     printf("Please input filename:\n");  
•     scanf("%s",filename);  
•     if((fp=fopen(filename,"w"))==NULL) {  
•         printf("cannot open file\n");  
•         exit(0);  
•     }  
•     for(k=0;k<mfre;k++) {  
•         faxis=fs*(k-m2);  
•         fprintf(fp,"%f,%f\n",faxis,10.*log10(psdr[k]));  
•     }  
•     fclose(fp);  
•     return;  
• }  
• 
  
• /*----------------------------------------------------------------------
•    Routine mar1psd: To compute the power spectum by AR-model parameters.
•    Input parameters:
•           ip : AR model order (integer)
•           ep   : White noise variance of model input (real)
•           ts   : Sample interval in seconds (real)
•           a    : Complex array of AR  parameters a(0) to a(ip)
•    Output parameters:
•           psdr : Real array of power spectral density values
•           psdi : Real work array
•                                         in chapter 12
• ---------------------------------------------------------------------*/
  
• void mar1psd(complex a[],int ip,int mfre,float *ep,float ts)  
• {  
•     static
  float psdr[4096];  
•     static
  float psdi[4096];  
•     int k;  
•     float p;  
• 
  
•     for(k=0;k<=ip;k++) {  
•         psdr[k]=a[k].real;  
•         psdi[k]=a[k].imag;  
•     }  
•     for(k=ip+1;k<mfre;k++) {  
•         psdr[k]=0.;  
•         psdi[k]=0.;  
•     }  
•     mrelfft(psdr,psdi,mfre,-1);  
•     for(k=0;k<mfre;k++) {  
•         p=pow(psdr[k],2)+pow(psdi[k],2);  
•         psdr[k]=(*ep)*ts/p;  
•     }  
• 
  
•     mpsplot(psdr,psdi,mfre,ts);  
• 
  
•     return;  
• }  
• 
  
• 
  
• /*
• * Below are examples for using @maryuwa and @mar1psd
• */
  
• #define PI            (3.1415926)
  
• #define N             (1024)
  
• #define AN            (10)
  
• complex x[N];  
• complex r[N];  
• complex a[AN];  
• 
  
• /*
• * generate random number which satify guass distribution
• */
  
• double guass_rand(void)  
• {  
•     static
  double V1, V2, S;  
•     static
  int phase = 0;  
•     double X;  
• 
  
•     if ( phase == 0 ) {  
•         do {  
•             double U1 = (double)rand() / RAND_MAX;  
•             double U2 = (double)rand() / RAND_MAX;  
• 
  
•             V1 = 2 * U1 - 1;  
•             V2 = 2 * U2 - 1;  
•             S = V1 * V1 + V2 * V2;  
•         } while(S >= 1 || S == 0);  
• 
  
•         X = V1 * sqrt(-2 * log(S) / S);  
•     } else {  
•         X = V2 * sqrt(-2 * log(S) / S);  
•     }  
• 
  
•     phase = 1 - phase;  
• 
  
•     return X;  
• }  
• 
  
• void zx_ar_model(void)  
• {  
•     int i=0;  
•     float ep = 0;  
•     int ierror = 0;  
• 
  
•     /*
•      * generate x[N]
•      */
  
•     srand(time(NULL));  
•     for (i=0; i<N; i++) {  
•         x.real = sin(2*PI*i/N) + guass_rand();  
•         x.imag = 0;  
•     }  
• 
  
•     /* Find parameters for AR model */
  
•     maryuwa(x, a, r, N, AN, &ep, &ierror);  
• 
  
•     /* Calculate power spectum using parameters of AR model */
  
•     mar1psd(a, AN, N, &ep, 1);  
• }  

• /*
• * main.c
• *
• *  Created on: 2013-8-11
• *      Author: monkeyzx
• */
  
• #include "ar_model.h"
  
• 
  
• int main(void)  
• {  
•     zx_ar_model();  
• 
  
•     return 0;  
• }  

上面的实例中给定输入信号为余弦信号，采样点数为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读取该文件，绘制出功率谱的图形

• data = load('output.txt');  
• plot(data(:,1),data(:,2));