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 }

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         }

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         }

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         }

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         }

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         }

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         }

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         }

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         }

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  }

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   }

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         }

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 }

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    }

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          }

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 }

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         } 

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 }

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*/

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 }

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*/

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*/

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         }

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      }

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 */

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 ()*/

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         }

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*/

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    }

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 }

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    }

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*/