00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 /*#define CLIPSECTIONS */
00008 /* Contrary to the good tradition, this .h file will include */
00009 /* function declarations and definitions. */
00010 /* The reason for this is that the plugin has to be a stand alone*/
00011 /* code and I must provide all the code for the functions. */
00012 /* This is just a way of splitting the code into smaller pieces */
00013 /* This file will hold all support functions for the plugin code */
00014 
00015 /*-------------------------------------------------------------------*/
00016 /* taken from #include "/usr/people/ziad/Programs/C/Z/Zlib/prototype.h" */
00017 
00018 #if defined(SCO) || defined(SOLARIS)
00019 #define drem remainder
00020 #endif
00021 #ifndef NOWAYXCORCOEF
00022         #define NOWAYXCORCOEF 0                                 /* A flag value indicating that something lethal went on */
00023 #endif
00024 
00025 /* do not change these three*/
00026 #define METH_SECONDS 0
00027 #define METH_DEGREES 1
00028 #define METH_RADIANS 2
00029 
00030 
00031 static int Read_file (float *x,char *f_name,int n_points);
00032 
00033 static void linear_reg (float *x,float *y,int size,float *a,float *b,int *err);
00034 
00035 static int equal_strings (char *s1,char *s2);
00036 
00037 static int float_file_size (char *f_name);
00038 
00039 static void error_message (char *origin,char *cause,int ext);
00040 
00041 static char **allocate2D (int rows,int cols,int element_size);
00042 
00043 static void free2D(char **a,int rows);
00044 
00045 static void linear_interp (float *x1,float *y1,float *x2,float *y2,float *x,float *y,int ln);
00046 
00047 static void float_to_complex (float *x,COMPLEX *xc,int ln);
00048 
00049 static void c_conj (COMPLEX *x,COMPLEX *y,int ln);
00050 
00051 static void c_get (COMPLEX *x,float *y,int p,int ln);
00052 
00053 static void c_scale (COMPLEX *x,COMPLEX *y,float scl,int ln);
00054 
00055 static void c_padd (COMPLEX *x,COMPLEX *y,COMPLEX pad_val,int ix,int lnx,int lny);
00056 
00057 static void c_mult (COMPLEX *x,COMPLEX *y,COMPLEX *z,int ln);
00058 
00059 static void detrend (float *y,float *yd,int lny,float *a,float *b);
00060 
00061 static void padd (float *x,float *y,float pad_val,int ix,int lnx,int lny);
00062 
00063 static void hanning (float *x,int l,int opt);
00064 
00065 static float punwrap (float p,int opt );
00066 
00067 static float Lagrange_interp (float *x,float *y,float xi,int ln);
00068 
00069 static void f_mult (float *x,float *y,float *z,int ln); 
00070         
00071 static int hilbertdelay_V2 (float *x,float *y,int lng_full,int Nseg,int Pover,\
00072                        int opt,int dtrnd, float Dtx, int biasrem,\
00073                        float *del,float *slp,float *xcor,float *xcorCoef,\
00074                        float *vx, float *vy);
00075                        
00076 static void hunwrap (float del, float fs, float T, float slp, int wrp, int unt, float *delu );
00077 
00078 static int isarg (int argc, char *argv[], char *probe);
00079 
00080 static float mean_array (float *ar,int size);
00081 
00082 static void zeromean (float *x, float *y, int ln );
00083 
00084 static void disp_comp_vect (COMPLEX *v,int l);
00085 
00086 static void disp_vect (float *v,int l);
00087 
00088 static int is_vect_null ( float * v , int npts );
00089 
00090 static int write_float (float *x,char *f_name,int n_points);
00091 
00092 
00093 /* definition and declaration part to suport the main algorithm */
00094 /* -----------------------END-----------------------------------*/
00095 
00096 
00097 /* -------------------------------------------------------------*/
00098 /* support functions declaration for main algorithm             */
00099 /* -----------------------START---------------------------------*/
00100 
00101 static void error_message (char *s1,char *s2,int ext)
00102  
00103  {
00104  
00105    printf ("\n\n\a\33[1mError: \33[0m%s\n",s2);
00106    printf ("\33[1mError origin:\33[0m %s\n\n",s1);
00107    if (ext == 1)
00108         {
00109                 printf ("Exiting Program ..\n\n");
00110                 exit (0);
00111         }
00112         else return;
00113    
00114   }
00115 
00116 /*-----------------------------------------------------------------------------------*/ 
00117 /**************************************************************************
00118 
00119 allocate2D.c - Make matrix of given size (rows x cols) and type
00120 
00121 The type is given by element_size (2 = ints, 4 = floats, 8 = doubles).
00122 Exits if the matrix could not be allocated.
00123 
00124     char **allocate2D(int rows,int cols,int element_size)
00125 SIZE might vary depending on platform used 
00126 
00127 This function was adapted from DSP_in_C library functions
00128 
00129                                 Ziad Saad       Oct_21_96
00130 *************************************************************************/
00131 
00132 static char **allocate2D (int rows,int cols,int element_size)
00133 
00134 {
00135     int i;
00136     char **A;
00137 
00138 /* try to allocate the request */
00139     switch(element_size) {
00140         case sizeof(short): {    /* integer matrix */
00141             short **int_matrix;
00142             int_matrix = (short **)calloc(rows,sizeof(short *));
00143             if(!int_matrix) {
00144                 printf("\nError making pointers in %dx%d int matrix\n"
00145                             ,rows,cols);
00146                 exit(1);
00147             }
00148             for(i = 0 ; i < rows ; i++) {
00149                 int_matrix[i] = (short *)calloc(cols,sizeof(short));
00150                 if(!int_matrix[i]) {
00151                     printf("\nError making row %d in %dx%d int matrix\n"
00152                             ,i,rows,cols);
00153                     exit(1);
00154                 }
00155             }
00156             A = (char **)int_matrix;
00157             break;
00158         }
00159         case sizeof(float): {    /* float matrix */
00160             float **float_matrix;
00161             float_matrix = (float **)calloc(rows,sizeof(float *));
00162             if(!float_matrix) {
00163                 printf("\nError making pointers in %dx%d float matrix\n"
00164                             ,rows,cols);
00165                 exit(1);
00166             }
00167             for(i = 0 ; i < rows ; i++) {
00168                 float_matrix[i] = (float *)calloc(cols,sizeof(float));
00169                 if(!float_matrix[i]) {
00170                     printf("\nError making row %d in %dx%d float matrix\n"
00171                             ,i,rows,cols);
00172                     exit(1);
00173                 }
00174             }
00175             A = (char **)float_matrix;
00176             break;
00177         }
00178         case sizeof(double): {   /* double matrix */
00179             double **double_matrix;
00180             double_matrix = (double **)calloc(rows,sizeof(double *));
00181             if(!double_matrix) {
00182                 printf("\nError making pointers in %dx%d double matrix\n"
00183                             ,rows,cols);
00184                 exit(1);
00185             }
00186             for(i = 0 ; i < rows ; i++) {
00187                 double_matrix[i] = (double *)calloc(cols,sizeof(double));
00188                 if(!double_matrix[i]) {
00189                     printf("\nError making row %d in %dx%d double matrix\n"
00190                             ,i,rows,cols);
00191                     exit(1);
00192                 }
00193             }
00194             A = (char **)double_matrix;
00195             break;
00196         }
00197         default:
00198             printf("\nERROR in matrix_allocate: unsupported type\n");
00199             exit(1);
00200     }
00201     return(A);
00202 }
00203 
00204 /*-----------------------------------------------------------------------------------*/ 
00205 
00206 /**************************************************************************
00207 
00208 free2D.c - Free all elements of matrix 
00209 
00210 Frees the 2D array (rows and cols) allocated using allocate2D
00211 
00212 Error message and exit if improper structure is
00213 passed to it (null pointers or zero size matrix).
00214 
00215     void free2D(char **a, int rows);
00216 
00217 This function was adapted from DSP_in_C library functions
00218 
00219                                 Ziad Saad                              Oct_22_96
00220 *************************************************************************/
00221 
00222 static void free2D(char **a,int rows)
00223     
00224 {
00225     int i;
00226     
00227 /* free each row of data */
00228     for(i = 0 ; i < rows ; i++) free(a[i]);
00229 
00230 /* free each row pointer */
00231     free((char *)a);
00232     a = NULL;           /* set to null for error */
00233     
00234         return;
00235 }
00236 
00237 /*-----------------------------------------------------------------------------------*/ 
00238 
00239 static void hanning (float *x,int l,int opt)
00240         {
00241                 int i;
00242                 float arg;
00243                 
00244                 arg = 2.0*3.14159/(l-1);
00245                 if (opt == 0)
00246                 {
00247                         for (i=0;i<l;++i)
00248                                 x[i]=0.5-0.5*cos(arg*(float)i);
00249                 } else if (opt == 1)
00250                         {
00251                          for (i=0;i<l;++i)
00252                                 x[i]=x[i] * (0.5-0.5*cos(arg*(float)i));
00253                         }
00254                 return;
00255         } 
00256 
00257 /*-----------------------------------------------------------------------------------*/ 
00258 
00259 static void detrend (float *y,float *yd,int lny,float *a,float *b)
00260 
00261         {
00262          int i;
00263          int err;
00264          float *x;
00265          
00266     x = (float *)calloc (lny+1,sizeof(float));
00267          if (x == NULL)
00268                                 {
00269                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
00270                                         printf ("Abandon Lab Immediately !\n\n");
00271                                         return;
00272                                 };
00273          
00274          for (i=0;i<lny;++i) /*creating x vector */
00275                 x[i] = (float)i;
00276          
00277          linear_reg (x,y,lny,a,b,&err);
00278          
00279          for (i=0;i<lny;++i)
00280                 {
00281                         yd[i] = y[i] - (*a * x[i] + *b);
00282                 }
00283          
00284           free (x);
00285           return;
00286         }
00287 
00288 /*-----------------------------------------------------------------------------------*/ 
00289 
00290 static void padd (float *x,float *y,float pad_val,int ix,int lnx,int lny)
00291         {
00292                 int i,di;
00293                 float *tmp;
00294                 
00295                 tmp = (float *) calloc (lnx+2,sizeof(float));
00296                 
00297                 if (tmp == NULL)
00298                                                 {
00299                                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
00300                                                         printf ("Abandon Lab Immediately !\n\n");
00301                                                         return;
00302                                                 };
00303 
00304                 
00305                 di = lny-lnx;
00306                 if (lny < lnx) 
00307                         {
00308                                 error_message ("padd","lny < lnx !",1);
00309                                 exit(1);
00310                         }
00311                 if (ix > lnx+1)
00312                         {
00313                                 error_message ("padd","ix > lnx+1 !",1);
00314                                 exit(1);
00315                         }
00316                         
00317                 for (i=0;i<lnx;++i)
00318                         {
00319                                 tmp[i] = x[i]; /* must use tmp to be safe when in client program function call is made with input and */
00320                         }
00321                 
00322                 for (i=0;i<(ix-1);++i)
00323                         {
00324                                 y[i] = tmp[i];
00325                         }
00326                 for (i=ix-1;i<ix+di-1;++i)
00327                         {
00328                                 y[i] = pad_val;
00329                         }
00330                 for (i=ix+di-1;i<lny;++i)
00331                         {
00332                                 y[i] = tmp[i-di];
00333                         }
00334                         
00335                 free (tmp);
00336                 return;
00337         }
00338 
00339 /*-----------------------------------------------------------------------------------*/ 
00340 
00341 /**************************************************************************
00342 
00343 fft - In-place radix 2 decimation in time FFT
00344 
00345 Requires pointer to complex array, x and power of 2 size of FFT, m
00346 (size of FFT = 2**m).  Places FFT output on top of input COMPLEX array.
00347 
00348 void fft(COMPLEX *x, int m)
00349 
00350 *************************************************************************/
00351 
00352 static void fft(COMPLEX *x,int m)
00353 {
00354     static COMPLEX *w;           /* used to store the w complex array */
00355     static int mstore = 0;       /* stores m for future reference */
00356     static int n = 1;            /* length of fft stored for future */
00357 
00358     COMPLEX u,temp,tm;
00359     COMPLEX *xi,*xip,*xj,*wptr;
00360 
00361     int i,j,k,l,le,windex;
00362 
00363     double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
00364 
00365     if(m != mstore) {
00366 
00367 /* free previously allocated storage and set new m */
00368 
00369         if(mstore != 0) free(w);
00370         mstore = m;
00371         if(m == 0) return;       /* if m=0 then done */
00372 
00373 /* n = 2**m = fft length */
00374 
00375         n = 1 << m;
00376         le = n/2;
00377 
00378 /* allocate the storage for w */
00379 
00380         w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
00381         if(!w) {
00382             printf("\nUnable to allocate complex W array\n");
00383             exit(1);
00384         }
00385 
00386 /* calculate the w values recursively */
00387 
00388         arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
00389         wrecur_real = w_real = cos(arg);
00390         wrecur_imag = w_imag = -sin(arg);
00391         xj = w;
00392         for (j = 1 ; j < le ; j++) {
00393             xj->real = (float)wrecur_real;
00394             xj->imag = (float)wrecur_imag;
00395             xj++;
00396             wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
00397             wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
00398             wrecur_real = wtemp_real;
00399         }
00400     }
00401 
00402 /* start fft */
00403 
00404     le = n;
00405     windex = 1;
00406     for (l = 0 ; l < m ; l++) {
00407         le = le/2;
00408 
00409 /* first iteration with no multiplies */
00410 
00411         for(i = 0 ; i < n ; i = i + 2*le) {
00412             xi = x + i;
00413             xip = xi + le;
00414             temp.real = xi->real + xip->real;
00415             temp.imag = xi->imag + xip->imag;
00416             xip->real = xi->real - xip->real;
00417             xip->imag = xi->imag - xip->imag;
00418             *xi = temp;
00419         }
00420 
00421 /* remaining iterations use stored w */
00422 
00423         wptr = w + windex - 1;
00424         for (j = 1 ; j < le ; j++) {
00425             u = *wptr;
00426             for (i = j ; i < n ; i = i + 2*le) {
00427                 xi = x + i;
00428                 xip = xi + le;
00429                 temp.real = xi->real + xip->real;
00430                 temp.imag = xi->imag + xip->imag;
00431                 tm.real = xi->real - xip->real;
00432                 tm.imag = xi->imag - xip->imag;             
00433                 xip->real = tm.real*u.real - tm.imag*u.imag;
00434                 xip->imag = tm.real*u.imag + tm.imag*u.real;
00435                 *xi = temp;
00436             }
00437             wptr = wptr + windex;
00438         }
00439         windex = 2*windex;
00440     }            
00441 
00442 /* rearrange data by bit reversing */
00443 
00444     j = 0;
00445     for (i = 1 ; i < (n-1) ; i++) {
00446         k = n/2;
00447         while(k <= j) {
00448             j = j - k;
00449             k = k/2;
00450         }
00451         j = j + k;
00452         if (i < j) {
00453             xi = x + i;
00454             xj = x + j;
00455             temp = *xj;
00456             *xj = *xi;
00457             *xi = temp;
00458         }
00459     }
00460 }
00461 
00462 /*-----------------------------------------------------------------------------------*/ 
00463 
00464 
00465 /**************************************************************************
00466 
00467 ifft - In-place radix 2 decimation in time inverse FFT
00468 
00469 Requires pointer to complex array, x and power of 2 size of FFT, m
00470 (size of FFT = 2**m).  Places inverse FFT output on top of input
00471 frequency domain COMPLEX array.
00472 
00473 void ifft(COMPLEX *x, int m)
00474 
00475 *************************************************************************/
00476 
00477 static void ifft(COMPLEX *x,int m)
00478 {
00479     static COMPLEX *w;           /* used to store the w complex array */
00480     static int mstore = 0;       /* stores m for future reference */
00481     static int n = 1;            /* length of ifft stored for future */
00482 
00483     COMPLEX u,temp,tm;
00484     COMPLEX *xi,*xip,*xj,*wptr;
00485 
00486     int i,j,k,l,le,windex;
00487 
00488     double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
00489     float scale;
00490 
00491     if(m != mstore) {
00492 
00493 /* free previously allocated storage and set new m */
00494 
00495         if(mstore != 0) free(w);
00496         mstore = m;
00497         if(m == 0) return;       /* if m=0 then done */
00498 
00499 /* n = 2**m = inverse fft length */
00500 
00501         n = 1 << m;
00502         le = n/2;
00503 
00504 /* allocate the storage for w */
00505 
00506         w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
00507         if(!w) {
00508             printf("\nUnable to allocate complex W array\n");
00509             exit(1);
00510         }
00511 
00512 /* calculate the w values recursively */
00513 
00514         arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
00515         wrecur_real = w_real = cos(arg);
00516         wrecur_imag = w_imag = sin(arg);  /* opposite sign from fft */
00517         xj = w;
00518         for (j = 1 ; j < le ; j++) {
00519             xj->real = (float)wrecur_real;
00520             xj->imag = (float)wrecur_imag;
00521             xj++;
00522             wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
00523             wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
00524             wrecur_real = wtemp_real;
00525         }
00526     }
00527 
00528 /* start inverse fft */
00529 
00530     le = n;
00531     windex = 1;
00532     for (l = 0 ; l < m ; l++) {
00533         le = le/2;
00534 
00535 /* first iteration with no multiplies */
00536 
00537         for(i = 0 ; i < n ; i = i + 2*le) {
00538             xi = x + i;
00539             xip = xi + le;
00540             temp.real = xi->real + xip->real;
00541             temp.imag = xi->imag + xip->imag;
00542             xip->real = xi->real - xip->real;
00543             xip->imag = xi->imag - xip->imag;
00544             *xi = temp;
00545         }
00546 
00547 /* remaining iterations use stored w */
00548 
00549         wptr = w + windex - 1;
00550         for (j = 1 ; j < le ; j++) {
00551             u = *wptr;
00552             for (i = j ; i < n ; i = i + 2*le) {
00553                 xi = x + i;
00554                 xip = xi + le;
00555                 temp.real = xi->real + xip->real;
00556                 temp.imag = xi->imag + xip->imag;
00557                 tm.real = xi->real - xip->real;
00558                 tm.imag = xi->imag - xip->imag;             
00559                 xip->real = tm.real*u.real - tm.imag*u.imag;
00560                 xip->imag = tm.real*u.imag + tm.imag*u.real;
00561                 *xi = temp;
00562             }
00563             wptr = wptr + windex;
00564         }
00565         windex = 2*windex;
00566     }            
00567 
00568 /* rearrange data by bit reversing */
00569 
00570     j = 0;
00571     for (i = 1 ; i < (n-1) ; i++) {
00572         k = n/2;
00573         while(k <= j) {
00574             j = j - k;
00575             k = k/2;
00576         }
00577         j = j + k;
00578         if (i < j) {
00579             xi = x + i;
00580             xj = x + j;
00581             temp = *xj;
00582             *xj = *xi;
00583             *xi = temp;
00584         }
00585     }
00586 
00587 /* scale all results by 1/n */
00588     scale = (float)(1.0/n);
00589     for(i = 0 ; i < n ; i++) {
00590         x->real = scale*x->real;
00591         x->imag = scale*x->imag;
00592         x++;
00593     }
00594 }
00595 
00596 /************************************************************
00597 
00598 rfft - trig recombination real input FFT
00599 
00600 Requires real array pointed to by x, pointer to complex
00601 output array, y and the size of real FFT in power of
00602 2 notation, m (size of input array and FFT, N = 2**m).
00603 On completion, the COMPLEX array pointed to by y 
00604 contains the lower N/2 + 1 elements of the spectrum.
00605 
00606 void rfft(float *x, COMPLEX *y, int m)
00607 
00608 ***************************************************************/
00609 
00610 static void rfft(float *x,COMPLEX *y,int m)
00611 {
00612     static    COMPLEX  *cf;
00613     static    int      mstore = 0;
00614     int       p,num,k,index;
00615     float     Realsum, Realdif, Imagsum, Imagdif;
00616     double    factor, arg;
00617     COMPLEX   *ck, *xk, *xnk, *cx;
00618 
00619 /* First call the fft routine using the x array but with
00620    half the size of the real fft */
00621 
00622     p = m - 1;
00623     cx = (COMPLEX *) x;
00624     fft(cx,p);
00625 
00626 /* Next create the coefficients for recombination, if required */
00627 
00628     num = 1 << p;    /* num is half the real sequence length.  */
00629 
00630     if (m!=mstore){
00631       if (mstore != 0) free(cf);
00632       cf = (COMPLEX *) calloc(num - 1,sizeof(COMPLEX));
00633       if(!cf){
00634         printf("\nUnable to allocate trig recomb coefficients.");
00635         exit(1);
00636       }
00637 
00638       factor = 4.0*atan(1.0)/num;
00639       for (k = 1; k < num; k++){
00640         arg = factor*k;
00641         cf[k-1].real = (float)cos(arg);
00642         cf[k-1].imag = (float)sin(arg);
00643       }
00644     }  
00645 
00646 /* DC component, no multiplies */
00647     y[0].real = cx[0].real + cx[0].imag;
00648     y[0].imag = 0.0;
00649 
00650 /* other frequencies by trig recombination */
00651     ck = cf;
00652     xk = cx + 1;
00653     xnk = cx + num - 1;
00654     for (k = 1; k < num; k++){
00655       Realsum = ( xk->real + xnk->real ) / 2;
00656       Imagsum = ( xk->imag + xnk->imag ) / 2;
00657       Realdif = ( xk->real - xnk->real ) / 2;
00658       Imagdif = ( xk->imag - xnk->imag ) / 2;
00659 
00660       y[k].real = Realsum + ck->real * Imagsum
00661                           - ck->imag * Realdif ;
00662 
00663       y[k].imag = Imagdif - ck->imag * Imagsum
00664                           - ck->real * Realdif ;
00665       ck++;
00666       xk++;
00667       xnk--;
00668     }
00669 }
00670 
00671 /*-----------------------------------------------------------------------------------*/ 
00672 
00673 static void float_to_complex (float *x,COMPLEX *xc,int ln)
00674 
00675         {
00676          int i;
00677       if ((ln-1) == 0)
00678          {       
00679               (*xc).real = (*x);
00680               (*xc).imag = 0.0;
00681               return;
00682            }
00683          else 
00684               {
00685                    for (i=0;i<ln;++i)
00686                        {
00687                           xc[i].real = x[i];
00688                           xc[i].imag = 0.0;
00689                         }
00690                      return;
00691                 }
00692          }
00693 
00694 /*-----------------------------------------------------------------------------------*/ 
00695 
00696 static void c_conj (COMPLEX *x,COMPLEX *y,int ln)
00697         {
00698                 int i;
00699                 if ((ln-1) == 0)
00700                         {       
00701                                 (*y).real = (*x).real;
00702                                 (*y).imag = -(*x).imag;
00703                                 return;
00704                         }
00705                 else 
00706                         {
00707                                 for (i=0;i<ln;++i)
00708                                         {
00709                                                 y[i].real = x[i].real;
00710                                                 y[i].imag = -x[i].imag;
00711                                         }
00712                                 return;
00713                         }
00714         }
00715 
00716 /*-----------------------------------------------------------------------------------*/ 
00717 
00718 static void c_mult (COMPLEX *x,COMPLEX *y,COMPLEX *z,int ln)
00719         {
00720                 int i;
00721                 COMPLEX t;
00722                 
00723                 if ((ln-1) == 0)
00724                         {       
00725                                 t.real = ((*x).real * (*y).real - (*x).imag * (*y).imag);
00726                                 t.imag = ((*x).real * (*y).imag + (*y).real * (*x).imag);
00727                                 (*z).real = t.real;
00728                                 (*z).imag = t.imag;
00729                                 return;
00730                         }
00731                 else 
00732                         {
00733                                 for (i=0;i<ln;++i)
00734                                         {
00735                                                 t.real = (x[i].real * y[i].real - x[i].imag * y[i].imag);
00736                                                 t.imag = (x[i].real * y[i].imag + y[i].real * x[i].imag);
00737                                                 z[i].real = t.real;
00738                                                 z[i].imag = t.imag;
00739                                         }
00740                                 return;
00741                         }
00742         }
00743 
00744 /*-----------------------------------------------------------------------------------*/ 
00745 
00746 static void c_padd (COMPLEX *x,COMPLEX *y,COMPLEX pad_val,int ix,int lnx,int lny)
00747         {
00748                 int i,di;
00749                 COMPLEX *tmp;
00750                 
00751                 tmp = (COMPLEX *) calloc (lnx+2,sizeof(COMPLEX));
00752                 
00753                 if (tmp == NULL)
00754                                                 {
00755                                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
00756                                                         printf ("Abandon Lab Immediately !\n\n");
00757                                                         return;
00758                                                 };
00759 
00760                 di = lny-lnx;
00761                 if (lny < lnx) 
00762                         {
00763                                 error_message ("c_padd","lny < lnx !",1);
00764                                 exit(1);
00765                         }
00766                 if (ix > lnx+1)
00767                         {
00768                                 error_message ("c_padd","ix > lnx+1 !",1);
00769                                 exit(1);
00770                         }
00771                         
00772                 for (i=0;i<lnx;++i)
00773                         {
00774                                 tmp[i].real = x[i].real; /* must do this to be safe when in client program function call is made with input and */
00775                                 tmp[i].imag = x[i].imag; /* output vectors being the same !!! */
00776                         }
00777                 
00778                 for (i=0;i<(ix-1);++i)
00779                         {
00780                                 y[i].real = tmp[i].real;
00781                                 y[i].imag = tmp[i].imag;
00782                         }
00783                 for (i=ix-1;i<ix+di-1;++i)
00784                         {
00785                                 y[i].real = pad_val.real;
00786                                 y[i].imag = pad_val.imag;
00787                         }
00788                 for (i=ix+di-1;i<lny;++i)
00789                         {
00790                                 y[i].real = tmp[i-di].real;
00791                                 y[i].imag = tmp[i-di].imag;
00792                         }
00793                         
00794                 free (tmp);
00795                 return;
00796         }
00797 
00798 /*-----------------------------------------------------------------------------------*/ 
00799 
00800 static void c_scale (COMPLEX *x,COMPLEX *y,float scl,int ln)   
00801         {
00802                 int i;
00803                 if ((ln-1) == 0)
00804                         {       
00805                                 (*y).real = (scl) * ((*x).real);
00806                                 (*y).imag = (scl) * ((*x).imag);
00807                                 return;
00808                         }
00809                 else 
00810                         {
00811                                 for (i=0;i<ln;++i)
00812                                         {
00813                                                 y[i].real = (scl) * (x[i].real);
00814                                                 y[i].imag = (scl) * (x[i].imag);
00815                                         }
00816                                 return;
00817                         }
00818         }
00819 
00820 /*-----------------------------------------------------------------------------------*/ 
00821 
00822 static void c_get (COMPLEX *x,float *y,int p,int ln)
00823 
00824         {
00825                         int i;
00826                   if ((ln-1) == 0)
00827                         {
00828                                 if (p == 0) { (*y) = (*x).real; }
00829                                         else {(*y) = (*x).imag;}
00830                                 return;
00831                         }
00832                 else 
00833                         {
00834                                 for (i=0;i<ln;++i)
00835                                         {
00836                                           if (p == 0) { y[i] = x[i].real; }
00837                                                 else { y[i] = x[i].imag; }
00838                                         }
00839                                 return;
00840                         }
00841         }
00842 
00843 /*-----------------------------------------------------------------------------------*/ 
00844 
00845 static void linear_interp (float *x1,float *y1,float *x2,float *y2,float *x,float *y,int ln)
00846 
00847    {int i;
00848         
00849     if ((ln -1) == 0)       
00850       {
00851         if ((*x2 - *x1) == 0.0) 
00852           {
00853               error_message ("linear_interp","identical X values in interpolation boundaries, causes division by zero !",1);
00854               exit (1);
00855            }
00856         *y = (*x - *x1)/(*x2 - *x1)*(*y2 - *y1) + *y1;
00857        }
00858          else 
00859         {
00860           for (i=0;i<ln;++i)
00861             {
00862               if ((x2[i] - x1[i]) == 0.0) 
00863                 {
00864                    error_message ("linear_interp","identical X values in interpolation boundaries, causes division by zero !",1);
00865                    exit (1);
00866                  }
00867            y[i] = (x[i] - x1[i])/(x2[i] - x1[i])*(y2[i] - y1[i]) + y1[i];
00868              }
00869          }
00870                 return;
00871      }
00872 
00873 /*-----------------------------------------------------------------------------------*/ 
00874 
00875 static int equal_strings (char *s1,char *s2)
00876 
00877  {
00878    int i=0;
00879    
00880    if (s1 == NULL && s2 == NULL) return (-2);
00881    
00882    if ((s1 == NULL && s2 != NULL) || (s1 != NULL && s2 == NULL)) return (-1);
00883    
00884    while (s1[i] == s2[i] 
00885                         && s1[i] != '\0' && s2[i] != '\0') ++i;
00886                         
00887         if (s1[i] == '\0' && s2[i] == '\0') return (1);
00888          else return (0);
00889  
00890  }
00891 
00892 /*-----------------------------------------------------------------------------------*/ 
00893 
00894 static int float_file_size (char *f_name)
00895    
00896     { 
00897       
00898 
00899      int cnt=0,ex;
00900      float buf;
00901      
00902      FILE*internal_file;
00903      
00904      internal_file = fopen (f_name,"r");
00905      if (internal_file == NULL) {
00906              printf ("\aCould not open %s \n",f_name);
00907              printf ("Exiting @ float_file_size function\n");
00908              exit (0);
00909                                  }
00910      ex = fscanf (internal_file,"%f",&buf);                                             
00911      while (ex != EOF)
00912       {
00913         ++cnt;
00914         ex = fscanf (internal_file,"%f",&buf);
00915       }
00916       
00917       
00918       fclose (internal_file);
00919       return (cnt);                                                          
00920    }
00921 
00922 /*-----------------------------------------------------------------------------------*/ 
00923 
00924 static int Read_file (float *x,char *f_name,int n_points)
00925    
00926     { /* pass a 0 to n_points if you want to read till EOF */
00927      int cnt=0,ex,dec;
00928      
00929      FILE*internal_file;
00930      
00931      internal_file = fopen (f_name,"r");
00932      if (internal_file == NULL) {
00933                                                                 printf ("\aCould not open %s \n",f_name);
00934                                                                 printf ("Exiting @ Read_file function\n");
00935                                                                 exit (0);
00936                                                         }
00937      ex = fscanf (internal_file,"%f",&x[cnt]);                                          
00938      while (ex != EOF)
00939       {
00940         ++cnt;
00941         /* NOT WORKING, RETURNS SIZEOF (FLOAT) .....
00942         if (sizeof(x) < cnt)
00943                 {
00944                   printf ("%d = sizeof(x)\n",sizeof(x));
00945                   printf ("\nNot Enough Memory Allocated \n\a");
00946                   printf ("Exiting @Read_file function\n");
00947                   exit (0);
00948                 }
00949         ............................................ */
00950         ex = fscanf (internal_file,"%f",&x[cnt]);
00951         
00952         if ((n_points != 0) && (cnt == n_points)) ex = EOF;
00953       }
00954       
00955       if (cnt < n_points) 
00956         {
00957          printf ("\a\nAttempt to read %d points failed,\n",n_points);
00958          printf (" file contains %d points only.\n",cnt);
00959          do {
00960          
00961          printf ("End Execution (Yes (1) No (0) ? : ");
00962          ex=scanf ("%d",&dec);
00963          } while (ex != 1 || (dec != 1 && dec !=0));
00964          if (dec)
00965           {
00966             printf ("Exiting @ Read_file function\n");
00967                 exit (0);
00968                 }
00969          else printf ("\nContinuing execution with %d points\n",cnt);
00970 
00971         }
00972       
00973       fclose (internal_file);
00974       return (cnt);                                                          
00975    }
00976 
00977 /*-----------------------------------------------------------------------------------*/ 
00978 
00979 static void linear_reg (float *x,float *y,int size,float *a,float *b,int *err)
00980 
00981         {/* linear_reg*/
00982           int i;
00983           float n,sum_x,sum_y,sum_x2,sum_xy;
00984           
00985           *err = 3;
00986           
00987           if (size <= 0) {
00988                                         *err = 1;
00989                                         return;
00990                                         }
00991                                 
00992        sum_x = 0.0;
00993        sum_y = 0.0;
00994        sum_xy = 0.0;       
00995            sum_x2 = 0.0;        
00996            n = (float)size;   
00997 
00998            for (i=0;i<size;++i)
00999              {
01000                sum_x = sum_x + x[i];
01001                sum_y = sum_y + y[i];
01002                sum_x2 = sum_x2 + x[i] * x[i];
01003                sum_xy = sum_xy + x[i] * y[i];
01004              }
01005            
01006            *a = ( n * sum_xy - sum_x * sum_y) / ( n * sum_x2 - sum_x * sum_x);
01007            
01008            *b = (sum_y - *a * sum_x) / n;
01009            
01010             *err = 0; 
01011              return ;
01012         }/* linear_reg */
01013 
01014 /*-----------------------------------------------------------------------------------*/ 
01015 
01016 static float punwrap (float p,int opt )
01017 {/*punwrap*/
01018  float topi,alf;
01019  
01020  if (opt == 0)
01021                 topi = 2.*3.14159;
01022         else if (opt == 1) topi = 360.0;
01023                 else {
01024                         error_message ("punwrap","wrong opt parameter",1);
01025                         exit (1);
01026                 }
01027  alf = (float) drem ((double)p,(double)topi);
01028  
01029  if (alf > topi/2.0) alf = topi/2.0-alf;
01030  
01031  return (alf);
01032 }/*punwrap*/
01033 
01034 /*-----------------------------------------------------------------------------------*/ 
01035 
01036 static float Lagrange_interp (float *x,float *y,float xi,int ln)
01037  
01038         {
01039                 int i,j;
01040                 float yi,p;
01041                 
01042                 yi = 0.0;
01043                 
01044                 for (i=0;i<ln;++i)
01045                 {
01046                         p = y[i];
01047                         for (j=0;j<ln;++j)
01048                                 {
01049                                         if (j != i) 
01050                                                 {
01051                                                         p = p * (xi - x[j]) / (x[i] - x[j]);
01052                                                 }
01053                                 }
01054                         yi=yi+p;
01055                 }       
01056                 return (yi);
01057         }
01058  
01059 /*-----------------------------------------------------------------------------------*/ 
01060 static int isarg (int argc, char *argv[], char *probe)
01061 {/*isarg*/
01062         int i=1;
01063                 
01064         while (i < argc)
01065                 {
01066                         if (equal_strings (argv[i],probe) == 1)
01067                                 {
01068                                         return (i);
01069                                 }
01070                         ++i;
01071                 }
01072         
01073         return (0);
01074         
01075 }/*isarg*/
01076 
01077 /*-----------------------------------------------------------------------------------*/ 
01078 static float mean_array (float *ar,int size)
01079          
01080         {/*Mean_array ()*/
01081          
01082          int i;
01083          float sum,mean;
01084          
01085          sum=0;
01086          
01087          for (i=0;i<size;++i)
01088           {
01089                 sum=sum+ar[i];
01090           }
01091           
01092           mean=sum/(float)size;
01093           
01094           
01095           return (mean);
01096         
01097         }/*Mean_array ()*/
01098 
01099 /*-----------------------------------------------------------------------------------*/ 
01100 static void zeromean (float *x, float *y, int ln )
01101 {/*zeromean*/
01102         int i;
01103         float meanx;
01104         
01105         meanx = mean_array (x,ln);
01106         
01107         for (i=0;i<ln;++i)
01108                 {
01109                         y[i] = x[i] - meanx;
01110                 }
01111         return;
01112 }/*zeromean*/
01113 
01114 /*-----------------------------------------------------------------------------------*/ 
01115 static void f_mult (float *x,float *y,float *z,int ln)  
01116         {
01117                 int i;
01118                 if ((ln-1) == 0)
01119                         {       
01120                                 (*z) = (*y) * (*x);
01121                                 return;
01122                         }
01123                 else 
01124                         {
01125                                 for (i=0;i<ln;++i)
01126                                         {
01127                                                 z[i] = y[i] * x[i];
01128                                         }
01129                                 return;
01130                         }
01131         }
01132 
01133 /*-----------------------------------------------------------------------------------*/ 
01134 /* Important :
01135         Before you replace this function by a new version, make note of the following
01136         changes to the function in here: maxdel is changed from lng /3 to lng /2, 
01137         and the funcion returns a 1, 2 or 3 in case of encoutered errors instead of 
01138         a regular 1 
01139         Also the function does not output any warning messages to the screen if the delay is
01140         larger than one 1/2 a segment length*/
01141 
01142 /* The difference between       hilbertdelay_V2 and hilbertdelay is that the parameter negslp is not used anymore in V2 version */
01143 static int hilbertdelay_V2 (float *x,float *y,int lng_full,int Nseg,int Pover,int opt,int dtrnd, float Dtx, int biasrem, float *del,float *slp,float *xcor,float *xcorCoef, float *vx, float *vy)
01144         {       
01145          static int i_call=0,olng,m,lng_use,lng,strt,nd,sg,cnt,maxdel = 0;
01146          static COMPLEX *fftx,*ffty,*Pxy,*fftyc,*fftxc,*Pxx,*Pyy,*Rxx,*Ryy,**fftyca,**Pxya,**Rxxa,**Ryya;
01147          static float *Px,*Py,*Rxy,*HRxy,*xp,*yp,*xcpy,*ycpy,*tmp_f_vect,*tmp_f_vect2,*ubias,a,\
01148                                 var_y,varu_y,stdv_y,stdvu_y,var_x,varu_x,stdv_x,stdvu_x;
01149          
01150                           int i,j,k;
01151                           char buf[100];
01152                           float *mPxy,tmp,sPx,sPy,alfx,betx,alfy,bety,\
01153                                                 f_i,f_i1,izero,reg_pnt,\
01154                                                 NoWayDelay = -100000.0,\
01155                                                 NoWayxcorCoef = NOWAYXCORCOEF;
01156                                                 
01157                           COMPLEX tmp_c;
01158                           char cbuf[30];
01159 
01160 if ((opt == 0) && (i_call == 0)) 
01161         {
01162                 error_message ("hilbertdelay_V2","Nothing to DO !",0);
01163                 return (1); 
01164         }
01165 
01166         
01167 *del = NoWayDelay;                                      /* setting the value of delay to an unlikely value ...*/
01168 *xcorCoef = NoWayxcorCoef;                      /* setting the cross correlation coefficient to an unlikely value ...*/
01169 
01170 
01171 if (opt > 0)                                                    /* Execution mode */
01172 {/* opt >0 */
01173 
01174 #ifdef DEBUG_ZEO_2              
01175         printf ("\nFunction call #%d\n",i_call);
01176 #endif  
01177 
01178 
01179 /*-----------------------------------------------------------------------------------*/         
01180 /* Steps that need to be perfromed the first time the function is called                                 */
01181 /*-----------------------------------------------------------------------------------*/ 
01182         
01183         if (i_call == 0)
01184                 {/*i_call == 0*/
01185                         if ((Nseg < 1) || (Nseg >= lng_full/5))
01186                                 {
01187                                         sprintf (buf,"Number of segments (%d) null or too large, or vector length too small (%d) for %d segments ",Nseg,lng_full,Nseg);
01188                                         error_message ("hilbertdelay_V2",buf,0);
01189                                         return (2); 
01190                                 } 
01191                         
01192                         lng = (int)((float)lng_full / (float)Nseg);     /* calculating individual segment length */
01193                         
01194                         maxdel = lng/2;                         /* delays should not exceed one third of the segment length (and that's pushing it !)*/
01195                         
01196                         m=0;
01197                         while ((int)pow((double)2,(double)m) < lng) ++m;        /* find closest power of two to the length of the series */
01198                         olng = lng;                                                                             /* save old length*/
01199                         lng = (int)pow((double)2,(double)m);    /* set new length as power of two actual padding length will double*/
01200                                                                                                                                 /* in order to correct for circular convolution effects */
01201                         lng_use = lng;                                                  /* useful length of spectrum after correction for circular convolution effects */
01202 
01203 #ifdef DEBUG_ZEO_2
01204         printf ("selected m=%d for a padded segment length of %d, old segment length was %d\nVector holds %d segments\n",m,lng,olng,Nseg);      
01205 #endif
01206                 
01207                         fftx = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01208                         fftxc = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01209                         ffty = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01210                         fftyc = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01211                         Pxy = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01212                         Px = (float *) calloc ((2*lng)+2,sizeof(float));
01213                         Pxx = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01214                         Py = (float *) calloc ((2*lng)+2,sizeof(float));
01215                         Pyy = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01216                         Rxx = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01217                         Ryy = (COMPLEX *) calloc ((2*lng)+2,sizeof(COMPLEX));
01218                         Rxy = (float *) calloc ((2*lng)+2,sizeof(float));
01219                         HRxy = (float *) calloc ((2*lng)+2,sizeof(float));
01220                         xcpy = (float *) calloc ((2*olng)+2,sizeof(float));
01221                         ycpy = (float *) calloc ((2*olng)+2,sizeof(float));
01222                         xp = (float *) calloc ((2*lng)+2,sizeof(float));
01223                         yp = (float *) calloc ((2*lng)+2,sizeof(float));
01224                         ubias = (float *) calloc ((2*lng)+2,sizeof(float));
01225                         tmp_f_vect = (float *) calloc ((2*lng)+2,sizeof(float));
01226                         tmp_f_vect2 = (float *) calloc ((2*lng)+2,sizeof(float));
01227                         fftyca = (COMPLEX **) allocate2D ((2*lng)+2,Nseg,sizeof(COMPLEX));
01228                         Pxya = (COMPLEX **) allocate2D ((2*lng)+2,Nseg,sizeof(COMPLEX));
01229                         Ryya = (COMPLEX **) allocate2D ((2*lng)+2,Nseg,sizeof(COMPLEX));
01230                         Rxxa = (COMPLEX **) allocate2D ((2*lng)+2,Nseg,sizeof(COMPLEX));        
01231                         
01232                         if (fftx == NULL ||  fftxc == NULL ||  ffty == NULL ||  fftyc == NULL ||  \
01233                             Pxy == NULL ||  Px == NULL ||  Py == NULL ||  xp == NULL ||  yp == NULL ||  \
01234                                  ubias == NULL ||  tmp_f_vect == NULL ||  Pxx == NULL ||  Pyy == NULL ||  \
01235                                  Rxx == NULL ||  Ryy == NULL ||  fftyca == NULL ||  Pxya == NULL ||  \
01236                                  Ryya == NULL ||  Rxxa == NULL ||  tmp_f_vect2 == NULL)
01237                                 {
01238                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
01239                                         printf ("Abandon Lab Immediately !\n\n");
01240                                         return(2);
01241                                 };
01242 
01243                         /* creating a vector to remove the bowtie artifact from the auto and cross correlation curves, and set to zero 
01244                                                         their irrelevant values */
01245                         if (biasrem == 1)
01246                         {
01247                                 for (i=0;i<(2*lng);++i)
01248                                         {               
01249                                                 if (i < olng)
01250                                                         {
01251                                                                 ubias[i] = (float)olng / (float)(olng - i);
01252                                                         }
01253                                                 else
01254                                                         {
01255                                                                 ubias[i] = 0.0;
01256                                                         }
01257                                         }
01258                         }
01259                         
01260                         a = (float)olng; /* Scaling parameter for spectral estimates  */
01261                         
01262                         strt = 0;                               /* setting up for segments loop */
01263                         nd = 0;
01264                         for (sg=0;sg<Nseg;++sg)                                 /* segments counter */
01265                         {/* for sg */
01266                                 strt = sg * olng;
01267                                 nd = strt + olng;
01268                                 
01269                                 cnt = 0;
01270                                 for (i=strt;i<nd;++i) 
01271                                         {
01272                                                 tmp_f_vect[cnt] = y[i];         /* copying segment number (sg+1) */
01273                                                 ++cnt;
01274                                         }
01275                                         
01276                                 if (dtrnd == 1)
01277                                         {/*dtrnd == 1*/
01278                                                 detrend (tmp_f_vect,ycpy,olng,&alfy,&bety);             /* removing linear trend alf? and bet? are the linear fit coefficients for the reference time course*/
01279                                                 padd (ycpy,yp,0.0,olng+1,olng,2*lng);                                   /* Zero padding of detrended reference time course*/
01280                                         }/*dtrnd == 1*/
01281                                 else
01282                                         {/*dtrnd != 1*/
01283                                                 zeromean (tmp_f_vect,ycpy,olng);                                        /* removing the mean only */
01284                                                 padd (ycpy,yp,0.0,olng+1,olng,2*lng);                   /* Zero padding of reference time course*/
01285                                         }/*dtrnd != 1*/
01286                                                         
01287                                 float_to_complex (yp,ffty,2*lng);                       /* reformatting reference vector into complex vector */
01288                 
01289                                 fft (ffty,m+1);                                                                 /* Radix 2 Butterfly fft computation of reference vector */
01290                                 
01291                                 c_conj (ffty,fftyc,2*lng);                                                      /* Computing the conjugate of ffty */
01292                                 
01293                                 c_mult (ffty,fftyc,Ryy,2*lng);          /* Powerspectrum of y (called Ryy) */
01294                                 
01295                                 for (i=0;i<2*lng;++i)                                           /* copying Power spectrum of y and conjugate of fft of individual segment into storage matrix */
01296                                         {
01297                                                 fftyca[i][sg].real = fftyc[i].real;
01298                                                 fftyca[i][sg].imag = fftyc[i].imag;
01299                                                 Ryya[i][sg].real = Ryy[i].real;
01300                                                 Ryya[i][sg].imag = Ryy[i].imag;
01301                                         }
01302                                 
01303                         }/* for sg */
01304                         
01305                         for (sg=0;sg<Nseg;++sg)                                 /* computing sum periodogram of y */
01306                                 {/* for sg */
01307                                         for (i=0;i<2*lng;++i)
01308                                                 {/* for i */
01309                                                         
01310                                                         if (sg == 0) 
01311                                                                 {
01312                                                                         Ryy[i].real = Ryya[i][sg].real;
01313                                                                         Ryy[i].imag = Ryya[i][sg].imag;
01314                                                                 }
01315                                                         else 
01316                                                                 {/* sg > 0 */
01317                                                                         Ryy[i].real = Ryy[i].real + Ryya[i][sg].real;
01318                                                                         Ryy[i].imag = Ryy[i].imag + Ryya[i][sg].imag;
01319                                                                 }/* sg > 0 */
01320                                                         
01321                                                 }/* for i */
01322                                                 
01323                                 }/* for sg */
01324                         
01325                         c_scale (Ryy,Ryy,1.0/((float)Nseg * a),2*lng);                                  /* scaling periodogram */
01326                         
01327                         ifft (Ryy,m+1);                                                                                                 /* calculating autocorrelation of y*/
01328                 
01329                 }/*i_call == 0*/
01330 
01331 /*-----------------------------------------------------------------------------------*/         
01332 /* Steps that need to be repeated each time the function is called with a new vector */
01333 /*-----------------------------------------------------------------------------------*/ 
01334 
01335                 strt = 0;                                                       /* setting up for segments loop */
01336                 nd =0;
01337                 for (sg=0;sg<Nseg;++sg)                                 /* segments counter */
01338                         {/* for sg */
01339                                 strt = sg * olng;
01340                                 nd = strt + olng;
01341                                 
01342                                 cnt = 0;
01343                                 for (i=strt;i<nd;++i) 
01344                                         {
01345                                                 tmp_f_vect[cnt] = x[i];         /* copying segment number (sg+1) */
01346                                                 ++cnt;
01347                                         }
01348 
01349                                 if (dtrnd == 1)
01350                                         {/*dtrnd == 1*/
01351                                                 detrend (tmp_f_vect,xcpy,olng,&alfx,&betx);             /* removing linear trend alf? and bet? are the linear fit coefficients*/
01352                                                 padd (xcpy,xp,0.0,olng+1,olng,2*lng);                                   /* Zero padding of detrended time course*/
01353                                         }/*dtrnd == 1*/
01354                                 else
01355                                         {/*dtrnd != 1*/
01356                                            zeromean (tmp_f_vect,xcpy,olng);                                                     /* removing the mean only*/
01357                                                 padd (xcpy,xp,0.0,olng+1,olng,2*lng);                   /* Zero padding of time course*/
01358                                         }/*dtrnd != 1*/
01359                                 
01360                                 float_to_complex (xp,fftx,2*lng);                                       /* reformatting vector into complex vector */
01361                 
01362                                 fft (fftx,m+1);                                                                                 /* Radix 2 Butterfly fft computation of vector */
01363                                 
01364                                 c_conj (fftx,fftxc,2*lng);                                                      /* Computing the conjugate of fftx */
01365                                 
01366                                 c_mult (fftx,fftxc,Rxx,2*lng);                                          /* Powerspectrum of x (called Rxx) */
01367                                 
01368 #ifdef DEBUG_ZEO_2      
01369                 write_float (xp,"dbg_xdp.txt",2*lng);
01370                 write_float (yp,"dbg_ydp.txt",2*lng);
01371                 printf ("a = %f\n",a);
01372 #endif
01373                                 
01374                                 for (i=0;i<2*lng;++i)                                                                   /* copying fftyc at segment sg from storage array */
01375                                         {
01376                                                 fftyc[i].real = fftyca[i][sg].real;
01377                                                 fftyc[i].imag = fftyca[i][sg].imag;
01378                                         }
01379                                 
01380                                 c_mult (fftx,fftyc,Pxy,2*lng);                                  /* Computing the cross power spectrum */
01381                                 
01382                                 for (i=0;i<2*lng;++i)                                                                   /* storing the power spectrum and the cross power spectrum at */
01383                                         {                                                                                                               /*different segments */
01384                                                 Pxya[i][sg].real = Pxy[i].real;
01385                                                 Pxya[i][sg].imag = Pxy[i].imag;
01386                                                 Rxxa[i][sg].real = Rxx[i].real;
01387                                                 Rxxa[i][sg].imag = Rxx[i].imag;
01388                                         }
01389 
01390                         }/* for sg */
01391                 
01392                 for (sg=0;sg<Nseg;++sg)                                                         /* calculating the sum of the periodograms */
01393                         {/* for sg */
01394                                 for (i=0;i<2*lng;++i)
01395                                         {/* for i*/
01396                                                 if (sg == 0)
01397                                                         {       
01398                                                                 Pxy[i].real = Pxya[i][sg].real;
01399                                                                 Pxy[i].imag = Pxya[i][sg].imag;
01400                                                                 Rxx[i].real = Rxxa[i][sg].real;
01401                                                                 Rxx[i].imag = Rxxa[i][sg].imag;
01402                                                         }
01403                                                 else 
01404                                                         {       
01405                                                                 Pxy[i].real = Pxy[i].real + Pxya[i][sg].real;
01406                                                                 Pxy[i].imag = Pxy[i].imag + Pxya[i][sg].imag;
01407                                                                 Rxx[i].real = Rxx[i].real + Rxxa[i][sg].real;
01408                                                                 Rxx[i].imag = Rxx[i].imag + Rxxa[i][sg].imag;
01409                                                         }
01410                                         }/* for i */
01411                         }/* for sg */
01412                 
01413                 
01414                 c_scale (Rxx,Rxx,1.0/((float)Nseg * a),2*lng); /* calculating average Rxx periodogram   */
01415                 
01416                 c_scale (Pxy,Pxy,2.0/((float)Nseg * a),2*lng);  /* calculating average Pxy and scaling it by 2 periodogram      */
01417                         
01418                 tmp_c.real = 0.0;                                                                       /*discarding half of the cross power spectrum and padding it back to lng */
01419                 tmp_c.imag = 0.0;
01420                 
01421                 c_padd (Pxy,Pxy,tmp_c,lng_use+1,lng_use,2*lng);
01422                                 
01423                 ifft (Pxy,m+1);                                                                 /* inverse FFT of the scaled power spectrum */
01424                 
01425                 ifft (Rxx,m+1);                                                                 /*calculating autocorrelation of x*/
01426           
01427                 c_get (Pxy,Rxy,0,2*lng);                                        /* seperation of real and imaginary parts, only extract meaningful segment */
01428                 c_get (Pxy,HRxy,1,2*lng);  
01429 
01430                 if (biasrem == 1)
01431                         {
01432                                 f_mult (ubias,Rxy,Rxy,2*lng);                           /* removing bowtie artifact and setting redundant values to zero */
01433                                 f_mult (ubias,HRxy,HRxy,2*lng);
01434                         }
01435                 
01436 #ifdef DEBUG_ZEO_2              
01437                 write_float (Rxy,"dbg_Rxy.txt",2*lng);
01438                 write_float (HRxy,"dbg_HRxy.txt",2*lng);
01439                 write_float (ubias,"dbg_ubias.txt",2*lng);
01440 #endif          
01441                 
01442                 for (i=0;i<lng-1;++i)                                                   /* searching for the Zero crossing      */
01443                         {/* for i */
01444                                 if (i > maxdel) 
01445                                         {/* i > maxdel */
01446                                                 /*sprintf (buf,"Delay larger than 1/2 segment length (%d)",maxdel);
01447                                                 error_message ("hilbertdelay_V2",buf,0);*/
01448                                                 return (3);     
01449                                         }/* i > maxdel */
01450                                 
01451                                 if (HRxy[i] == 0.0)
01452                                         {/* HRxy[i] == 0.0 */   
01453                                                 
01454                                                 if (Rxy[i] > 0.0) 
01455                                                                                 *slp = 1.0;
01456                                                                 else 
01457                                                                                 *slp = -1.0;
01458                                                                                 
01459                                                 izero = (float) i;
01460                                                 *xcor = Rxy[i];                                 /* storing value of max covariance */
01461                                                 i=lng;
01462                                         }/* HRxy[i] == 0.0 */
01463                                 else 
01464                                         {/* HRxy[i] != 0.0 */
01465                                                 if ((HRxy[i] * HRxy[i+1]) < 0.0)
01466                                                         {
01467                                                                 /* the sign of slp was used to determine the sign of the cross correlation  */
01468                                                                 /* the sign of slp should be the same as that of Rxy[i] close to izero */
01469                                                                 /* moreover, I think the use of the sign of Rxy[i] is better since with no */
01470                                                                 /*subtraction performed, I am less sensitive to high freq. noise */
01471                                                                 if (Rxy[i] >= 0.0) 
01472                                                                                 *slp = 1.0;
01473                                                                 else 
01474                                                                                 *slp = -1.0;
01475                                                                 
01476                                                                 if ((*slp > 0.0)) 
01477                                                                         {
01478                                                                                 f_i = (float) (i);
01479                                                                                 f_i1 = (float) (i+1);
01480                                                                                 reg_pnt = 0.0;
01481                                                                                 linear_interp (&HRxy[i],&f_i,&HRxy[i+1],&f_i1,&reg_pnt,&izero,1);
01482 
01483                                                                                 /* find the peak of the cross correlation curve */
01484                                                                                 k = 0;
01485                                         for (j=-3;j<=3;++j)
01486                                                 {
01487                                          if (((i+j) >= 0) && ((i+j) < lng))
01488                                   {
01489                                     tmp_f_vect[k] = (float) (i+j);
01490                                     tmp_f_vect2[k] = Rxy[i+j];
01491                                     ++k;
01492                                   }
01493                                         }
01494                                         if (k > 1)
01495                                         {/* at least a 1st order interpolation must be performed */
01496                                 *xcor = Lagrange_interp (tmp_f_vect,tmp_f_vect2,izero,k);
01497                                         }
01498 
01499 
01500                                                                                 i = lng;
01501                                                                         }
01502                                                         }
01503                                         }/* HRxy[i] != 0.0 */
01504                         
01505                         }/* for i */
01506                 
01507                                 *del = izero + Dtx;             /* delay is in sample units corrected by the sampling time difference*/ 
01508                                         
01509                                 if (Rxx[0].real && Ryy[0].real)
01510                                         *xcorCoef = *xcor / sqrt (Rxx[0].real * Ryy[0].real) * *slp; /*correction for sign of cross correlation coefficient (slp = 1.0 or -1.00*/
01511                                 else
01512                                         {       
01513                                                 #ifdef DEBUG_ZEO_3              
01514                                                         printf ("\nZero Variance...\n");
01515                                                 #endif
01516                                         }
01517                                 
01518                                 /* set vx and vy */
01519                                 
01520                                 *vx = Rxx[0].real;
01521                                 *vy = Ryy[0].real;
01522                                 
01523                 ++i_call;
01524                 return (0);
01525                 
01526 }/* opt > 0 */                                                                          /* Execution mode */
01527 else if (opt == 0)
01528         {/* opt == 0 */ 
01529 
01530 #ifdef DEBUG_ZEO_3              
01531         printf ("\nCleaning Up...\n");
01532 #endif
01533 
01534                 free (fftx);                                                            /*Cleaning up used space*/
01535                 free (fftxc);
01536                 free (ffty);    
01537                 free (fftyc);
01538                 free (Px);
01539                 free (Py);
01540                 free (Pxx);
01541                 free (Pyy);
01542                 free (Pxy);
01543                 free (Rxx);
01544                 free (Ryy);
01545                 free (Rxy);
01546                 free (HRxy);
01547                 free (xp);
01548                 free (yp);
01549                 free (ubias);
01550                 free (xcpy);
01551                 free (ycpy);
01552                 free (tmp_f_vect);
01553                 free (tmp_f_vect2);
01554                 
01555                 free2D ((char **)fftyca,lng+2);
01556                 free2D ((char **)Pxya,lng+2);
01557                 free2D ((char **)Rxxa,lng+2);
01558                 free2D ((char **)Ryya,lng+2);
01559 
01560                 i_call = 0;
01561                 
01562                 *del = NoWayDelay;                      /* setting variables to out of bound values for safety */
01563                 *xcorCoef = NoWayxcorCoef;
01564                 
01565                 return (0);                                             
01566         }/* opt == 0 */
01567 
01568  return(0) ;
01569 }
01570 
01571 
01572 /*-----------------------------------------------------------------------------------*/ 
01573 
01574 static void hunwrap (float del, float fs, float T, float slp, int wrp, int unt, float *delu )
01575 {/*hunwrap*/
01576 
01577         float pi = 3.1416, tmp;
01578         
01579         
01580         if (fs > 0.0)
01581                 {       /* delay should be in seconds */
01582                         del = del / fs;
01583                 }
01584         
01585         if (T > 0.0)
01586                 {/* Period unwrapping possible */
01587                         if (slp < 0.0)
01588                                 {/* Unwrapping required to determine correct delay */
01589 
01590                                                         tmp = del * 360.0 / T;          /* from time to polar angle */
01591                                                         tmp = punwrap (tmp,1);          /* unwrap */
01592                                                         tmp = tmp + 180.0;                      /* augment by pi */
01593                                                         del = tmp * T / 360.0;          /* from polar to time */
01594                                 }
01595                         
01596                         /* Now for the case where we get negative delays although no wrapping has been
01597                                 done, namely because of the sampling time correction. */
01598                         
01599                         if (del < 0.0 || del > T)
01600                                 {
01601                                         tmp = del * 360.0 / T;          /* from time to polar angle */
01602                                         if (del < 0.0)
01603                                                 { tmp = tmp + 360.0; }
01604                                         else 
01605                                                 {
01606                                                         if (del > T)
01607                                                                 tmp = tmp - 360.0;
01608                                                 }
01609                                         del = tmp * T / 360.0;  /* from polar to time */
01610                                 }  
01611                         
01612                         if (wrp == 1)
01613                                                                 {/* map of (0-pi) to (0-pi) and (pi-2pi) to (pi-0) */
01614                                                                         tmp = del * 360.0 / T;          /* from time to polar angle */
01615                                                                         tmp = punwrap (tmp,1);          /* unwrap */
01616                                                                         del = tmp * T / 360.0;          /* from polar to time */
01617                                                                 }/* map of (0-pi) to (0-pi) and (pi-2pi) to (pi-0) */
01618 
01619                         if (unt == METH_DEGREES) del = del * 360.0 / T;         /* from time to polar angle in degrees*/
01620                         if (unt == METH_RADIANS) del = del * 2 * pi / T;                /* from time to polar angle in degrees*/        
01621         }/* Period unwrapping possible */
01622         
01623         *delu = del;
01624         
01625         return;
01626 }/*hunwrap*/
01627  
01628 /*-----------------------------------------------------------------------------------*/ 
01629 
01630 static void disp_comp_vect (COMPLEX *v,int l)
01631         {
01632                 int i;
01633 
01634                 printf ("\n");
01635                 
01636                 if ((l-1) == 0)
01637                         {
01638                                 printf ("V = (%f,%f)\n",(*v).real,(*v).imag);
01639                         }
01640                 else 
01641                 {
01642                         for (i=0;i<l;++i)
01643                         {
01644                                 printf ("V[%d] = (%f,%f)\t",i,v[i].real,v[i].imag);
01645                         }
01646                 printf ("\n");
01647                 }
01648                 return;
01649 
01650         }
01651 
01652 /*-----------------------------------------------------------------------------------*/ 
01653 
01654 static void disp_vect (float *v,int l)
01655         {
01656                 int i;
01657 
01658                 printf ("\n");
01659                 if ((l-1) == 0)
01660                         {
01661                                 printf ("V = %f\n",*v);
01662                         }
01663                 else 
01664                 {
01665                         for (i=0;i<l;++i)
01666                         {
01667                                 printf ("V[%d] = %f\t",i,v[i]);
01668                         }
01669                         printf ("\n");
01670                 }
01671                 return;
01672 
01673         }
01674 
01675 /*-----------------------------------------------------------------------------------*/ 
01676 
01677 static int is_vect_null ( float * v , int npts )
01678         {/*is_vect_null*/
01679         
01680         int i , is_ts_null;
01681                 
01682         is_ts_null = 1;                 /* Start loop in bad faith */
01683         
01684         for (i=0;i<npts;++i) 
01685                   {
01686                         if (v[i] != 0.0) 
01687                                 {       
01688                                         is_ts_null = 0; /* vector is not null */
01689                                         break;          
01690                                 }
01691                   }
01692                  
01693         return (is_ts_null);    
01694         
01695         }/*is_vect_null*/
01696  
01697 /*-----------------------------------------------------------------------------------*/ 
01698 
01699 
01700 static int write_float (float *x,char *f_name,int n_points)
01701         
01702    
01703     { /*  */
01704      int i;
01705      
01706      FILE*internal_file;
01707      
01708      internal_file = fopen (f_name,"w");
01709      if (internal_file == NULL) {
01710                                                                 printf ("\aCould not open %s \n",f_name);
01711                                                                 printf ("Exiting program\n");
01712                                                                 exit (0);
01713                                                         }
01714    
01715    for (i=0;i<n_points;++i) fprintf (internal_file,"%f\n",x[i]);  
01716      fclose (internal_file);
01717       return (i);                                                            
01718    }
01719 
01720 
01721 
01722 /* support functions declaration for main algorithm             */
01723 /* -----------------------END-----------------------------------*/