diff filter/recurse.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/filter/recurse.c	Fri May 20 15:19:45 2011 +0100
@@ -0,0 +1,235 @@
+/*
+    Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
+    ===========================================================================
+
+    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.
+*/
+
+/*
+             ================================================
+              recurse.c - start-up code for recursive filter
+             ================================================
+
+
+    Copyright (c), 1989  The Medical Research Council, Applied Psychology Unit.
+
+
+    Author  : John Holdsworth
+    Written : 22th March, 1989.
+
+    Edited  :
+
+*/
+
+
+#include <math.h>
+
+#ifndef  _STITCH_H_
+#include "stitch.h"
+#endif
+
+#ifndef  _RECURSE_H_
+#include "recurse.h"
+#endif
+
+#define   Pi                 ( 3.1415926535 )
+#define   TwoPi              ( 2*Pi )
+
+#ifndef  lint
+static char *sccs_id = "@(#)recurse.c	1.8 John Holdsworth (MRC-APU) 11/8/90" ;
+#endif
+
+static void generateSineTable() ;
+static double factorial() ;
+
+#ifdef DSP32
+int filterSineTableBits  =  8 ;
+int filterSineTableShift =  0 ;
+#else
+int filterSineTableBits  = 12 ;
+int filterSineTableShift = 14 ;
+#endif
+unsigned long filterSineTableSize, filterSineTableMask ;
+Table filterSineTable, filterCosineTable ;
+
+char *bstart, *bend ;
+
+
+RecursiveFilterState *NewRecursiveFilter( samplerate, center_frequency, bandwidth, output_scale, order, phase_compensation, input_bits, sample_delay )
+double samplerate, center_frequency, bandwidth, output_scale ;
+int order, phase_compensation, input_bits ;
+double *sample_delay ;
+{
+    DeclareNew( RecursiveFilterState *, filter_state ) ;
+    double omega ;
+
+    if( filterSineTableSize != 1 << filterSineTableBits )
+	generateSineTable() ;
+
+    filter_state->order        = order        ;
+    filter_state->input_bits   = input_bits   ;
+
+    /* compensation for nominal 6dB loss in filter */
+
+    filter_state->output_scale = output_scale * 2. ;
+
+    omega = TwoPi * bandwidth / bandwidth_normalisation( filter_state->order ) ;
+
+    /* delta phi of carrier between sample ticks */
+
+    filter_state->output_delta_phi = center_frequency / samplerate * ( filterSineTableSize << 16 ) + 0.5 ;
+    filter_state->delta_phi        = filter_state->output_delta_phi ;
+
+    /* initialise phases to give rounding */
+
+    ALL( filter_state->phi ) = ALL( filter_state->output_phi ) = 1l << 15 ;
+
+    /* calculate time delay required for phase compensation selected */
+    /* time advance is currently implemented as delay of       */
+    /* maximum advance minus desired advance                   */
+
+    if( center_frequency < samplerate / 4 )
+	filter_state->over_sample = 1 ;
+    else {
+	/* if center frequency is over half the nyquist rate filtering process */
+	/* is performed at twice the sampling rate to avoid aliases */
+
+	filter_state->over_sample = 2 ;
+
+	ALL( filter_state->phi ) -= filter_state->delta_phi / 4 ;
+#ifndef DSP32
+	ALL( filter_state->phi ) &= filterSineTableMask ;
+#endif
+	filter_state->delta_phi /= 2 ;
+
+	/* must amplitude compensate for attenuation of high frequencies introduced by over sampling */
+
+	filter_state->output_scale /= ( filter_state->over_sample * cos( TwoPi * center_frequency / samplerate / 4 ) ) ;
+    }
+
+    filter_state->k = ( double ) 1. - exp( -omega / samplerate / filter_state->over_sample ) ;
+
+    /* fiddle delay time to take into acount half sample  */
+    /* interval advance introduced by difference equation */
+
+#if 0 /* not needed any more as alternate forward and backward differences */
+
+    *sample_delay += filter_state->order / filter_state->over_sample / 2 ;
+
+    /* compensate carrier phase for delay introduced above */
+
+    ALL( filter_state->output_phi ) += filter_state->order / filter_state->over_sample / 2 *
+			    filter_state->delta_phi * filter_state->over_sample ;
+#ifndef DSP32
+    ALL( filter_state->output_phi ) &= filterSineTableMask ;
+#endif
+#endif
+
+    /* munge to sin phase gammatone for mfsais */
+
+    ALL( filter_state->output_phi ) += 3 * filterSineTableMask / 4  ;
+
+    /* perform required phase compensation */
+
+    if( phase_compensation > 0 )
+	*sample_delay -= phase_compensation / center_frequency * samplerate ;
+    else if( phase_compensation < 0 )
+	*sample_delay -= ( filter_state->order - 1. ) / omega * samplerate ;
+
+    /* if phase compensation of type -2 or -4 is asked for the carrier phase */
+    /* is shifted to be aligned with the envelope maxima */
+
+    if( phase_compensation == FINE_ALIGNMENT || phase_compensation <= ACUASAL + FINE_ALIGNMENT ) {
+	ALL( filter_state->phi ) += ( filter_state->order - 1. ) / omega * center_frequency * ( filterSineTableSize << 16 ) + 0.5 ;
+#ifndef DSP32
+	ALL( filter_state->phi ) &= filterSineTableMask ;
+#endif
+    }
+
+    /* new recurse filter coefts */
+
+    filter_state->gain = pow(  2 * omega / samplerate, (double) order ) * output_scale ;
+    filter_state->k1   = exp(     -omega / samplerate ) * cos( TwoPi * center_frequency / samplerate ) ;
+    filter_state->k2   = exp( -2 * omega / samplerate ) ;
+
+    /* flag channel not fully initialised */
+    /* leave filter specific initialisation to specific filter code */
+
+    filter_state->states = ( char * ) 0 ;
+
+    /* generate sin lookup table if necessary */
+
+    return ( filter_state ) ;
+}
+
+double bandwidth_normalisation( order )
+int order ;
+{
+    return( Pi * factorial( 2*order - 2 ) / factorial( order - 1 ) / factorial( order - 1 ) / ( 1 << ( 2 * order - 2 ) ) ) ;
+}
+
+static double factorial( n )
+int n ;
+{
+    double fact ;
+    int i ;
+
+    fact = 1. ;
+
+    for( i = n ; i > 1 ; i-- )
+	fact *= i ;
+
+    return( fact ) ;
+}
+
+
+static void generateSineTable()
+{
+    static double filterTwoPi ;
+    register int i ;
+
+    if( filterSineTableSize != 1 << filterSineTableBits )
+    {
+	if( filterSineTable != ( Table ) 0 )
+	    stitch_free( ( Pointer ) filterSineTable ) ;
+
+	if( filterTwoPi == 0. )
+	    filterTwoPi = atan( 1. ) * 8. ;
+
+	filterSineTableSize = 1l << filterSineTableBits ;
+	filterSineTableMask = ( filterSineTableSize << 16 ) - 1 ;
+
+	/* sin table 5/4 times larger than required for cosine values */
+
+	filterSineTable = ( Table ) stitch_malloc( ( unsigned ) ( filterSineTableSize * 5 / 4 + 1 ) * sizeof( *filterSineTable ), "recurse.c for log table" ) ;
+
+	for( i=0 ; i <= filterSineTableSize/4 ; i++ ) {
+	    filterSineTable[ i ] = sin( filterTwoPi * i / filterSineTableSize ) * ( 1<<filterSineTableShift ) ;
+	    filterSineTable[ filterSineTableSize/2 - i ] =  filterSineTable[ i ] ;
+	    filterSineTable[ filterSineTableSize/2 + i ] = -filterSineTable[ i ] ;
+	    filterSineTable[ filterSineTableSize   - i ] = -filterSineTable[ i ] ;
+	    filterSineTable[ filterSineTableSize   + i ] =  filterSineTable[ i ] ;
+	}
+
+	filterCosineTable = filterSineTable + filterSineTableSize / 4 ;
+    }
+
+    return ;
+}
+