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