tomwalters@0: /*
tomwalters@0:   filt1.c       1st order LP filter using an exponential smoother
tomwalters@0: 		(recursive "leaky integrator" filter)
tomwalters@0: 
tomwalters@0:   The n'th recursive update is:
tomwalters@0: 
tomwalters@0:     a)    y[n] = a.y[n-1] + x[n]                recursive sum
tomwalters@0:     b)    y[n] = a.y[n-1] + (1-a).x[n]          recursive mean
tomwalters@0: 
tomwalters@0:   where decay constant a = exp(-Ts/T)
tomwalters@0:   and   Ts is sample interval in seconds
tomwalters@0: 	T  is decay time-constant in seconds
tomwalters@0: 
tomwalters@0:   The 1st order LP filter is an exponential window which decays into the past.
tomwalters@0:   With decay constant a=1 there is infinite memory because it does not decay.
tomwalters@0:   The result is either a steadily accumulating sum or mean.
tomwalters@0:   (A zero mean is assumed).
tomwalters@0:   With 0 < a < 1 there is a finite decaying memory.
tomwalters@0:   The smaller `a', the smaller the memory (it decays faster).
tomwalters@0:   With a = 0 the filter output is identical with its input.
tomwalters@0: 
tomwalters@0:   The window width parameter is the time-constant parameter.
tomwalters@0:   The half life of the exponential window is given by -T.ln(0.5) = 0.693T
tomwalters@0:   ie. for a given time constant, T secs, the window decays to half of its
tomwalters@0:   current value in a period of 0.693T secs.
tomwalters@0:   (Thus T can be set to accomodate an expected period).
tomwalters@0: 
tomwalters@0:   When the time-constant parameter is given with time units (s or ms) then
tomwalters@0:   it is converted to a decay constant using  a = exp(-Ts/T).
tomwalters@0:   Otherwise it is taken to be the decay constant directly.
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[]     = "1st order low-pass filter using exponential smoother." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *sampstr, *tcstr,  *destr  ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                  , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"      , DEBUG   },
tomwalters@0:     {   "samplerate",   "20kHz"     ,  &sampstr     ,   "Samplerate "           , VAL     },
tomwalters@0:     {   "tc"        ,   "3ms"       ,  &tcstr       ,   "Decay time constant"   , VAL     },
tomwalters@0:     {   "de"        ,   "mean"      ,  &destr       ,   "Difference equation"   , VAL     },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     samplerate  ;
tomwalters@0: double  a           ;   /* decay constant */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char  **argv;
tomwalters@0: {
tomwalters@0:     FILE  *fp  ;
tomwalters@0:     short  x, y ;
tomwalters@0:     int    morehelp() ;
tomwalters@0: 
tomwalters@0:     fp = openopts( option,argc,argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts1( helpstr, argv[0], applic, option, morehelp ) ;
tomwalters@0: 
tomwalters@0:     samplerate = to_Hz( sampstr ) ;
tomwalters@0: 
tomwalters@0:     if ( Units( tcstr ) )
tomwalters@0: 	a = exp( -(double)( 1. / ( samplerate * to_s( tcstr, samplerate ) ) ) ) ;
tomwalters@0:     else
tomwalters@0: 	a = (double)atof( tcstr ) ;
tomwalters@0: 
tomwalters@0:     y = 0 ;     /* initialization */
tomwalters@0: 
tomwalters@0:     if ( iststr( destr, "mean" ) )
tomwalters@0: 	while ( fread( &x, sizeof(short), 1, fp ) ) {
tomwalters@0: 	    y = a * y  +  ( 1 - a ) * x ;
tomwalters@0: 	    fwrite( &y, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     else if ( iststr( destr, "sum" ) )
tomwalters@0: 	while ( fread( &x, sizeof(short), 1, fp ) ) {
tomwalters@0: 	    y = a * y  +  x ;
tomwalters@0: 	    fwrite( &y, sizeof(short), 1, stdout ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     else {
tomwalters@0: 	fprintf( stderr,"unknown de [%s]\n", destr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     fclose( fp );
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* Return 1 if the str has appended units. Otherwise return 0.             */
tomwalters@0: 
tomwalters@0: Units(str)
tomwalters@0: char *str ;
tomwalters@0: {
tomwalters@0:     char *eptr = str + strlen( str ) ;
tomwalters@0: 
tomwalters@0:     if( isdigit( *--eptr ) ) return 0 ; /* last char is digit, so no units */
tomwalters@0:     else return 1 ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: morehelp()
tomwalters@0: {
tomwalters@0:     fprintf(stderr,"\nde: \n");
tomwalters@0:     fprintf(stderr,"    sum     y[n] = a.y[n-1] + x[n]   \n");
tomwalters@0:     fprintf(stderr,"    mean    y[n] = a.y[n-1] + (1-a).x[n]    (unity steady-state gain)\n");
tomwalters@0:     exit( 1 ) ;
tomwalters@0: }