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: THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING tomwalters@0: ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE tomwalters@0: A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY tomwalters@0: DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN tomwalters@0: AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF tomwalters@0: OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. tomwalters@0: */ tomwalters@0: tomwalters@0: /* tomwalters@0: ================================================ tomwalters@0: recurse.c - start-up code for recursive filter tomwalters@0: ================================================ 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 : 22th March, 1989. tomwalters@0: tomwalters@0: Edited : tomwalters@0: tomwalters@0: */ tomwalters@0: tomwalters@0: tomwalters@0: #include tomwalters@0: tomwalters@0: #ifndef _STITCH_H_ tomwalters@0: #include "stitch.h" tomwalters@0: #endif tomwalters@0: tomwalters@0: #ifndef _RECURSE_H_ tomwalters@0: #include "recurse.h" tomwalters@0: #endif tomwalters@0: tomwalters@0: #define Pi ( 3.1415926535 ) tomwalters@0: #define TwoPi ( 2*Pi ) tomwalters@0: tomwalters@0: #ifndef lint tomwalters@0: static char *sccs_id = "@(#)recurse.c 1.8 John Holdsworth (MRC-APU) 11/8/90" ; tomwalters@0: #endif tomwalters@0: tomwalters@0: static void generateSineTable() ; tomwalters@0: static double factorial() ; tomwalters@0: tomwalters@0: #ifdef DSP32 tomwalters@0: int filterSineTableBits = 8 ; tomwalters@0: int filterSineTableShift = 0 ; tomwalters@0: #else tomwalters@0: int filterSineTableBits = 12 ; tomwalters@0: int filterSineTableShift = 14 ; tomwalters@0: #endif tomwalters@0: unsigned long filterSineTableSize, filterSineTableMask ; tomwalters@0: Table filterSineTable, filterCosineTable ; tomwalters@0: tomwalters@0: char *bstart, *bend ; tomwalters@0: tomwalters@0: tomwalters@0: RecursiveFilterState *NewRecursiveFilter( samplerate, center_frequency, bandwidth, output_scale, order, phase_compensation, input_bits, sample_delay ) tomwalters@0: double samplerate, center_frequency, bandwidth, output_scale ; tomwalters@0: int order, phase_compensation, input_bits ; tomwalters@0: double *sample_delay ; tomwalters@0: { tomwalters@0: DeclareNew( RecursiveFilterState *, filter_state ) ; tomwalters@0: double omega ; tomwalters@0: tomwalters@0: if( filterSineTableSize != 1 << filterSineTableBits ) tomwalters@0: generateSineTable() ; tomwalters@0: tomwalters@0: filter_state->order = order ; tomwalters@0: filter_state->input_bits = input_bits ; tomwalters@0: tomwalters@0: /* compensation for nominal 6dB loss in filter */ tomwalters@0: tomwalters@0: filter_state->output_scale = output_scale * 2. ; tomwalters@0: tomwalters@0: omega = TwoPi * bandwidth / bandwidth_normalisation( filter_state->order ) ; tomwalters@0: tomwalters@0: /* delta phi of carrier between sample ticks */ tomwalters@0: tomwalters@0: filter_state->output_delta_phi = center_frequency / samplerate * ( filterSineTableSize << 16 ) + 0.5 ; tomwalters@0: filter_state->delta_phi = filter_state->output_delta_phi ; tomwalters@0: tomwalters@0: /* initialise phases to give rounding */ tomwalters@0: tomwalters@0: ALL( filter_state->phi ) = ALL( filter_state->output_phi ) = 1l << 15 ; tomwalters@0: tomwalters@0: /* calculate time delay required for phase compensation selected */ tomwalters@0: /* time advance is currently implemented as delay of */ tomwalters@0: /* maximum advance minus desired advance */ tomwalters@0: tomwalters@0: if( center_frequency < samplerate / 4 ) tomwalters@0: filter_state->over_sample = 1 ; tomwalters@0: else { tomwalters@0: /* if center frequency is over half the nyquist rate filtering process */ tomwalters@0: /* is performed at twice the sampling rate to avoid aliases */ tomwalters@0: tomwalters@0: filter_state->over_sample = 2 ; tomwalters@0: tomwalters@0: ALL( filter_state->phi ) -= filter_state->delta_phi / 4 ; tomwalters@0: #ifndef DSP32 tomwalters@0: ALL( filter_state->phi ) &= filterSineTableMask ; tomwalters@0: #endif tomwalters@0: filter_state->delta_phi /= 2 ; tomwalters@0: tomwalters@0: /* must amplitude compensate for attenuation of high frequencies introduced by over sampling */ tomwalters@0: tomwalters@0: filter_state->output_scale /= ( filter_state->over_sample * cos( TwoPi * center_frequency / samplerate / 4 ) ) ; tomwalters@0: } tomwalters@0: tomwalters@0: filter_state->k = ( double ) 1. - exp( -omega / samplerate / filter_state->over_sample ) ; tomwalters@0: tomwalters@0: /* fiddle delay time to take into acount half sample */ tomwalters@0: /* interval advance introduced by difference equation */ tomwalters@0: tomwalters@0: #if 0 /* not needed any more as alternate forward and backward differences */ tomwalters@0: tomwalters@0: *sample_delay += filter_state->order / filter_state->over_sample / 2 ; tomwalters@0: tomwalters@0: /* compensate carrier phase for delay introduced above */ tomwalters@0: tomwalters@0: ALL( filter_state->output_phi ) += filter_state->order / filter_state->over_sample / 2 * tomwalters@0: filter_state->delta_phi * filter_state->over_sample ; tomwalters@0: #ifndef DSP32 tomwalters@0: ALL( filter_state->output_phi ) &= filterSineTableMask ; tomwalters@0: #endif tomwalters@0: #endif tomwalters@0: tomwalters@0: /* munge to sin phase gammatone for mfsais */ tomwalters@0: tomwalters@0: ALL( filter_state->output_phi ) += 3 * filterSineTableMask / 4 ; tomwalters@0: tomwalters@0: /* perform required phase compensation */ tomwalters@0: tomwalters@0: if( phase_compensation > 0 ) tomwalters@0: *sample_delay -= phase_compensation / center_frequency * samplerate ; tomwalters@0: else if( phase_compensation < 0 ) tomwalters@0: *sample_delay -= ( filter_state->order - 1. ) / omega * samplerate ; tomwalters@0: tomwalters@0: /* if phase compensation of type -2 or -4 is asked for the carrier phase */ tomwalters@0: /* is shifted to be aligned with the envelope maxima */ tomwalters@0: tomwalters@0: if( phase_compensation == FINE_ALIGNMENT || phase_compensation <= ACUASAL + FINE_ALIGNMENT ) { tomwalters@0: ALL( filter_state->phi ) += ( filter_state->order - 1. ) / omega * center_frequency * ( filterSineTableSize << 16 ) + 0.5 ; tomwalters@0: #ifndef DSP32 tomwalters@0: ALL( filter_state->phi ) &= filterSineTableMask ; tomwalters@0: #endif tomwalters@0: } tomwalters@0: tomwalters@0: /* new recurse filter coefts */ tomwalters@0: tomwalters@0: filter_state->gain = pow( 2 * omega / samplerate, (double) order ) * output_scale ; tomwalters@0: filter_state->k1 = exp( -omega / samplerate ) * cos( TwoPi * center_frequency / samplerate ) ; tomwalters@0: filter_state->k2 = exp( -2 * omega / samplerate ) ; tomwalters@0: tomwalters@0: /* flag channel not fully initialised */ tomwalters@0: /* leave filter specific initialisation to specific filter code */ tomwalters@0: tomwalters@0: filter_state->states = ( char * ) 0 ; tomwalters@0: tomwalters@0: /* generate sin lookup table if necessary */ tomwalters@0: tomwalters@0: return ( filter_state ) ; tomwalters@0: } tomwalters@0: tomwalters@0: double bandwidth_normalisation( order ) tomwalters@0: int order ; tomwalters@0: { tomwalters@0: return( Pi * factorial( 2*order - 2 ) / factorial( order - 1 ) / factorial( order - 1 ) / ( 1 << ( 2 * order - 2 ) ) ) ; tomwalters@0: } tomwalters@0: tomwalters@0: static double factorial( n ) tomwalters@0: int n ; tomwalters@0: { tomwalters@0: double fact ; tomwalters@0: int i ; tomwalters@0: tomwalters@0: fact = 1. ; tomwalters@0: tomwalters@0: for( i = n ; i > 1 ; i-- ) tomwalters@0: fact *= i ; tomwalters@0: tomwalters@0: return( fact ) ; tomwalters@0: } tomwalters@0: tomwalters@0: tomwalters@0: static void generateSineTable() tomwalters@0: { tomwalters@0: static double filterTwoPi ; tomwalters@0: register int i ; tomwalters@0: tomwalters@0: if( filterSineTableSize != 1 << filterSineTableBits ) tomwalters@0: { tomwalters@0: if( filterSineTable != ( Table ) 0 ) tomwalters@0: stitch_free( ( Pointer ) filterSineTable ) ; tomwalters@0: tomwalters@0: if( filterTwoPi == 0. ) tomwalters@0: filterTwoPi = atan( 1. ) * 8. ; tomwalters@0: tomwalters@0: filterSineTableSize = 1l << filterSineTableBits ; tomwalters@0: filterSineTableMask = ( filterSineTableSize << 16 ) - 1 ; tomwalters@0: tomwalters@0: /* sin table 5/4 times larger than required for cosine values */ tomwalters@0: tomwalters@0: filterSineTable = ( Table ) stitch_malloc( ( unsigned ) ( filterSineTableSize * 5 / 4 + 1 ) * sizeof( *filterSineTable ), "recurse.c for log table" ) ; tomwalters@0: tomwalters@0: for( i=0 ; i <= filterSineTableSize/4 ; i++ ) { tomwalters@0: filterSineTable[ i ] = sin( filterTwoPi * i / filterSineTableSize ) * ( 1<