comparison tools/filt1.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

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