1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 /*
8   This file contains routines to initialize and implement the
9   Linear+Ort noise model.
10 
11   File:     model_linplusort.c
12   Author:   RW Cox
13   Date:     24 Jul 2006
14 
15 */
16 
17 
18 /*---------------------------------------------------------------------------*/
19 
20 #include <math.h>
21 #include "NLfit_model.h"
22 
23 
24 void noise_model
25 (
26   float * gn,                /* parameters for noise model */
27   int ts_length,             /* length of time series data */
28   float ** x_array,          /* independent variable matrix */
29   float * ts_array           /* estimated noise model time series */
30 );
31 
32 /*---------------------------------------------------------------------------*/
33 /*
34   Routine to initialize the noise model by defining the number of parameters
35   in the noise model, the name of the noise model, and the default values
36   for the minimum and maximum parameter constraints.
37 */
38 
39 DEFINE_MODEL_PROTOTYPE
40 
initialize_model()41 MODEL_interface * initialize_model ()
42 {
43   MODEL_interface * mi = NULL;
44 
45   /*----- allocate memory space for model interface -----*/
46   mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
47 
48   /*----- define Linear+Ort noise model -----*/
49 
50   /*----- name of this model -----*/
51   strcpy (mi->label, "Linear+Ort");
52 
53   /*----- this is a noise model -----*/
54   mi->model_type = MODEL_NOISE_TYPE;
55 
56   /*----- number of parameters in the model -----*/
57   mi->params = 3;
58 
59   /*----- parameter labels -----*/
60   strcpy (mi->plabel[0], "constant");
61   strcpy (mi->plabel[1], "linear");
62   strcpy (mi->plabel[2], "Ort");
63 
64   /*----- minimum and maximum parameter constraints -----*/
65   mi->min_constr[0] = -100.0;   mi->max_constr[0] = 100.0;
66   mi->min_constr[1] =   -1.0;   mi->max_constr[1] =   1.0;
67   mi->min_constr[2] =   -1.0;   mi->max_constr[2] =   1.0;
68 
69   /*----- function which implements the model -----*/
70   mi->call_func = noise_model;
71 
72 
73   /*----- return pointer to the model interface -----*/
74   return (mi);
75 }
76 
77 
78 /*---------------------------------------------------------------------------*/
79 /*
80   Routine to calculate the time series which results from the linear trend
81   noise model and model parameters.
82 
83   Definition of model parameters:
84 
85      gn[0] = constant coefficient
86      gn[1] = linear coefficient
87      gn[2] = quadratic coefficient
88 */
89 
noise_model(float * gn,int ts_length,float ** x_array,float * ts_array)90 void noise_model
91 (
92   float * gn,                /* parameters for noise model */
93   int ts_length,             /* length of time series data */
94   float ** x_array,          /* independent variable matrix */
95   float * ts_array           /* estimated noise model time series */
96 )
97 
98 {
99   int it;                           /* time index */
100   float t;                          /* time */
101   float fval;                       /* time series value at time t */
102 
103   int ib = ts_length % 4 , nt = ts_length ;
104   float g0=gn[0] , g1=gn[1] , g2=gn[2] ;
105 
106   switch( ib ){
107     case 3: ts_array[2] = g1*x_array[2][1] + g0 + g2*x_array[2][2]; /* fall thru */
108     case 2: ts_array[1] = g1*x_array[1][1] + g0 + g2*x_array[1][2]; /* fall thru */
109     case 1: ts_array[0] = g1*x_array[0][1] + g0 + g2*x_array[0][2]; break ;
110   }
111   for( it=ib ; it < nt ; it+=4 ){
112     ts_array[it  ] = g1*x_array[it  ][1] + g0 + g2*x_array[it  ][2];
113     ts_array[it+1] = g1*x_array[it+1][1] + g0 + g2*x_array[it+1][2];
114     ts_array[it+2] = g1*x_array[it+2][1] + g0 + g2*x_array[it+2][2];
115     ts_array[it+3] = g1*x_array[it+3][1] + g0 + g2*x_array[it+3][2];
116   }
117 }
118