comparison 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 ) ;
+
+}
+
+
+