libs/dtgf/multi_burg_code_C99.h

// $Header: $

#include "ctypes_def.h"

// type specific definitions
#if TYPE == FLOAT
#define MULTI_BURG smulti_burg
#define EPSILON (10*FLT_EPSILON)   
#elif TYPE == DOUBLE 
#define MULTI_BURG dmulti_burg
#define EPSILON (10*DBL_EPSILON)
#elif TYPE == COMPLEXFLOAT
#define MULTI_BURG cmulti_burg
#define EPSILON (10*FLT_EPSILON)
#elif TYPE == COMPLEXDOUBLE
#define MULTI_BURG zmulti_burg
#define EPSILON (10*DBL_EPSILON)
#endif


#ifdef NOBLAS

int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1], 
                RTYPE E[order+1])
{
  int k, i, j, skip, N, n_total, effective_order;
  CTYPE num, *atmp, *parcor, **ef, **eb, pk, cpk;
  RTYPE den, err, err0, err_old, alpha;

  effective_order = order;
  skip = 0;
  N = 0;
  n_total=0;
  for (i=0; i<n; i++)
  {
    n_total += m[i];
    if ( m[i] < order ) 
      skip++;
    else
      N += m[i];
  }
#ifdef DEBUGOUT
  printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
#endif

  if (skip >= n)
  {
    fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
            " Cannot compute AR model.\n",__func__,order);
    abort();
  }

  parcor = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
  atmp = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
  ef = (CTYPE **)malloc(n*sizeof(CTYPE*));
  eb = (CTYPE **)malloc(n*sizeof(CTYPE*));  
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
    {
      ef[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
      eb[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
      memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
      memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
    }
  }
  
  memset(a,0,(order+1)*sizeof(CTYPE));
  a[0] = CONE;

  // Initial error.
  err = 0;
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
    {
      for (j=0; j<m[i]; j++)
        err += SQR(x[i][j]); 
    }
  }
  err /= N; 
  err0 = err;

  if (E != NULL)
    E[0] = err;

  for (k=1; k<=order; k++)
  {
    num = CZERO;
    den = RZERO;
    for (i=0; i<n; i++)
    {
      if ( m[i] >= order ) 
      {
        for (j=0; j<m[i]-k; j++)
        {
          CTYPE f = ef[i][j+k];
          CTYPE b = eb[i][j];
          
          num -= CONJ(b) * f;
          den += SQR(f) + SQR(b);
        } 
      }
    }
    if (den == RZERO)
    {
      fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
      effective_order = k-1;
      break;
    }
    parcor[k] = (2*num)/den;
    pk = parcor[k];
    cpk = CONJ(pk);

    for (i=0; i<n; i++)
    {
      if ( m[i] >= order ) 
      {
        for (j=0; j<m[i]-k; j++)
        {
          CTYPE tmp = ef[i][j+k];         
          
          ef[i][j+k] +=  pk * eb[i][j];
          eb[i][j] += cpk * tmp; 
        }
      }
    }

    alpha =  SQR(pk);
    /*    if ( alpha < EPSILON )
          {
          fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
          " %e.\n",k-1,alpha);
          effective_order = k-1;
          break;
          } */
    err_old = err;
    err *= RONE - alpha;

    memcpy(atmp, a, k*sizeof(CTYPE));
    atmp[k] = RZERO;
    for (i=1; i<=k; i++)
      a[i] = atmp[i] + pk * CONJ(atmp[k-i]);

    if (E != NULL)
      E[k] = err;


    if (err < EPSILON*err0 )
    {
      fprintf(stderr,"Model order set to %d since prediction error was"
              " reduced to %f x variance.\n",k-1,EPSILON);
      effective_order = k;
      break;
    } 
  }

  free(parcor);
  free(atmp);
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
    {
      free(ef[i]);
      free(eb[i]);
    }
  }
  free(ef);
  free(eb);
  return effective_order;
}

#else

int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1], 
                RTYPE E[order+1])
{
  int k, i, j, skip, N, n_total, effective_order;
  CTYPE num, *atmp, *parcor, **ef, **eb, pk, cpk;
  RTYPE den, err, err0, err_old, alpha;

  effective_order = order;
  skip = 0;
  N = 0;
  n_total=0;
  for (i=0; i<n; i++)
  {
    n_total += m[i];
    if ( m[i] < order ) 
      skip++;
    else
      N += m[i];
  }
#ifdef DEBUGOUT
    printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
#endif

  if (skip >= n)
  {
    fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
            " Cannot compute AR model.\n",__func__,order);
    abort();
  }

  parcor = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
  atmp = (CTYPE *)malloc((order+1)*sizeof(CTYPE));
  ef = (CTYPE **)malloc(n*sizeof(CTYPE*));
  eb = (CTYPE **)malloc(n*sizeof(CTYPE*));  
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
    {
      ef[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
      eb[i] = (CTYPE *)malloc(m[i]*sizeof(CTYPE));
      memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
      memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
    }
  }
  
  memset(a,0,(order+1)*sizeof(CTYPE));
  a[0] = CONE;

  // Initial error.
  err = 0;
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
      err += DOTC(m[i], x[i], 1, x[i],1);
  }
  err /= N; 
  err0 = err;

  if (E != NULL)
    E[0] = err;

  for (k=1; k<=order; k++)
  {
    num = CZERO;
    den = RZERO;
    for (i=0; i<n; i++)
    {
      if ( m[i] >= order ) 
      {
        num -= DOTC(m[i]-k, eb[i], 1, &(ef[i][k]), 1);
        den += DOTC(m[i]-k, &(ef[i][k]), 1, &(ef[i][k]), 1) + 
          DOTC(m[i]-k, eb[i], 1, eb[i], 1); 
        
      }
    }
    if (den == RZERO)
    {
      fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
      effective_order = k-1;
      break;
    }
    parcor[k] = (2*num)/den;
    pk = parcor[k];
    cpk = CONJ(pk);

    for (i=0; i<n; i++)
    {
      if ( m[i] >= order ) 
      {
        for (j=0; j<m[i]-k; j++)
        {
          CTYPE tmp = ef[i][j+k];         

          ef[i][j+k] +=  pk * eb[i][j];
          eb[i][j] += cpk * tmp; 
        }
      }
    }

    alpha =  SQR(pk);
    /*    if ( alpha < EPSILON )
          {
          fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
          " %e.\n",k-1,alpha);
          effective_order = k-1;
          break;
          } */
    err_old = err;
    err *= RONE - alpha;

    memcpy(atmp, a, k*sizeof(CTYPE));
    atmp[k] = RZERO;
    for (i=1; i<=k; i++)
      a[i] = atmp[i] + pk * CONJ(atmp[k-i]);

    if (E != NULL)
      E[k] = err;


    if (err < EPSILON*err0 )
    {
      fprintf(stderr,"Model order set to %d since prediction error was"
              " reduced to %f x variance.\n",k-1,EPSILON);
      effective_order = k;
      break;
    } 
  }

  free(parcor);
  free(atmp);
  for (i=0; i<n; i++)
  {
    if ( m[i] >= order ) 
    {
      free(ef[i]);
      free(eb[i]);
    }
  }
  free(ef);
  free(eb);
  return effective_order;
}

#endif /* NOBLAS */

#undef DOTC
#undef MULTI_BURG
#undef EPSILON
#undef SQR
#include "ctypes_undef.h"

Revision:	1.1
Committed:	Sun Apr 28 07:46:58 2013 UTC (10 years, 4 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	globalhs_version_5, Ver_8-7, Ver_8-5, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_8-3, globalhs_version_8, globalhs_version_9, globalhs_version_0, globalhs_version_1, globalhs_version_2, globalhs_version_3, globalhs_version_4, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-2, Ver_8-10, Ver_8-1, Ver_8-6, Ver_9-1, Ver_8-4, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Log Message:	functions needed for detrending and gapfilling
#	Content
1	// $Header: $
2
3	#include "ctypes_def.h"
4
5	// type specific definitions
6	#if TYPE == FLOAT
7	#define MULTI_BURG smulti_burg
8	#define EPSILON (10*FLT_EPSILON)
9	#elif TYPE == DOUBLE
10	#define MULTI_BURG dmulti_burg
11	#define EPSILON (10*DBL_EPSILON)
12	#elif TYPE == COMPLEXFLOAT
13	#define MULTI_BURG cmulti_burg
14	#define EPSILON (10*FLT_EPSILON)
15	#elif TYPE == COMPLEXDOUBLE
16	#define MULTI_BURG zmulti_burg
17	#define EPSILON (10*DBL_EPSILON)
18	#endif
19
20
21	#ifdef NOBLAS
22
23	int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1],
24	RTYPE E[order+1])
25	{
26	int k, i, j, skip, N, n_total, effective_order;
27	CTYPE num, atmp, parcor, ef, eb, pk, cpk;
28	RTYPE den, err, err0, err_old, alpha;
29
30	effective_order = order;
31	skip = 0;
32	N = 0;
33	n_total=0;
34	for (i=0; i<n; i++)
35	{
36	n_total += m[i];
37	if ( m[i] < order )
38	skip++;
39	else
40	N += m[i];
41	}
42	#ifdef DEBUGOUT
43	printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
44	#endif
45
46	if (skip >= n)
47	{
48	fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
49	" Cannot compute AR model.\n",__func__,order);
50	abort();
51	}
52
53	parcor = (CTYPE )malloc((order+1)sizeof(CTYPE));
54	atmp = (CTYPE )malloc((order+1)sizeof(CTYPE));
55	ef = (CTYPE *)malloc(nsizeof(CTYPE*));
56	eb = (CTYPE *)malloc(nsizeof(CTYPE*));
57	for (i=0; i<n; i++)
58	{
59	if ( m[i] >= order )
60	{
61	ef[i] = (CTYPE )malloc(m[i]sizeof(CTYPE));
62	eb[i] = (CTYPE )malloc(m[i]sizeof(CTYPE));
63	memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
64	memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
65	}
66	}
67
68	memset(a,0,(order+1)*sizeof(CTYPE));
69	a[0] = CONE;
70
71	// Initial error.
72	err = 0;
73	for (i=0; i<n; i++)
74	{
75	if ( m[i] >= order )
76	{
77	for (j=0; j<m[i]; j++)
78	err += SQR(x[i][j]);
79	}
80	}
81	err /= N;
82	err0 = err;
83
84	if (E != NULL)
85	E[0] = err;
86
87	for (k=1; k<=order; k++)
88	{
89	num = CZERO;
90	den = RZERO;
91	for (i=0; i<n; i++)
92	{
93	if ( m[i] >= order )
94	{
95	for (j=0; j<m[i]-k; j++)
96	{
97	CTYPE f = ef[i][j+k];
98	CTYPE b = eb[i][j];
99
100	num -= CONJ(b) * f;
101	den += SQR(f) + SQR(b);
102	}
103	}
104	}
105	if (den == RZERO)
106	{
107	fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
108	effective_order = k-1;
109	break;
110	}
111	parcor[k] = (2*num)/den;
112	pk = parcor[k];
113	cpk = CONJ(pk);
114
115	for (i=0; i<n; i++)
116	{
117	if ( m[i] >= order )
118	{
119	for (j=0; j<m[i]-k; j++)
120	{
121	CTYPE tmp = ef[i][j+k];
122
123	ef[i][j+k] += pk * eb[i][j];
124	eb[i][j] += cpk * tmp;
125	}
126	}
127	}
128
129	alpha = SQR(pk);
130	/* if ( alpha < EPSILON )
131	{
132	fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
133	" %e.\n",k-1,alpha);
134	effective_order = k-1;
135	break;
136	} */
137	err_old = err;
138	err *= RONE - alpha;
139
140	memcpy(atmp, a, k*sizeof(CTYPE));
141	atmp[k] = RZERO;
142	for (i=1; i<=k; i++)
143	a[i] = atmp[i] + pk * CONJ(atmp[k-i]);
144
145	if (E != NULL)
146	E[k] = err;
147
148
149	if (err < EPSILON*err0 )
150	{
151	fprintf(stderr,"Model order set to %d since prediction error was"
152	" reduced to %f x variance.\n",k-1,EPSILON);
153	effective_order = k;
154	break;
155	}
156	}
157
158	free(parcor);
159	free(atmp);
160	for (i=0; i<n; i++)
161	{
162	if ( m[i] >= order )
163	{
164	free(ef[i]);
165	free(eb[i]);
166	}
167	}
168	free(ef);
169	free(eb);
170	return effective_order;
171	}
172
173	#else
174
175	int MULTI_BURG( int n, int m[n], CTYPE *x[n], int order, CTYPE a[order+1],
176	RTYPE E[order+1])
177	{
178	int k, i, j, skip, N, n_total, effective_order;
179	CTYPE num, atmp, parcor, ef, eb, pk, cpk;
180	RTYPE den, err, err0, err_old, alpha;
181
182	effective_order = order;
183	skip = 0;
184	N = 0;
185	n_total=0;
186	for (i=0; i<n; i++)
187	{
188	n_total += m[i];
189	if ( m[i] < order )
190	skip++;
191	else
192	N += m[i];
193	}
194	#ifdef DEBUGOUT
195	printf("\nMultiburg using %d out of %d samples.\n",N,n_total);
196	#endif
197
198	if (skip >= n)
199	{
200	fprintf(stderr,"%s: All data segments are shorter than model order (%d)."
201	" Cannot compute AR model.\n",__func__,order);
202	abort();
203	}
204
205	parcor = (CTYPE )malloc((order+1)sizeof(CTYPE));
206	atmp = (CTYPE )malloc((order+1)sizeof(CTYPE));
207	ef = (CTYPE *)malloc(nsizeof(CTYPE*));
208	eb = (CTYPE *)malloc(nsizeof(CTYPE*));
209	for (i=0; i<n; i++)
210	{
211	if ( m[i] >= order )
212	{
213	ef[i] = (CTYPE )malloc(m[i]sizeof(CTYPE));
214	eb[i] = (CTYPE )malloc(m[i]sizeof(CTYPE));
215	memcpy(ef[i],x[i],m[i]*sizeof(CTYPE));
216	memcpy(eb[i],x[i],m[i]*sizeof(CTYPE));
217	}
218	}
219
220	memset(a,0,(order+1)*sizeof(CTYPE));
221	a[0] = CONE;
222
223	// Initial error.
224	err = 0;
225	for (i=0; i<n; i++)
226	{
227	if ( m[i] >= order )
228	err += DOTC(m[i], x[i], 1, x[i],1);
229	}
230	err /= N;
231	err0 = err;
232
233	if (E != NULL)
234	E[0] = err;
235
236	for (k=1; k<=order; k++)
237	{
238	num = CZERO;
239	den = RZERO;
240	for (i=0; i<n; i++)
241	{
242	if ( m[i] >= order )
243	{
244	num -= DOTC(m[i]-k, eb[i], 1, &(ef[i][k]), 1);
245	den += DOTC(m[i]-k, &(ef[i][k]), 1, &(ef[i][k]), 1) +
246	DOTC(m[i]-k, eb[i], 1, eb[i], 1);
247
248	}
249	}
250	if (den == RZERO)
251	{
252	fprintf(stderr,"Data is exactly modelled by an AR(%d) model\n",k-1);
253	effective_order = k-1;
254	break;
255	}
256	parcor[k] = (2*num)/den;
257	pk = parcor[k];
258	cpk = CONJ(pk);
259
260	for (i=0; i<n; i++)
261	{
262	if ( m[i] >= order )
263	{
264	for (j=0; j<m[i]-k; j++)
265	{
266	CTYPE tmp = ef[i][j+k];
267
268	ef[i][j+k] += pk * eb[i][j];
269	eb[i][j] += cpk * tmp;
270	}
271	}
272	}
273
274	alpha = SQR(pk);
275	/* if ( alpha < EPSILON )
276	{
277	fprintf(stderr,"Model order set to %d since error reduction factor - 1 was"
278	" %e.\n",k-1,alpha);
279	effective_order = k-1;
280	break;
281	} */
282	err_old = err;
283	err *= RONE - alpha;
284
285	memcpy(atmp, a, k*sizeof(CTYPE));
286	atmp[k] = RZERO;
287	for (i=1; i<=k; i++)
288	a[i] = atmp[i] + pk * CONJ(atmp[k-i]);
289
290	if (E != NULL)
291	E[k] = err;
292
293
294	if (err < EPSILON*err0 )
295	{
296	fprintf(stderr,"Model order set to %d since prediction error was"
297	" reduced to %f x variance.\n",k-1,EPSILON);
298	effective_order = k;
299	break;
300	}
301	}
302
303	free(parcor);
304	free(atmp);
305	for (i=0; i<n; i++)
306	{
307	if ( m[i] >= order )
308	{
309	free(ef[i]);
310	free(eb[i]);
311	}
312	}
313	free(ef);
314	free(eb);
315	return effective_order;
316	}
317
318	#endif /* NOBLAS */
319
320	#undef DOTC
321	#undef MULTI_BURG
322	#undef EPSILON
323	#undef SQR
324	#include "ctypes_undef.h"
325