JSOC_Documentation: proj/globalhs/libs/dtgf/multi_burg_code

00001 // $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/dtgf/multi_burg_code_C99.h,v 1.1 2013/04/28 07:46:58 tplarson Exp $
00002 
00003 #include "ctypes_def.h"
00004 
00005 // type specific definitions
00006 #if TYPE == FLOAT
00007 #define MULTI_BURG smulti_burg
00008 #define EPSILON (10*FLT_EPSILON)   
00009 #elif TYPE == DOUBLE 
00010 #define MULTI_BURG dmulti_burg
00011 #define EPSILON (10*DBL_EPSILON)
00012 #elif TYPE == COMPLEXFLOAT
00013 #define MULTI_BURG cmulti_burg
00014 #define EPSILON (10*FLT_EPSILON)
00015 #elif TYPE == COMPLEXDOUBLE
00016 #define MULTI_BURG zmulti_burg
00017 #define EPSILON (10*DBL_EPSILON)
00018 #endif
00019 
00020 
00021 #ifdef NOBLAS
00022 
00023 int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1], 
00024         RTYPE E[order+1])
00025 {
00026   int k, i, j, skip, N, n_total, effective_order;
00027   CTYPE num, *atmp, *parcor, **ef, **eb, pk, cpk;
00028   RTYPE den, err, err0, err_old, alpha;
00029 
00030   effective_order = order;
00031   skip = 0;
00032   N = 0;
00033   n_total=0;
00034   for (i=0; i<n; i++)
00035   {
00036     n_total += m[i];
00037     if ( m[i] < order ) 
00038       skip++;
00039     else
00040       N += m[i];
00041   }
00042 #ifdef DEBUGOUT
00043   printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
00044 #endif
00045 
00046   if (skip >= n)
00047   {
00048     fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
00049         " Cannot compute AR model.\n",__func__,order);
00050     abort();
00051   }
00052 
00053   parcor = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
00054   atmp = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
00055   ef = (CTYPE **)malloc(n*sizeof(CTYPE*));
00056   eb = (CTYPE **)malloc(n*sizeof(CTYPE*));  
00057   for (i=0; i<n; i++)
00058   {
00059     if ( m[i] >= order ) 
00060     {
00061       ef[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
00062       eb[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
00063       memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
00064       memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
00065     }
00066   }
00067   
00068   memset(a,0,(order+1)*sizeof(CTYPE));
00069   a[0] = CONE;
00070 
00071   // Initial error.
00072   err = 0;
00073   for (i=0; i<n; i++)
00074   {
00075     if ( m[i] >= order ) 
00076     {
00077       for (j=0; j<m[i]; j++)
00078     err += SQR(x[i][j]); 
00079     }
00080   }
00081   err /= N; 
00082   err0 = err;
00083 
00084   if (E != NULL)
00085     E[0] = err;
00086 
00087   for (k=1; k<=order; k++)
00088   {
00089     num = CZERO;
00090     den = RZERO;
00091     for (i=0; i<n; i++)
00092     {
00093       if ( m[i] >= order ) 
00094       {
00095     for (j=0; j<m[i]-k; j++)
00096     {
00097       CTYPE f = ef[i][j+k];
00098       CTYPE b = eb[i][j];
00099       
00100       num -= CONJ(b) * f;
00101       den += SQR(f) + SQR(b);
00102     } 
00103       }
00104     }
00105     if (den == RZERO)
00106     {
00107       fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
00108       effective_order = k-1;
00109       break;
00110     }
00111     parcor[k] = (2*num)/den;
00112     pk = parcor[k];
00113     cpk = CONJ(pk);
00114 
00115     for (i=0; i<n; i++)
00116     {
00117       if ( m[i] >= order ) 
00118       {
00119     for (j=0; j<m[i]-k; j++)
00120     {
00121       CTYPE tmp = ef[i][j+k];     
00122       
00123       ef[i][j+k] +=  pk * eb[i][j];
00124       eb[i][j] += cpk * tmp; 
00125     }
00126       }
00127     }
00128 
00129     alpha =  SQR(pk);
00130     /*    if ( alpha < EPSILON )
00131       {
00132       fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
00133       " %e.\n",k-1,alpha);
00134       effective_order = k-1;
00135       break;
00136       } */
00137     err_old = err;
00138     err *= RONE - alpha;
00139 
00140     memcpy(atmp, a, k*sizeof(CTYPE));
00141     atmp[k] = RZERO;
00142     for (i=1; i<=k; i++)
00143       a[i] = atmp[i] + pk * CONJ(atmp[k-i]);
00144 
00145     if (E != NULL)
00146       E[k] = err;
00147 
00148 
00149     if (err < EPSILON*err0 )
00150     {
00151       fprintf(stderr,"Model order set to %d since prediction error was"
00152           " reduced to %f x variance.\n",k-1,EPSILON);
00153       effective_order = k;
00154       break;
00155     } 
00156   }
00157 
00158   free(parcor);
00159   free(atmp);
00160   for (i=0; i<n; i++)
00161   {
00162     if ( m[i] >= order ) 
00163     {
00164       free(ef[i]);
00165       free(eb[i]);
00166     }
00167   }
00168   free(ef);
00169   free(eb);
00170   return effective_order;
00171 }
00172 
00173 #else
00174 
00175 int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1], 
00176         RTYPE E[order+1])
00177 {
00178   int k, i, j, skip, N, n_total, effective_order;
00179   CTYPE num, *atmp, *parcor, **ef, **eb, pk, cpk;
00180   RTYPE den, err, err0, err_old, alpha;
00181 
00182   effective_order = order;
00183   skip = 0;
00184   N = 0;
00185   n_total=0;
00186   for (i=0; i<n; i++)
00187   {
00188     n_total += m[i];
00189     if ( m[i] < order ) 
00190       skip++;
00191     else
00192       N += m[i];
00193   }
00194 #ifdef DEBUGOUT
00195     printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
00196 #endif
00197 
00198   if (skip >= n)
00199   {
00200     fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
00201         " Cannot compute AR model.\n",__func__,order);
00202     abort();
00203   }
00204 
00205   parcor = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
00206   atmp = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
00207   ef = (CTYPE **)malloc(n*sizeof(CTYPE*));
00208   eb = (CTYPE **)malloc(n*sizeof(CTYPE*));  
00209   for (i=0; i<n; i++)
00210   {
00211     if ( m[i] >= order ) 
00212     {
00213       ef[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
00214       eb[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
00215       memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
00216       memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
00217     }
00218   }
00219   
00220   memset(a,0,(order+1)*sizeof(CTYPE));
00221   a[0] = CONE;
00222 
00223   // Initial error.
00224   err = 0;
00225   for (i=0; i<n; i++)
00226   {
00227     if ( m[i] >= order ) 
00228       err += DOTC(m[i], x[i], 1, x[i],1);
00229   }
00230   err /= N; 
00231   err0 = err;
00232 
00233   if (E != NULL)
00234     E[0] = err;
00235 
00236   for (k=1; k<=order; k++)
00237   {
00238     num = CZERO;
00239     den = RZERO;
00240     for (i=0; i<n; i++)
00241     {
00242       if ( m[i] >= order ) 
00243       {
00244     num -= DOTC(m[i]-k, eb[i], 1, &(ef[i][k]), 1);
00245     den += DOTC(m[i]-k, &(ef[i][k]), 1, &(ef[i][k]), 1) + 
00246       DOTC(m[i]-k, eb[i], 1, eb[i], 1); 
00247     
00248       }
00249     }
00250     if (den == RZERO)
00251     {
00252       fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
00253       effective_order = k-1;
00254       break;
00255     }
00256     parcor[k] = (2*num)/den;
00257     pk = parcor[k];
00258     cpk = CONJ(pk);
00259 
00260     for (i=0; i<n; i++)
00261     {
00262       if ( m[i] >= order ) 
00263       {
00264     for (j=0; j<m[i]-k; j++)
00265     {
00266       CTYPE tmp = ef[i][j+k];     
00267 
00268       ef[i][j+k] +=  pk * eb[i][j];
00269       eb[i][j] += cpk * tmp; 
00270     }
00271       }
00272     }
00273 
00274     alpha =  SQR(pk);
00275     /*    if ( alpha < EPSILON )
00276       {
00277       fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
00278       " %e.\n",k-1,alpha);
00279       effective_order = k-1;
00280       break;
00281       } */
00282     err_old = err;
00283     err *= RONE - alpha;
00284 
00285     memcpy(atmp, a, k*sizeof(CTYPE));
00286     atmp[k] = RZERO;
00287     for (i=1; i<=k; i++)
00288       a[i] = atmp[i] + pk * CONJ(atmp[k-i]);
00289 
00290     if (E != NULL)
00291       E[k] = err;
00292 
00293 
00294     if (err < EPSILON*err0 )
00295     {
00296       fprintf(stderr,"Model order set to %d since prediction error was"
00297           " reduced to %f x variance.\n",k-1,EPSILON);
00298       effective_order = k;
00299       break;
00300     } 
00301   }
00302 
00303   free(parcor);
00304   free(atmp);
00305   for (i=0; i<n; i++)
00306   {
00307     if ( m[i] >= order ) 
00308     {
00309       free(ef[i]);
00310       free(eb[i]);
00311     }
00312   }
00313   free(ef);
00314   free(eb);
00315   return effective_order;
00316 }
00317 
00318 #endif /* NOBLAS */
00319 
00320 #undef DOTC
00321 #undef MULTI_BURG
00322 #undef EPSILON
00323 #undef SQR
00324 #include "ctypes_undef.h"
00325