view filter/recurse.c @ 0:5242703e91d3 tip

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

/*
    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 ;
}