/* Previously called mempar() Originally in FORTRAN, hence the array offsets of 1, Yuk. Original code from Kay, 1988, appendix 8D. Perform Burg's Maximum Entropy AR parameter estimation outputting (or not) successive models en passant. Sourced from Alex Sergejew Two small changes made by NH in November 1998: tstarz.h no longer included, just say "typedef double REAL" instead Declare ar by "REAL **ar" instead of "REAL ar[MAXA][MAXA] Further "cleaning" by Paul Bourke.....for personal style only. Converted to zero-based arrays by Paul Sanders, June 2007 */ #define GENERATE_INTERMEDIATES void ARMaxEntropy (double *inputseries, int length, int degree, #ifdef GENERATE_INTERMEDIATES double **ar, #else double *coef, #endif double *per, double *pef, double *h, double *g) { double t1, t2; int n; for (n = 1; n <= degree; n++) { double sn = 0.0; double sd = 0.0; int j; int jj = length - n; for (j = 0; j < jj; j++) { t1 = inputseries [j + n] + pef [j]; t2 = inputseries [j] + per [j]; sn -= 2.0 * t1 * t2; sd += (t1 * t1) + (t2 * t2); } t1 = g [n] = sn / sd; if (n != 1) { for (j = 1; j < n; j++) h [j] = g [j] + t1 * g [n - j]; for (j = 1; j < n; j++) g [j] = h [j]; jj--; } for (j = 0; j < jj; j++) { per [j] += t1 * pef [j] + t1 * inputseries [j + n]; pef [j] = pef [j + 1] + t1 * per [j + 1] + t1 * inputseries [j + 1]; } #ifdef GENERATE_INTERMEDIATES for (j = 0; j < n; j++) ar [n][j] = g [j + 1]; #endif } #ifndef GENERATE_INTERMEDIATES for (n = 0; n < degree; n++) coef [n] = g [n + 1]; #endif }