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

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

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