view tools/conv.c @ 0:5242703e91d3 tip

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

/*
  conv.c  Convolution.  [NR: 407-413].
  ------

  Usage:  conv [options]  signal_file  response_file

  The response of a linear filter to an arbitiary input signal is the
  convolution of the signal with the filter's impulse response.
  The `response_file' contains the impulse response which characterises the
  filter, and the result of the convolution operation is the response of
  that filter to the signal contained in the `signal_file'.

  The program uses the FFT to transform both input files, and then obtains
  the convolution by multiplying in the frequency domain, and inverse
  transforming the result.

Eg:

ptrain > foo1                               -signal
gauss -b5 | op type=float mult=500 > foo3   -response
conv factor=.0005 foo1 foo3 | x11plot

Deconv doesn't seem to work !
*/


#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <sys/dir.h>
#include <sys/stat.h>
#include "options.h"
#include "units.h"
#include "strmatch.h"
#include "sigproc.h"

char applic[]     = "Convolution. \n i/p and o/p signal in binary shorts, impulse response file in floats." ;
char usage[]      = "conv [options] signal_file  response_file" ;

static char *helpstr, *debugstr, *sampstr, *lenstr, *domstr,  *scalestr ;

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     },
    {   "domain"    ,   "time"        ,  &domstr      ,   "Convolve in time or frequency domain"  , VAL     },
    {   "scale"     ,   "1.0"         ,  &scalestr    ,   "scale factor for output"               , VAL     },
   ( char * ) 0 } ;


int     samplerate  ;
float   scale       ;

short  *signal      ;
float  *response    ;
float  *out         ;

main (argc, argv)
int     argc;
char  **argv;
{
    FILE     *fp1, *fp2  ;
    struct    stat  stbuf ;
    int       i, n, m ;
    float     f ;

    if ( getopts( option, argc, argv ) != 2 || !isoff( helpstr ) )
	helpopts3( helpstr, argv[0], applic, usage, option ) ;

    samplerate = to_Hz( sampstr ) ;
    scale      = atof( scalestr ) ;

    if ( ( fp1 = fopen( argv[argc-2], "r" ) ) == (FILE *)0 ) {
	fprintf( stderr,"%s: can't open %s\n", argv[0], argv[argc-2] ) ;
	exit( 1 ) ;
    }
    stat( argv[argc-2], &stbuf ) ;
    if ( ismax( lenstr ) )
	n = stbuf.st_size / sizeof( short ) ;       /* size of signal file */
    else if ( ( n = to_p( lenstr,  samplerate ) ) > stbuf.st_size / sizeof( short ) ) {
	fprintf( stderr,"conv: insufficient signal\n");
	exit( 1 ) ;
    }


    if ( ( fp2 = fopen( argv[argc-1], "r" ) ) == (FILE *)0 ) {
	fprintf( stderr,"%s: can't open %s\n", argv[0], argv[argc-1] ) ;
	exit( 1 ) ;
    }
    stat( argv[argc-1], &stbuf ) ;
    m = stbuf.st_size / sizeof( float ) ;       /* size of response file */

    if ( m > n ) {      /* size of response file must <= signal file */
	fprintf( stderr,"conv: size of response file must be <= size of signal file\n" ) ;
	exit( 1 ) ;
    }

    signal   = (short *)malloc( n * sizeof(short) ) ;
    response = (float *)malloc( m * sizeof(float) ) ;
    out      = (float *)malloc( n * sizeof(float) ) ;

    fread( signal,   sizeof(short), n, fp1 ) ;
    fread( response, sizeof(float), m, fp2 ) ;
    fclose( fp1 ) ; fclose( fp2 ) ;


    if ( iststr( domstr, "time" ) )
	convolve_time( signal, n, response, m, out ) ;

    else if ( iststr( domstr, "frequency" ) )
	convolve_freq( signal, n, response, m, out ) ;

    else {
	fprintf( stderr,"conv: unknown domain [%s]\n", domstr ) ;
	exit( 1 ) ;
    }


    if ( ( f = ftos( out, signal, n, scale ) ) < 1. )
	fprintf( stderr,"Warning: 16-bit overflow during convolution. Try scale factor < %f\n", f ) ;

    fwrite( signal, sizeof(short), n, stdout ) ;

}