view 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 ) ;
}