annotate wdf/bank_tl.c @ 0:5242703e91d3 tip

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

/*
    Copyright (c) Applied Psychology Unit, Medical Research Council. 1994
    ===========================================================================

    Permission to use, copy, modify, and distribute this software without fee 
    is hereby granted for research purposes, provided that this copyright 
    notice appears in all copies and in all supporting documentation, and that 
    the software is not redistributed for any fee (except for a nominal shipping 
    charge). Anyone wanting to incorporate all or part of this software in a
    commercial product must obtain a license from the Medical Research Council.

    The MRC makes no representations about the suitability of this 
    software for any purpose.  It is provided "as is" without express or implied 
    warranty.
 
    THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
    ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
    A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY 
    DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
    AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 
    OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/

/*  
    Acknowledgment:
    ==============

    The source code provided in this file was originally developed by 
    Christian Giguere as part of a Ph.D degree at the Department of
    Engineering of the University of Cambridge from April 1990 to
    November 1993. The code was subsequently adapted under a grant
    from the Hearing Research Trust for full compatibility with 
    AIM Release 6.15.

    Christian Giguere 25/03/94

*/

/*
    ===========================================================
    bank_tl.c
    ===========================================================

    Design of cochlear transmission line filterbank (TLF) with
    coupled outer ear and middle ear (EAR) filter.

    Author        : Christian Giguere
    First written : 01st April, 1991
    Last edited   : 07th March, 1994
    
    Reference:
    (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
    ===========================================================
*/

/***** includes *****/

#include <math.h>
#include <stdio.h>
#include "stitch.h"
#include "source.h"
#include "calc.h"
#include "calc_tl.h"
#include "wdf_tl.h"
#include "formulae.h"
#include "formulae_tl.h"
#include "scales.h"
#include "scales_tl.h"
#include "bank.h"
#include "bank_tl.h"
#include "ear.h"     
#include "wdf_ear.h"
     
/***** defines *****/

#if 0
#define  _DEBUG_ 
#endif

#define   MIN_CF        (     16.0 )     /* Minimum BM segment CF allowed (Hz) */
#define   MAX_CF        (  15000.0 )     /* Maximum BM segment CF allowed (Hz) */
#define   REF_CF        (   1000.0 )     /* Reference CF (Hz) for scaling BM output and Q factor */
#define   L             (      3.5 )     /* Length of BM (cm) */
#define   MinLength     (   1.0E-6 )     /* Mininum length of each BM segment (cm) */
#define   B0            (    0.015 )     /* BM width at basal end (cm) */
#define   B_exp         (     0.30 )     /* Exponential constant for BM width (1/cm) */
#define   A0            (     0.03 )     /* Cross-sectional area of each scala at basal end (cm2) */
#define   A_exp         (    -0.60 )     /* Exponential constant for scalea cross-sectional area (1/cm) */
#define   RHO_fluid     (      1.0 )     /* Density of cochlear fluids [g/cm3] */
#define   MASS_PER_AREA (    0.015 )     /* Transerval mass per area of basilar membrane (g/cm2) */
#define   Q_CHANNEL     (        0 )     /* Introduce channel-dependent Q factor */
#define   RATIO         (     30.0 )     /* Transformer ratio between oval window and eardrum */
#define   WARPING       (        0 )     /* Compensate frequency warping of bilinear transform */


/***** externals */

extern Source InterpSource() ;

/***** functions *****/

/**************************************************************************
* name:                    function:
*
* tlf_bank_callback()      Callable procedure returning pointer to filtered
*                          data (BM velocity or displacement).
*
* TLF_GenBank()            Set-up and initialization function for the design 
*                          of the cochlear filterbank. Returns pointer to
*                          a new source.  
***************************************************************************/

typedef struct _bank_segment  BankSegment ;

struct _bank_segment {
  double  center_frequency ;
  double  seglength, position ;
  int     active ;
  } ;

struct _tlf_bank_state {
  struct           _fillable_source parent ;
  Source           source ; 
  Pointer          input ; 
  int              inout_size ;
  void             (*proc)(), (*closeEar)(), (*closeTLF)() ;
  int              total_chans, display_chans ;
  TLF_BankInfo     *bank ;
  WDFilterState    **states ;
  int              Ntube ;
  WaveWDFstate     *wave_states ;
  EartubeWDFstate  **eartube_states ;
  EarmiddleWDFstate *earmiddle_states ;
  } ;


/************************* tlf_bank_callback() ****************************/

static Pointer tlf_bank_callback( state, bytes, buffer )
struct _tlf_bank_state *state ;
ByteCount *bytes ;
Pointer   buffer ;
{
    register TLF_BankInfo *bankInfo = state->bank ;
    register int          last = *bytes == 0 ;
    register int          points ;
    register ByteCount    bytecount ;

   /***** process *****/
    if( !last ) {

	points = *bytes * bankInfo->decimate_factor / 
		 ( state->display_chans * state->inout_size ) ;

	state->input = Pull( state->source, points * state->inout_size ) ;

	state->proc( bankInfo, state->states, state->input, buffer, points, 
		     state->total_chans, state->wave_states, state->eartube_states,
		     state->earmiddle_states, state->Ntube ) ;

	return ( buffer ) ;
    }

   /***** close *****/
    else {

	Pull( state->source, 0 ) ;
	state->closeEar( state->wave_states, state->eartube_states, state->earmiddle_states, state->Ntube ) ;
	state->closeTLF( state->states, state->total_chans, state->bank ) ;

	return ( DeleteFillableSource( state ) ) ;  
     }
}


/************************** TLF_GenBank() *****************************/

Source TLF_GenBank( source, interps, info, chans, decimate_factor, samplerate, center_frequencies, 
		    output_gain, outdens, qbase, OHC_gain, OHC_sat, motionstr, concha, canal )
Source source ;
int    interps, info, *chans, *decimate_factor ;
double samplerate, *center_frequencies, output_gain, outdens ;
double qbase, OHC_gain, OHC_sat ;
char   *motionstr ;
TubeInfo *concha, *canal ;
{
    DeclareNew( struct _tlf_bank_state *, state ) ;
    BankSegment      **bankSeg, *segmentPointer, **GenerateBankSegments() ;
    double           density, freq, scale_disp, scale_vel ;
    double           qfactor, bn, xn, delta_xn, an, Lt, zov ;
    double           length, diam, attn ;
    Source           outputSource ;
    int              chan, total_chans, segment, counter = 0 ;

   /***** initialise input-related parameters *****/
     state->inout_size = sizeof ( DataType ) ;
     state->input  = ( char * ) 0 ;
     state->source = source ;

   /***** set WDF-TLF filterbank parameters *****/
     density = DensityOnScale( ErbScale( center_frequencies[ 0 ]), ErbScale( center_frequencies[ *chans - 1]),
			       *chans ) ;
     if( density == 0. ) {
	if( outdens <= 0. )
	   density = 1. ;   /* arbitrary */
	else
	   density = outdens ;
     }

     if( interps > 0 )
	interps = MIN( interps, ( int ) floor( log( ( double ) *chans ) / log( 2. ) ) ) ;

     bankSeg = GenerateBankSegments( interps, samplerate, center_frequencies, &density, &outdens, 
				     chans, &total_chans ) ;
     state->display_chans = *chans ;
     state->total_chans  = total_chans ;
     state->states = NewArray( WDFilterState *, state->total_chans, "bank_tl.c for states" ) ;

     if( info )
        fprintf( stderr, "\nTLF filterbank information:\n" ) ; 

     xn = bankSeg[0]->position ;    
     delta_xn = bankSeg[0]->seglength ;    
     Lt = -2. * RHO_fluid / ( A0 * A_exp ) * ( exp( -A_exp * L ) - exp( -A_exp * ( xn + delta_xn ) ) ) ;

     for( chan = 0 ; chan < state->total_chans  ; chan++ ) {
	  segmentPointer = bankSeg[ chan ] ;
	  freq = segmentPointer->center_frequency ;
	  xn = segmentPointer->position ;       
	  delta_xn = MAX( segmentPointer->seglength, MinLength ) ;

	/*** channel-dependent scalea cross-sectional area ***/
	  an = A0 * exp( A_exp * xn ) ;        

	/*** channel-dependent BM width ***/
	  bn = B0 * exp( B_exp * xn ) ;

	/*** channel-dependent Q-factor ***/
	  if ( Q_CHANNEL != 0 )          
	       qfactor = qbase * ( freq / REF_CF ) * ( Erb( REF_CF ) / Erb(freq) ) ;
	  else
	       qfactor = qbase ;

	 /*** print filterbank configuration ***/
          if( info ) {
            counter += 1 * segmentPointer->active ;
	    fprintf( stderr, "%3d -- active=%3d -- cf:%7.1f Hz =%6.2f ERBs -- x=%.3fcm  delta=%.3fcm  b=%.3fcm  A=%.3fcm2\n",
		   state->total_chans - chan, counter * segmentPointer->active, freq, ErbScale( freq ), xn, delta_xn, bn, an ) ; 
	  }

	/*** scale output ***/
	 scale_vel  = output_gain * FILTERBANK_SCALE ;
	 scale_disp = TwoPi * REF_CF * scale_vel ;

	/*** get filter states ***/
	 state->states[ chan ] = WDFilter( samplerate, freq, scale_vel, scale_disp, RHO_fluid, an, bn, 
					   qfactor, MASS_PER_AREA, delta_xn, Lt, WARPING, 
					   segmentPointer->active, &zov, OHC_gain, OHC_sat ) ;

	 Delete( segmentPointer ) ;
     }
   
     Delete( bankSeg ) ;

 
     /*** set parameters common to all channels ***/

     state->bank = New( TLF_BankInfo * ) ;
     state->bank->output_chans = state->display_chans ;

     /* set decimation parameters */
     state->bank->decimateCount = 0 ;
     *decimate_factor = MAX( 1, *decimate_factor ) ;
     state->bank->decimate_factor = *decimate_factor ;    
    

     /***** set output and nonlinearity variables, and filter procedure *****/

     if( strncmp( motionstr, "velocity", 3 ) == 0 ) {
	state->proc = DoWDFdataArray ;
	state->bank->output_var = VELOCITY ;
	state->bank->nl_var = DISPLACEMENT ;
     }
     else {
	state->proc = DoWDFdataArray ;
	state->bank->output_var = DISPLACEMENT ;
	state->bank->nl_var = DISPLACEMENT ;
     } 


   /***** Set WDF-EAR filter design parameters *****/
    state->wave_states = FreefieldWDF( samplerate, RHO_air, C, As, concha->diameter/2. ) ;

    state->Ntube = concha->Nsegments + canal->Nsegments ;
    state->eartube_states = NewArray( EartubeWDFstate *, state->Ntube, "ear.c for eartube states" ) ;

    length = concha->length / concha->Nsegments ;
    diam   = concha->diameter ;
    attn   = concha->att_factor ;
    for( segment = 0 ; segment < concha->Nsegments ; segment++ ) 
	state->eartube_states[ segment ] = EartubeWDF( samplerate, RHO_air, C, diam, length, attn ) ;

    length = canal->length / canal->Nsegments ;
    diam   = canal->diameter ;
    attn   = canal->att_factor ;
    for( ; segment < state->Ntube ; segment++ ) 
	state->eartube_states[ segment ] = EartubeWDF( samplerate, RHO_air, C, diam, length, attn ) ;

    state->earmiddle_states = EarmiddleWDF( samplerate, zov, 1.0, RATIO ) ;


   /***** specify procedures upon closing *****/
    state->closeEar = CloseEarWDF ;
    state->closeTLF = CloseWDF ;

   /***** initialise output-related parameters and return *****/
    outputSource = SetFillableSource( state, tlf_bank_callback, "bank_tl.c filter" ) ;

    for( *chans = state->display_chans ; interps-- > 0 ; *chans = *chans * 2 - 1, density = density * 2. )
	outputSource = InterpSource( outputSource, *chans ) ;

    return ( outputSource ) ;
}


/************************** lower level functions *************************/

BankSegment **GenerateBankSegments( interps, samplerate, display_frequencies, display_dens, 
				    out_dens, display_chans, total_chans )
int    interps ;
double samplerate ;
double *display_frequencies ;
double *display_dens, *out_dens ;
int    *display_chans, *total_chans ;
{
    BankSegment **bankSeg, *segmentPointer ;
    double      min_cf, max_cf ;
    double      apical_cf, basal_cf ;
    double      *apical_frequencies, *basal_frequencies ;
    double      apical_dens, basal_dens, central_dens ;
    int         apical_chans, basal_chans, central_chans ;
    int         ichan, chan ;
    int         remainder ; 
    double      cm_per_CB = GetERBscaling ( ) / 10. ;


   /***** set number and density of central channels *****/
    if( interps >= 0 ) {

	central_chans = *display_chans = *display_chans >> interps ;
	central_dens  = *display_dens  = ldexp( *display_dens, -interps ) ;
	min_cf        = FofErbScale( ErbScale( display_frequencies[0] ) - 1. / central_dens ) ;      
	max_cf        = display_frequencies[ ( central_chans - 1 ) << interps ] ;
    }

    else {
	central_chans = *display_chans << -interps ;
	central_dens  = ldexp( *display_dens, -interps ) ;
	min_cf        = FofErbScale( ErbScale( display_frequencies[0] ) 
				     - ( 1 << -interps ) / central_dens ) ;
	max_cf        = display_frequencies[*display_chans - 1] ; 
    }


   /***** set number and density of apical and basal channels outside display range *****/ 

    if( *out_dens <= 0. ) {

	apical_cf = min_cf ;
	apical_dens  = 0. ;
	basal_cf = max_cf ;
	basal_dens = 0. ;
    }
    else {
	apical_cf = adjustCF( MIN_CF, samplerate ) ;
	apical_dens = *out_dens ;
	basal_cf  = adjustCF( MAX_CF, samplerate ) ;
	basal_dens = *out_dens ;
    }

	apical_chans = 0 ;
	basal_chans = 0 ;
    

   /***** generate apical and basal center frequencies outside display range *****/ 

    apical_frequencies = GenerateFrequencies( apical_cf, min_cf, min_cf, &apical_dens, &apical_chans ) ;
    basal_frequencies = GenerateFrequencies( max_cf, basal_cf, max_cf, &basal_dens, &basal_chans ) ;

   /***** fill in array of BanKSegment structures with BM segment info *****/
    *total_chans = apical_chans + central_chans + basal_chans - 2 ;   
    bankSeg = NewArray( BankSegment *, *total_chans + 1, "tl_bank.c for segments" ) ;

    ichan = 0 ;
    for( chan = 1 ; chan < apical_chans ; chan++ ) {

	 bankSeg[ ichan++ ] = segmentPointer = New( BankSegment * ) ; 
	 segmentPointer->center_frequency = apical_frequencies[ chan ] ;
	 segmentPointer->seglength = cm_per_CB / apical_dens ;  
	 segmentPointer->position = L - cm_per_CB * ErbScale( segmentPointer->center_frequency ) ; 
	 segmentPointer->active = 0 ;
    }

    for( chan = 0 ; chan < central_chans ; chan++ ) {

	 bankSeg[ ichan++ ] = segmentPointer = New( BankSegment * ) ;         

	 if( interps >= 0 ) {
	     segmentPointer->center_frequency = display_frequencies[ chan << interps ] ;
	     segmentPointer->active = 1 ;
	 }
	 else {
	     remainder =  chan % ( 1 << -interps ) ;
	     segmentPointer->center_frequency = FofErbScale( ErbScale( display_frequencies[ chan >> -interps ] )
						+ ( remainder - ( 1 << -interps ) + 1 ) / central_dens ) ;
	     segmentPointer->active = ( remainder == ( ( 1 << -interps ) - 1 ) ) ;
	 }

	 segmentPointer->seglength = cm_per_CB / central_dens ;  
	 segmentPointer->position = L - cm_per_CB * ErbScale( segmentPointer->center_frequency ) ; 
    }


    for( chan = 1 ; chan < basal_chans ;  chan++ ) {

	 bankSeg[ ichan++ ] = segmentPointer = New( BankSegment * ) ;
	 segmentPointer->center_frequency = basal_frequencies[ chan ] ;
	 segmentPointer->seglength = cm_per_CB / basal_dens ;  
	 segmentPointer->position = L - cm_per_CB * ErbScale( segmentPointer->center_frequency ) ; 
	 segmentPointer->active = 0 ;
    }

    segmentPointer = bankSeg[ 0 ] ;
    if( ( segmentPointer->position + segmentPointer->seglength ) > L ) 
	segmentPointer->seglength = L - MIN( L, segmentPointer->position ) ;

   /***** return and deallocate dynamic memory *****/
    Delete( apical_frequencies ) ;
    Delete(  basal_frequencies ) ;

    return( bankSeg ) ;
}


double *GenerateFrequencies( min_cf, max_cf, base_cf, density, channels )
double  min_cf, max_cf, base_cf ;
double  *density ;
int     *channels ;
{
    double freq, *frequencies ;

   /*** map characteristic frequencies onto specified scale ***/ 
    min_cf = ErbScale( min_cf ) ;
    max_cf = ErbScale( max_cf ) ;
    max_cf = MAX( min_cf, max_cf ) ;                 /* max_cf cannot be smaller than min_cf */ 
    base_cf = ErbScale( base_cf ) ;


   /*** call appropriate generating functions ***/
    if( *channels <= 0 ) {

       if( *density <= 0. ) {
	  *channels = 1 ;                           /* there must be at least one channel */
	  freq = FofErbScale( min_cf ) ;            /* arbitrary */
	  frequencies = &freq ;
	  *density = 1. ;                           /* arbitrary */

       }
       else {
	  frequencies = GenerateScale( min_cf, max_cf, *density, base_cf, FofErbScale ) ;
	  *channels   = NumberOnScale( min_cf, max_cf, *density, base_cf ) ;
       }
    }

    else {
       frequencies = NumberedScale( min_cf, max_cf, *channels, FofErbScale ) ;
       *density    = DensityOnScale( min_cf, max_cf, *channels ) ;
       if( *density == 0. )
	  *density = 1. ;    /* arbitrary */
    }

   /*** return pointer to array of center frequencies ***/
    return ( frequencies ) ;
}


double adjustCF( cf, samplerate )
double cf, samplerate ;
{
    cf = MAX( cf, MIN_CF ) ;                      
    cf = MAX( cf, FofErbScale( 0. ) ) ;                    

    cf = MIN( cf, samplerate / 2. ) ;
    cf = MIN( cf, MAX_CF ) ;
    cf = MIN( cf, FofErbScale( 10. * L / GetERBscaling() ) );

    return ( cf ) ;
}