tomwalters@0: /*
tomwalters@0:   conv.c  Convolution.  [NR: 407-413].
tomwalters@0:   ------
tomwalters@0: 
tomwalters@0:   Usage:  conv [options]  signal_file  response_file
tomwalters@0: 
tomwalters@0:   The response of a linear filter to an arbitiary input signal is the
tomwalters@0:   convolution of the signal with the filter's impulse response.
tomwalters@0:   The `response_file' contains the impulse response which characterises the
tomwalters@0:   filter, and the result of the convolution operation is the response of
tomwalters@0:   that filter to the signal contained in the `signal_file'.
tomwalters@0: 
tomwalters@0:   The program uses the FFT to transform both input files, and then obtains
tomwalters@0:   the convolution by multiplying in the frequency domain, and inverse
tomwalters@0:   transforming the result.
tomwalters@0: 
tomwalters@0: Eg:
tomwalters@0: 
tomwalters@0: ptrain > foo1                               -signal
tomwalters@0: gauss -b5 | op type=float mult=500 > foo3   -response
tomwalters@0: conv factor=.0005 foo1 foo3 | x11plot
tomwalters@0: 
tomwalters@0: Deconv doesn't seem to work !
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include <sys/types.h>
tomwalters@0: #include <sys/dir.h>
tomwalters@0: #include <sys/stat.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[]     = "Convolution. \n i/p and o/p signal in binary shorts, impulse response file in floats." ;
tomwalters@0: char usage[]      = "conv [options] signal_file  response_file" ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *sampstr, *lenstr, *domstr,  *scalestr ;
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:     {   "domain"    ,   "time"        ,  &domstr      ,   "Convolve in time or frequency domain"  , VAL     },
tomwalters@0:     {   "scale"     ,   "1.0"         ,  &scalestr    ,   "scale factor for output"               , VAL     },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: int     samplerate  ;
tomwalters@0: float   scale       ;
tomwalters@0: 
tomwalters@0: short  *signal      ;
tomwalters@0: float  *response    ;
tomwalters@0: float  *out         ;
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char  **argv;
tomwalters@0: {
tomwalters@0:     FILE     *fp1, *fp2  ;
tomwalters@0:     struct    stat  stbuf ;
tomwalters@0:     int       i, n, m ;
tomwalters@0:     float     f ;
tomwalters@0: 
tomwalters@0:     if ( getopts( option, argc, argv ) != 2 || !isoff( helpstr ) )
tomwalters@0: 	helpopts3( helpstr, argv[0], applic, usage, option ) ;
tomwalters@0: 
tomwalters@0:     samplerate = to_Hz( sampstr ) ;
tomwalters@0:     scale      = atof( scalestr ) ;
tomwalters@0: 
tomwalters@0:     if ( ( fp1 = fopen( argv[argc-2], "r" ) ) == (FILE *)0 ) {
tomwalters@0: 	fprintf( stderr,"%s: can't open %s\n", argv[0], argv[argc-2] ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     stat( argv[argc-2], &stbuf ) ;
tomwalters@0:     if ( ismax( lenstr ) )
tomwalters@0: 	n = stbuf.st_size / sizeof( short ) ;       /* size of signal file */
tomwalters@0:     else if ( ( n = to_p( lenstr,  samplerate ) ) > stbuf.st_size / sizeof( short ) ) {
tomwalters@0: 	fprintf( stderr,"conv: insufficient signal\n");
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( ( fp2 = fopen( argv[argc-1], "r" ) ) == (FILE *)0 ) {
tomwalters@0: 	fprintf( stderr,"%s: can't open %s\n", argv[0], argv[argc-1] ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     stat( argv[argc-1], &stbuf ) ;
tomwalters@0:     m = stbuf.st_size / sizeof( float ) ;       /* size of response file */
tomwalters@0: 
tomwalters@0:     if ( m > n ) {      /* size of response file must <= signal file */
tomwalters@0: 	fprintf( stderr,"conv: size of response file must be <= size of signal file\n" ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     signal   = (short *)malloc( n * sizeof(short) ) ;
tomwalters@0:     response = (float *)malloc( m * sizeof(float) ) ;
tomwalters@0:     out      = (float *)malloc( n * sizeof(float) ) ;
tomwalters@0: 
tomwalters@0:     fread( signal,   sizeof(short), n, fp1 ) ;
tomwalters@0:     fread( response, sizeof(float), m, fp2 ) ;
tomwalters@0:     fclose( fp1 ) ; fclose( fp2 ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( iststr( domstr, "time" ) )
tomwalters@0: 	convolve_time( signal, n, response, m, out ) ;
tomwalters@0: 
tomwalters@0:     else if ( iststr( domstr, "frequency" ) )
tomwalters@0: 	convolve_freq( signal, n, response, m, out ) ;
tomwalters@0: 
tomwalters@0:     else {
tomwalters@0: 	fprintf( stderr,"conv: unknown domain [%s]\n", domstr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     if ( ( f = ftos( out, signal, n, scale ) ) < 1. )
tomwalters@0: 	fprintf( stderr,"Warning: 16-bit overflow during convolution. Try scale factor < %f\n", f ) ;
tomwalters@0: 
tomwalters@0:     fwrite( signal, sizeof(short), n, stdout ) ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: