tomwalters@0: /*
tomwalters@0:   smooth.c      Low-pass filter by windowing in the time or frequency domain.
tomwalters@0:   --------
tomwalters@0: 
tomwalters@0: In the time domain the input is convolved with a unit area Gaussian window.
tomwalters@0: In the frequency domain the spectrum is multiplied by a Gaussian mask which
tomwalters@0: zeroes higher frequency components. The transform between domains uses the FFT.
tomwalters@0: 
tomwalters@0: One parameter controls the amount of smoothing. This parameter controls the
tomwalters@0: size of the window.
tomwalters@0: 
tomwalters@0: Caveats for the frequency-domain approach are descibed in NR p436 (such as
tomwalters@0: the damped ringing effect which results when when sharp edges are smoothed).
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: #include "sigproc.h"
tomwalters@0: 
tomwalters@0: char applic[]     = "Data smoothing by convolution with Gaussian window." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *sampstr, *lenstr, *varstr, *domstr  ;
tomwalters@0: static char *ranstr,  *typestr,  *sizestr ;
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:     {   "length"    ,   "max"       ,  &lenstr      ,   "Size of input."                        , VAL     },
tomwalters@0:     {   "variance"  ,   "20"        ,  &varstr      ,   "Variance of Gaussian window"           , VAL     },
tomwalters@0:     {   "range"     ,   "4"         ,  &ranstr      ,   "Range in standard deviations about mean", VAL     },
tomwalters@0:     {   "domain"    ,   "time"      ,  &domstr      ,   "Convolve in time or frequency domain"  , VAL     },
tomwalters@0:     {   "type"      ,   "short"     ,  &typestr     ,   "Data type input and output"            , VAL     },
tomwalters@0:     {   "SIZE"      ,   "262144"    ,  &sizestr     ,   "buffer size (s, ms, samples)"          , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     samplerate  ;
tomwalters@0: int     length      ;
tomwalters@0: 
tomwalters@0: float  *window ;
tomwalters@0: int     points ;
tomwalters@0: 
tomwalters@0: short  *buf1   ;
tomwalters@0: float  *buf2   ;
tomwalters@0: float  *out    ;
tomwalters@0: FILE   *fp     ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char  **argv;
tomwalters@0: {
tomwalters@0:     int   maxlen, i ;
tomwalters@0: 
tomwalters@0:     fp = openopts( option,argc,argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts( helpstr, argv[0], applic, option ) ;
tomwalters@0: 
tomwalters@0:     samplerate = to_Hz( sampstr ) ;
tomwalters@0: 
tomwalters@0:     if ( ismax( lenstr ) )  maxlen = to_p( sizestr, samplerate ) ;
tomwalters@0:     else                    maxlen = to_p( lenstr,  samplerate ) ;
tomwalters@0: 
tomwalters@0:     out    = (float *)malloc( maxlen * sizeof(float) ) ;
tomwalters@0:     window = gauss_window( to_p( sqrt_units( varstr ), samplerate ), atof( ranstr ), &points ) ;
tomwalters@0:     normalize_area( window, points ) ;
tomwalters@0: 
tomwalters@0:     if ( isstr( typestr, "short" ) ) {
tomwalters@0: 	buf1 = (short *)malloc( maxlen * sizeof(short) ) ;
tomwalters@0: 	length = fread( buf1, sizeof(short), maxlen, fp ) ;
tomwalters@0:     }
tomwalters@0:     else if ( isstr( typestr, "float" ) ) {
tomwalters@0: 	buf2 = (float *)malloc( maxlen * sizeof(float) ) ;
tomwalters@0: 	length = fread( buf2, sizeof(float), maxlen, fp ) ;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	fprintf( stderr,"smooth: unknown data type [%s]\n", typestr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( length < maxlen ) {
tomwalters@0: 	if ( length == 0 ) {
tomwalters@0: 	    fprintf( stderr,"smooth: missing input data\n");
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 	if ( !ismax( lenstr ) )
tomwalters@0: 	    fprintf( stderr,"smooth: warning: %d samples found\n", length ) ;
tomwalters@0:     }
tomwalters@0:     fclose( fp ) ;
tomwalters@0: 
tomwalters@0:     if ( points <= 0 ) {
tomwalters@0: 	fprintf( stderr,"smooth: window size = %d\n", points ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( isstr( typestr, "short" ) ) {
tomwalters@0: 
tomwalters@0: 	if ( points == 1 )  /* no smoothing */
tomwalters@0: 	    fwrite( buf1, sizeof(short), length, stdout ) ;
tomwalters@0: 
tomwalters@0: 	else if ( iststr( domstr, "time" ) )
tomwalters@0: 	    convolve_time( buf1, length, window, points, out ) ;
tomwalters@0: 
tomwalters@0: 	else if ( iststr( domstr, "frequency" ) )
tomwalters@0: 	    convolve_freq( buf1, length, window, points, out ) ;
tomwalters@0: 
tomwalters@0: 	else {
tomwalters@0: 	    fprintf( stderr,"smooth: unknown domain [%s]\n", domstr ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	ftos( out, buf1, length, 1. ) ;
tomwalters@0: 	fwrite( buf1, sizeof(short), length, stdout ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     else {
tomwalters@0: 
tomwalters@0: 	if ( points == 1 )  /* no smoothing */
tomwalters@0: 	    fwrite( buf2, sizeof(float), length, stdout ) ;
tomwalters@0: 
tomwalters@0: 	else if ( iststr( domstr, "time" ) )
tomwalters@0: 	    convolve_time2( buf2, length, window, points, out ) ;
tomwalters@0: 
tomwalters@0: 	else if ( iststr( domstr, "frequency" ) )
tomwalters@0: 	    convolve_freq2( buf2, length, window, points, out ) ;
tomwalters@0: 
tomwalters@0: 	else {
tomwalters@0: 	    fprintf( stderr,"smooth: unknown domain [%s]\n", domstr ) ;
tomwalters@0: 	    exit( 1 ) ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	fwrite( out, sizeof(float), length, stdout ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: convolve_time2( y, n, x, m, z )
tomwalters@0: float  *y ;             /* input signal                                  */
tomwalters@0: int     n ;             /* length of array y                             */
tomwalters@0: float  *x ;             /* convolution window or filter impulse response */
tomwalters@0: int     m ;             /* length of array x                             */
tomwalters@0: float  *z ;             /* output signal of length n (same as input)     */
tomwalters@0: {
tomwalters@0:     register int  i, j, k, lim1, lim2 ;
tomwalters@0:     register float sum ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ ) {
tomwalters@0: 
tomwalters@0: 	k = m - 1 ;
tomwalters@0: 	if ( ( lim1 = i - (m-1)/2 ) < 0 ) {
tomwalters@0: 	    k += lim1 ;
tomwalters@0: 	    lim1 = 0 ;
tomwalters@0: 	}
tomwalters@0: 	if ( ( lim2 = i + (m-1)/2 ) > n )
tomwalters@0: 	    lim2 = n - 1 ;
tomwalters@0: 
tomwalters@0: 	sum = 0 ;
tomwalters@0: 	for ( j = lim1 ; j <= lim2 ; j++, --k )
tomwalters@0: 	    sum += x[k] * y[j] ;
tomwalters@0: 
tomwalters@0: 	z[i] = sum ;
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Frequency-domain convolution.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: convolve_freq2( y, n, x, m, z )
tomwalters@0: float  *y ;             /* input signal                                  */
tomwalters@0: int     n ;             /* length of array y                             */
tomwalters@0: float  *x ;             /* convolution window or filter impulse response */
tomwalters@0: int     m ;             /* length of array x                             */
tomwalters@0: float  *z ;             /* output signal of length n (same as input)     */
tomwalters@0: {
tomwalters@0:     int     i, j, n1 ;
tomwalters@0:     float  *data, *respns, *ans ;
tomwalters@0: 
tomwalters@0:     n1 = getpower( n + m ) ;    /* padded length of input signal */
tomwalters@0: 
tomwalters@0:     data   = (float *)malloc( ( n1     + 1 ) * sizeof(float) ) ;
tomwalters@0:     respns = (float *)malloc( ( n1     + 1 ) * sizeof(float) ) ;
tomwalters@0:     ans    = (float *)malloc( ( n1 * 2 + 1 ) * sizeof(float) ) ;
tomwalters@0: 
tomwalters@0:     /* copy and pad signal */
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	data[i] = y[i] ;
tomwalters@0:     for (  ; i < n1 ; i++ )
tomwalters@0: 	data[i] = (float)0 ;
tomwalters@0: 
tomwalters@0:     /* copy window into wrap-around order */
tomwalters@0: 
tomwalters@0:     for ( i = m/2, j = 0 ; i < m ; i++, j++ )
tomwalters@0: 	respns[i] = x[j] ;
tomwalters@0:     for ( i = 0 ; i < m/2 ; i++, j++ )
tomwalters@0: 	respns[i] = x[j] ;
tomwalters@0: 
tomwalters@0:     /* Convolve and copy output */
tomwalters@0: 
tomwalters@0:     convlv( data-1, n1, respns-1, m, 1, ans-1 ) ;
tomwalters@0: 
tomwalters@0:     for ( i = 0 ; i < n ; i++ )
tomwalters@0: 	z[i] = ans[i] ;
tomwalters@0: 
tomwalters@0:     free( data   ) ;
tomwalters@0:     free( respns ) ;
tomwalters@0:     free( ans    ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: