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

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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/wdf/bank_tl.c	Fri May 20 15:19:45 2011 +0100
@@ -0,0 +1,499 @@
+/*
+    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 ) ;
+}
+