JSOC_Documentation: proj/libs/gapfiller/multi_burg_code

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