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:     Acknowledgment:
tomwalters@0:     ==============
tomwalters@0: 
tomwalters@0:     The source code provided in this file was originally developed by 
tomwalters@0:     Christian Giguere as part of a Ph.D degree at the Department of
tomwalters@0:     Engineering of the University of Cambridge from April 1990 to
tomwalters@0:     November 1993. The code was subsequently adapted under a grant
tomwalters@0:     from the Hearing Research Trust for full compatibility with 
tomwalters@0:     AIM Release 6.15.
tomwalters@0: 
tomwalters@0:     Christian Giguere 25/03/94
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:     ===========================================================
tomwalters@0:     fir.c
tomwalters@0:     ===========================================================
tomwalters@0: 
tomwalters@0:     Finite Impulse Response (FIR) filtering module.
tomwalters@0:     
tomwalters@0:     Author        : Christian Giguere
tomwalters@0:     First written : 12th November, 1990
tomwalters@0:     Last edited   : 22nd September, 1991
tomwalters@0:     
tomwalters@0:     Reference:
tomwalters@0:     (1) A.V.Oppenheim and R.W.Schafer (1975). Digital Signal
tomwalters@0:         Processing (Prentice-Hall), Sections 4.5.5 and 5.5.
tomwalters@0:     ===========================================================
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** includes *****/
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: #include <stdio.h>
tomwalters@0: #include "stitch.h"
tomwalters@0: #include "calc.h"
tomwalters@0: #include "calc_tl.h"
tomwalters@0: #include "fir.h"
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** defines *****/
tomwalters@0: 
tomwalters@0: #define NUMBER_OF_COEFFS     (               51 )                 /* odd order */
tomwalters@0: #define NUMBER_OF_STATES     ( NUMBER_OF_COEFFS ) 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** functions *****/
tomwalters@0: 
tomwalters@0: /*******************************************************************************
tomwalters@0: * name:                 function:  
tomwalters@0: *
tomwalters@0: * NewLpFIRfilter()      Design of LowPass FIR filter using windowing technique
tomwalters@0: *                       (Hamming window). It generates the multiplier coefficients
tomwalters@0: *                       and initializes the states of the filter. It returns a 
tomwalters@0: *                       pointer to a structure containing all relevant information 
tomwalters@0: *                       for the realisation of the filter.
tomwalters@0: *
tomwalters@0: * DoFIRfilter()         Direct form realisation for FIR filter of odd order linear
tomwalters@0: *                       phase. It is called by function ``upsample_callback()'' 
tomwalters@0: *                       in file ``upsample.c'' once for a specified number of input 
tomwalters@0: *                       points.  It computes the filter output for each input point
tomwalters@0: *                       and keeps track of the filter states. 
tomwalters@0: *
tomwalters@0: * CloseFIRfilter()      Deletes all private data structures and arrays of the FIR 
tomwalters@0: *                       filter upon comletion of filtering.  It is called by function
tomwalters@0: *                       ``upsample_callback()'' in file ``upsample.c''. 
tomwalters@0: ********************************************************************************/
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /************************* NewLpFIRfilter() *******************************/
tomwalters@0: 
tomwalters@0: FIRfilterState *NewLpFIRfilter( wc, upsample_factor, delay )
tomwalters@0:   double wc ;
tomwalters@0:   int    upsample_factor, *delay ;
tomwalters@0: {
tomwalters@0:   DeclareNew( FIRfilterState *, filter_state ) ;
tomwalters@0:   CoeffType *hn ;
tomwalters@0:   double    filterGain, hamming() ;
tomwalters@0:   int       n ;
tomwalters@0: 
tomwalters@0:  /***** generate and store FIR filter coefficients (linear phase, odd order, window desing technique) *****/
tomwalters@0:   filter_state->number_of_coeffs = ( NUMBER_OF_COEFFS + 1 ) / 2 ;  
tomwalters@0:   filter_state->coeffs = ( char * ) NewArray( CoeffType, filter_state->number_of_coeffs, "fir.c for coefficients\n" ) ;
tomwalters@0: 
tomwalters@0:   *delay = ( NUMBER_OF_COEFFS - 1 ) / 2 ;
tomwalters@0:   hn = ( CoeffType * ) filter_state->coeffs ;
tomwalters@0:   filterGain = ( double ) upsample_factor ;
tomwalters@0: 
tomwalters@0:   for( n = 0 ; n < filter_state->number_of_coeffs - 1 ; n++ ) 
tomwalters@0:       *hn++ = sin( wc * ( n - *delay ) ) / ( Pi * ( n - *delay ) ) * hamming( n, NUMBER_OF_COEFFS ) * filterGain ;
tomwalters@0:   *hn = wc / Pi * filterGain ;
tomwalters@0:     
tomwalters@0: 
tomwalters@0:  /***** initialise filter states *****/
tomwalters@0:   filter_state->number_of_states = NUMBER_OF_STATES ;
tomwalters@0:   filter_state->states = ( char * ) NewZeroedArray( StateType, NUMBER_OF_STATES, "fir.c for filter states\n" ) ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0:  /***** return *****/   
tomwalters@0:   return( filter_state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /****************************** DoFIRfilter() *******************************/
tomwalters@0: 
tomwalters@0: void DoFIRfilter( filter_state, inoutput, points )
tomwalters@0:   FIRfilterState *filter_state ;
tomwalters@0:   DataType       *inoutput ;
tomwalters@0:   int            points ;
tomwalters@0: {
tomwalters@0:   register StateType  *st = ( StateType * ) filter_state->states ;
tomwalters@0:   register CoeffType  *cf = ( CoeffType * ) filter_state->coeffs ;
tomwalters@0:   register DataType   *inout_ptr = inoutput ;
tomwalters@0:   register double     acc ;
tomwalters@0:   register int        n, Ns = filter_state->number_of_states ;
tomwalters@0:   register int        Nc = filter_state->number_of_coeffs ;
tomwalters@0: 
tomwalters@0:  /***** for all points *****/
tomwalters@0:   for( ; points > 0 ; points-- ) {
tomwalters@0:       acc = 0. ;
tomwalters@0: 
tomwalters@0:      /*** shift the states ***/
tomwalters@0:       for( n = Ns - 1 ; n > 0 ; n-- ) 
tomwalters@0:           st[n] = st[n-1] ;
tomwalters@0:       st[0] = *inout_ptr ;
tomwalters@0:    
tomwalters@0:      /*** filter current point ***/
tomwalters@0:       for( n = 0  ; n < Nc - 1 ; n++ ) 
tomwalters@0:           acc += cf[n] * ( st[n] + st[Ns-1-n] ) ; 
tomwalters@0:       acc += cf[n] * st[n] ;
tomwalters@0: 
tomwalters@0:      /*** store output point ***/
tomwalters@0:       *inout_ptr++ = acc ;
tomwalters@0:   }  
tomwalters@0: 
tomwalters@0:  /***** return *****/
tomwalters@0:   return ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /************************* CloseFIRfilter() *******************************/
tomwalters@0: 
tomwalters@0: void CloseFIRfilter( filter_state )
tomwalters@0:   FIRfilterState *filter_state ;
tomwalters@0: {
tomwalters@0:   Delete( filter_state->states ) ;
tomwalters@0:   Delete( filter_state->coeffs ) ;
tomwalters@0:   Delete( filter_state ) ;
tomwalters@0: 
tomwalters@0:   return ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*********************** low level functions ***********************/
tomwalters@0: 
tomwalters@0: double hamming( n, N )
tomwalters@0:   int n, N ;
tomwalters@0: { 
tomwalters@0:   return ( 0.54 - 0.46 * cos( TwoPi * n / ( N - 1. ) ) ) ;
tomwalters@0: }
tomwalters@0: