tomwalters@0: /*
tomwalters@0:     Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
tomwalters@0:     ===========================================================================
tomwalters@0: 
tomwalters@0:     Permission to use, copy, modify, and distribute this software without fee
tomwalters@0:     is hereby granted for research purposes, provided that this copyright 
tomwalters@0:     notice appears in all copies and in all supporting documentation, and that 
tomwalters@0:     the software is not redistributed for any fee (except for a nominal shipping 
tomwalters@0:     charge). Anyone wanting to incorporate all or part of this software in a
tomwalters@0:     commercial product must obtain a license from the Medical Research Council.
tomwalters@0: 
tomwalters@0:     The MRC makes no representations about the suitability of this 
tomwalters@0:     software for any purpose.  It is provided "as is" without express or implied 
tomwalters@0:     warranty.
tomwalters@0:  
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:     bank.c
tomwalters@0:     ======
tomwalters@0: 
tomwalters@0:     long awaited proper generic multiplexed filter bank with phase compensation!
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     Copyright (c), 1989  The Medical Research Council, Applied Psychology Unit.
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     Author  : John Holdsworth
tomwalters@0:     Written : 10th March, 1989.
tomwalters@0: 
tomwalters@0:     Edited  :
tomwalters@0:               Christian Giguere, March 1994
tomwalters@0:               - Provided facilities to print filterbank information to stderr
tomwalters@0:               - To locate changes, search for "CG"
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: #include <stdio.h>       /* CG */
tomwalters@0: 
tomwalters@0: #include "stitch.h"
tomwalters@0: #include "source.h"
tomwalters@0: #include "recurse.h"
tomwalters@0: #include "formulae.h"
tomwalters@0: 
tomwalters@0: #include "bank.h"
tomwalters@0: 
tomwalters@0: #ifndef  lint
tomwalters@0: static char *sccs_id = "@(#)bank.c	1.16 John Holdsworth (MRC-APU) 6/6/91" ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* interpolation */
tomwalters@0: 
tomwalters@0: struct _interp_state { struct _fillable_source parent ; Source source ; unsigned ichans ; int (*interp)() ; } ;
tomwalters@0: 
tomwalters@0: static Pointer interp_callback( state, bytes, buffer )
tomwalters@0: struct _interp_state *state ;
tomwalters@0: ByteCount *bytes ;
tomwalters@0: short *buffer ;
tomwalters@0: {
tomwalters@0:     register int last = *bytes == 0 ;
tomwalters@0:     register int ipointer, points = ToPoints( short, *bytes ) ;
tomwalters@0:     register int ipoints = points / (state->ichans*2-1) * state->ichans ;
tomwalters@0:     short *input = PullShorts( state->source, ipoints ) ;
tomwalters@0: 
tomwalters@0:     if( !last ) {
tomwalters@0: 	for( ipointer=0 ; ipointer < ipoints ; ipointer += state->ichans )
tomwalters@0: 	    state->interp( input+ipointer, buffer+ipointer/state->ichans*(state->ichans*2-1), state->ichans, state->ichans*2-1 ) ;
tomwalters@0: 
tomwalters@0: 	return ( (Pointer) buffer ) ;
tomwalters@0:     }
tomwalters@0:     else
tomwalters@0: 	return ( DeleteFillableSource( state ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: Source InterpSource( source, ichans )
tomwalters@0: Source source ;
tomwalters@0: int ichans ;
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _interp_state *, state ) ;
tomwalters@0:     extern int interp() ;
tomwalters@0: 
tomwalters@0:     state->source = source ;
tomwalters@0:     state->ichans = ichans ;
tomwalters@0:     state->interp = interp ;
tomwalters@0: 
tomwalters@0:     return ( SetFillableSource( state, interp_callback, "bank.c interpolating" ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: Source LinterpSource( source, ichans )
tomwalters@0: Source source ;
tomwalters@0: int ichans ;
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _interp_state *, state ) ;
tomwalters@0:     extern int linterp() ;
tomwalters@0: 
tomwalters@0:     state->source = source ;
tomwalters@0:     state->ichans = ichans ;
tomwalters@0:     state->interp = linterp ;
tomwalters@0: 
tomwalters@0:     return ( SetFillableSource( state, interp_callback, "bank.c linterpolating" ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: static int segsize = 50 ;
tomwalters@0: 
tomwalters@0: #ifdef DSP32
tomwalters@0: #define BUFFSIZE 256
tomwalters@0: #else
tomwalters@0: #define BUFFSIZE 2048
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: /* multiplexed filter bank */
tomwalters@0: 
tomwalters@0: struct _bank_state {
tomwalters@0:  struct _fillable_source parent ;
tomwalters@0:  Source source ; char *input ; unsigned chans ;
tomwalters@0:  int active, pointer, iptr, *delays, (*proc)(), input_size, output_size ;
tomwalters@0:  RecursiveFilterState **states ;
tomwalters@0:  } ;
tomwalters@0: 
tomwalters@0: static Pointer bank_callback( state, bytes, buffer )
tomwalters@0: struct _bank_state *state ;
tomwalters@0: ByteCount *bytes ;
tomwalters@0: char *buffer ;
tomwalters@0: {
tomwalters@0: extern char *bstart, *bend ;
tomwalters@0:     register int last = *bytes == 0 ;
tomwalters@0:     static int segment = BUFFSIZE ;
tomwalters@0:     int points, bytecount ;
tomwalters@0: 
tomwalters@0:     if( state->input == ( char * ) 0 ) {
tomwalters@0: 
tomwalters@0: 	if( state->delays[0] < 0 )
tomwalters@0: 	    state->pointer = state->delays[0] ;
tomwalters@0: 	else
tomwalters@0: 	    state->pointer = 0 ;
tomwalters@0: 
tomwalters@0: 	if( segment < state->delays[state->chans-1] - state->delays[0] )
tomwalters@0: 	    segment = state->delays[state->chans-1] - state->delays[0] << 1 ;
tomwalters@0: 
tomwalters@0: 	state->input = Pull( state->source, segment * state->input_size ) ;
tomwalters@0: 	state->iptr   = state->pointer ;
tomwalters@0: 	state->active = 0 ;
tomwalters@0: 
tomwalters@0: 	bstart = state->input - segment * state->input_size ;
tomwalters@0: 	bend   = state->input + segment * state->input_size ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if( !last ) {
tomwalters@0: 	for( bytecount=0 ; bytecount<*bytes ; ) {
tomwalters@0: 
tomwalters@0: 	    if( state->iptr + segsize - state->delays[0] >= segment ) {
tomwalters@0: 		state->input = Pull( state->source, segment * state->input_size ) ;
tomwalters@0: 		state->iptr -= segment ;
tomwalters@0: 
tomwalters@0: 	    bstart = state->input - segment * state->input_size ;
tomwalters@0: 	    bend   = state->input + segment * state->input_size ;
tomwalters@0: 
tomwalters@0: 	    }
tomwalters@0: 
tomwalters@0: 	    while( state->active < state->chans && state->delays[ state->active ] <= state->pointer )
tomwalters@0: 		state->active++ ;
tomwalters@0: 
tomwalters@0: 	    if( state->pointer < 0 || state->active < state->chans || state->output_size == 0 ) {
tomwalters@0: 
tomwalters@0: 		points = 1 ;
tomwalters@0: 
tomwalters@0: 		state->output_size = state->proc( state->states, state->delays, state->input + state->iptr * state->input_size, buffer+bytecount, points, state->active ) ;
tomwalters@0: 
tomwalters@0: 		stitch_bzero( buffer+bytecount+state->active*state->output_size, (int) ( ( state->chans - state->active ) * state->output_size ) ) ;
tomwalters@0: 	    }
tomwalters@0: 	    else {
tomwalters@0: 		points = ( *bytes - bytecount ) / state->chans / state->output_size ;
tomwalters@0: 
tomwalters@0: 		if( points > segsize )
tomwalters@0: 		    points = segsize ;
tomwalters@0: 
tomwalters@0: 		state->output_size = state->proc( state->states, state->delays, state->input + state->iptr * state->input_size, buffer+bytecount, points, state->active ) ;
tomwalters@0: 	    }
tomwalters@0: 
tomwalters@0: 	    state->iptr += points ;
tomwalters@0: 
tomwalters@0: 	    if( ++state->pointer > 0 )
tomwalters@0: 		bytecount += points*state->chans*state->output_size ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	return ( buffer ) ;
tomwalters@0:     }
tomwalters@0:     else {
tomwalters@0: 	Pull( state->source, 0 ) ;
tomwalters@0: 
tomwalters@0: 	for( last=0 ; last<state->chans ; last++ )
tomwalters@0: 	    Delete( state->states[last] ) ;
tomwalters@0: 
tomwalters@0: 	Delete( state->states ) ;
tomwalters@0: 	Delete( state->delays ) ;
tomwalters@0: 
tomwalters@0: 	return ( DeleteFillableSource( state ) ) ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: Source GenericFilterBank( source, interps, chans, samplerate, center_frequencies, bandwidth_function, scale_function, output_scale, audio_power, order, phase_compensation, input_bits, frames, proc, input_size, info )
tomwalters@0: Source source ;
tomwalters@0: int interps, *chans ;
tomwalters@0: double samplerate, *center_frequencies, (*bandwidth_function)(),  (*scale_function)() ;    /* CG: 23-03-94 */
tomwalters@0: double output_scale, audio_power ;                                                         /* CG: 23-03-94 */
tomwalters@0: int order, phase_compensation, input_bits ;
tomwalters@0: long frames ;
tomwalters@0: int (*proc)(), input_size, info ;                                                          /* CG: 23-03-94 */
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _bank_state *, state ) ;
tomwalters@0:     double sample_delay, max_delay, freq ;
tomwalters@0:     extern Source Interp() ;
tomwalters@0:     Source out ;
tomwalters@0:     int chan ;
tomwalters@0: 
tomwalters@0:     state->chans  = *chans >> interps ;
tomwalters@0: 
tomwalters@0:     state->states = NewArray( RecursiveFilterState *, state->chans, "bank.c for states" ) ;
tomwalters@0:     state->delays = NewArray( int,                    state->chans, "bank.c for delays" ) ;
tomwalters@0: 
tomwalters@0:     state->input_size = input_size ;
tomwalters@0: 
tomwalters@0:     if( info )                                                        /* CG: 23-03-94 */
tomwalters@0:        fprintf( stderr, "\nGTF filterbank information:\n" ) ;
tomwalters@0: 
tomwalters@0:     for( chan=0 ; chan < state->chans ; chan++ ) {
tomwalters@0: 	sample_delay = 0 ;
tomwalters@0: 
tomwalters@0: 	freq = center_frequencies[ chan<<interps ] ;
tomwalters@0:  
tomwalters@0:         if( info )                                                   /* CG: 23-03-94 */
tomwalters@0:            fprintf( stderr, "%3d -- cf:%7.1f Hz =%6.2f ERBs -- bwidth=%5.1f Hz\n",
tomwalters@0: 		   chan + 1 , freq, scale_function( freq ), bandwidth_function( freq ) ) ; 
tomwalters@0: 
tomwalters@0: 	state->states[ chan ] = NewRecursiveFilter( samplerate, freq, bandwidth_function( freq ),
tomwalters@0: 					    output_scale * pow( Audiogram( freq ), audio_power ),
tomwalters@0: 					    order, phase_compensation, input_bits, &sample_delay ) ;
tomwalters@0: 
tomwalters@0: 	if( chan == 0 )
tomwalters@0: 	    max_delay = sample_delay ;
tomwalters@0: 
tomwalters@0: 	if( phase_compensation == FINE_ALIGNMENT || phase_compensation == ENVELOPE_ALIGNMENT )
tomwalters@0: 	    state->delays[ chan ] = sample_delay - max_delay ;
tomwalters@0: 	else
tomwalters@0: 	    state->delays[ chan ] = sample_delay ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     state->output_size = 0 ;
tomwalters@0: 
tomwalters@0:     state->source = NewRetainingSource( source, BUFFSIZE * state->input_size ) ;
tomwalters@0: 
tomwalters@0:     state->input  = ( char * ) 0 ;
tomwalters@0:     state->proc   = proc ;
tomwalters@0: 
tomwalters@0:     out = SetFillableSource( state, bank_callback, "bank.c filter" ) ;
tomwalters@0: 
tomwalters@0:     for( chan = state->chans ; interps-- > 0 ; chan = chan*2-1 )
tomwalters@0: 	out = InterpSource( out, chan ) ;
tomwalters@0: 
tomwalters@0:     *chans = chan ;
tomwalters@0: 
tomwalters@0:     return ( out ) ;
tomwalters@0: }
tomwalters@0: