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 #include <stdio.h>
00008 #include <math.h>
00009 #include <stdlib.h>
00010 #include <string.h>
00011 
00012 #include "ts.h"
00013 
00014 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00015 
00016 /*** free a time_series structure ***/
00017 
00018 void RWC_free_time_series(ts)
00019      time_series *ts;
00020 {
00021    if( ts != NULL ){
00022       if( ts->fname != NULL ) free( ts->fname ) ;
00023       if( ts->ts    != NULL ) free( ts->ts ) ;
00024       free( ts ) ;
00025    }
00026 }
00027 
00028 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00029 
00030 /*** make a blank time_series ***/
00031 
00032 time_series * RWC_blank_time_series( npt )
00033   int npt ;
00034 {
00035    int ii ;
00036    time_series *vec ;
00037 
00038    if( npt <= 0 ) return NULL ;
00039 
00040    vec = (time_series *) malloc( sizeof(time_series) ) ;
00041    if( vec == NULL ) MALLOC_ERR("new time_series") ;
00042 
00043    vec->ts = (float *) malloc( sizeof(float) * npt ) ;
00044    if( vec->ts == NULL ) MALLOC_ERR("new time_series vector") ;
00045 
00046    vec->fname = NULL ;   /* user must supply name */
00047    vec->len   = npt ;
00048 
00049    for( ii=0 ; ii < npt ; ii++ ) vec->ts[ii] = 0.0 ;
00050    return vec ;
00051 }
00052 
00053 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00054 
00055 /*** median filter a time series ***/
00056 
00057 void RWC_medfilt_time_series( vec )
00058    time_series * vec ;
00059 {
00060    int ii , npt ;
00061    float aa , bb , cc , mtemp ;
00062    float * id ;
00063 
00064 #define MEDIAN(a,b,c) ( mtemp = 4*((a)<(b)) + 2*((a)<(c)) + ((b)<(c)) , \
00065                         (mtemp==3||mtemp==4) ? (a) :                    \
00066                         (mtemp==7||mtemp==0) ? (b) : (c) )
00067 
00068    if( vec == NULL || (npt = vec->len) < 3 ) return ;
00069 
00070    id = vec->ts ;
00071 
00072    bb = id[0] ; cc = id[1] ;
00073    for( ii=1 ; ii < npt-1 ; ii++ ){
00074       aa = bb ; bb = cc ; cc = id[ii+1] ;
00075       id[ii] = MEDIAN(aa,bb,cc) ;
00076    }
00077 
00078    return ;
00079 }
00080 
00081 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00082 
00083 /*** read a time_series structure from a file ***/
00084 
00085 #define INC_TSARSIZE 256
00086 
00087 time_series * RWC_read_time_series(fname)
00088      char *fname;
00089 {
00090    int ii , used_tsar , alloc_tsar ;
00091    float ftemp ;
00092    time_series *vec ;
00093    FILE *fts ;
00094    char *cpt , linbuf[1024] ; /* 07 Dec 2003 */
00095 
00096    vec = (time_series *) malloc( sizeof(time_series) ) ;
00097    if( vec == NULL ) MALLOC_ERR("new time_series") ;
00098 
00099    fts = fopen( fname , "r" ) ;
00100    if( fts == NULL ){
00101       fprintf( stderr , "cannot open time series file %s\n" , fname ) ;
00102       free( vec ) ;
00103       return NULL ;
00104    }
00105 
00106    ii = strlen( fname ) + 1 ;
00107    vec->fname = (char *) malloc( sizeof(char) * ii ) ;
00108    if( vec->fname == NULL ) MALLOC_ERR("new time_series fname") ;
00109    strcpy( vec->fname , fname ) ;
00110 
00111    used_tsar  = 0 ;
00112    alloc_tsar = INC_TSARSIZE ;
00113    vec->ts    = (float *) malloc( sizeof(float) * alloc_tsar ) ;
00114    if( vec->ts == NULL ) MALLOC_ERR("new time_series initial ts") ;
00115 
00116    while( 1 ){
00117 
00118 #if 1
00119       cpt = fgets(linbuf,1023,fts) ;     /* 07 Dec 2003: read line buffer */
00120       if( cpt == NULL ) break ;                     /* bad read */
00121       cpt = linbuf ;
00122       while( isspace(*cpt) ) cpt++ ;                /* skip leading blanks */
00123       if( *cpt == '#' || *cpt == '\0' ) continue ;  /* skip comments */
00124       ii = sscanf( cpt , "%f" , &ftemp ) ;          /* read value from buffer */
00125 #else
00126       ii = fscanf( fts , "%f" , &ftemp ) ;          /* Ye Olde Way */
00127 #endif
00128       if( ii != 1 ) break ;
00129 
00130       if( used_tsar == alloc_tsar ){
00131          alloc_tsar += INC_TSARSIZE ;
00132          vec->ts     = (float *)realloc( vec->ts,sizeof(float)*alloc_tsar );
00133          if( vec->ts == NULL ) MALLOC_ERR("expanding new time_series ts") ;
00134       }
00135 
00136       vec->ts[used_tsar++] = ftemp ;
00137 
00138    }
00139 
00140    fclose(fts) ;  /* Oopsie - 09 Apr 2004! */
00141 
00142    if( used_tsar <= 0 ){
00143       fprintf( stderr , "cannot read data from time_series %s\n" , fname ) ;
00144       free( vec->fname ) ;
00145       free( vec->ts ) ;
00146       free( vec ) ;
00147       return NULL ;
00148    }
00149 
00150    vec->len = used_tsar ;
00151    vec->ts  = (float *) realloc( vec->ts , sizeof(float) * used_tsar ) ;
00152    if( vec->ts == NULL ) MALLOC_ERR("shrinking new time_series ts") ;
00153 
00154    return vec ;
00155 }
00156 
00157 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00158 
00159 /*** compute the L1 norm of a time_series ***/
00160 
00161 float RWC_norm_ts(nn, vec)
00162      int nn;
00163      time_series *vec;
00164 {
00165    int ii , itop = MIN(nn,vec->len) ;
00166    float *var = vec->ts ;
00167    float vsum = 0.0 ;
00168 
00169    for( ii=0 ; ii < itop ; ii++ ) vsum += fabs(var[ii]) ;
00170 
00171    return vsum ;
00172 }
00173 
00174 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00175 
00176 /*** compute the max of a time_series ***/
00177 
00178 float RWC_max_ts(nn, vec)
00179      int nn;
00180      time_series *vec;
00181 {
00182    int ii , itop = MIN(nn,vec->len) ;
00183    float *var = vec->ts ;
00184    float vsum ;
00185 
00186    vsum = var[0] ;
00187    for( ii=1 ; ii < itop ; ii++ ) vsum = MAX(vsum,var[ii]) ;
00188    return vsum ;
00189 }
00190 
00191 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00192 
00193 /*** compute the min of a time_series ***/
00194 
00195 float RWC_min_ts(nn, vec)
00196      int nn;
00197      time_series *vec;
00198 {
00199    int ii , itop = MIN(nn,vec->len) ;
00200    float *var = vec->ts ;
00201    float vsum ;
00202 
00203    vsum = var[0] ;
00204    for( ii=1 ; ii < itop ; ii++ ) vsum = MIN(vsum,var[ii]) ;
00205    return vsum ;
00206 }