/* 
 *  =====================================================================
 *  test_convolution.c
 * 
 *      - implementation of the convolution sum of two arrays
 *      - compute the convolution of input x with the impulse response h
 * 
 *                m=n
 *         y(n) = SUM [ x(m).h(n-m) ]
 *                n=0
 * 
 *  Output : Columns 1,2, and 3 are the input function, impulse response
 *           function, and computed convolution, respectively.
 * 
 *  =====================================================================
 * 
 *  Copyright (C) 1994-96 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *  =====================================================================
 */

int main( void ) {
int iCount, iCount2;

  /* [a] : Define Input Function and Lowpass Impulse response */

  float faInput[] =  {0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0};
  float faImpRes[]=  {1.00000,  0.36000, 0.14000, 0.05000, 0.01800,
                      0.00670,  0.00250, 0.00090, 0.00034, 0.00012};  

  /* [b] : compute size of working parameters and arrays */

  float faOutput[(sizeof(faInput)+sizeof(faImpRes))/sizeof(faImpRes[0])-1];

  int iOutLength = sizeof(faOutput)/sizeof(faOutput[0]);
  int iInpLength = sizeof(faInput)/sizeof(faInput[0]);
  int iResLength = sizeof(faImpRes)/sizeof(faImpRes[0]);

  /* [c] : Initialize the output array */

  for (iCount = 0; iCount < iOutLength; iCount++)
       faOutput[ iCount ] = 0;

  /* [d] : Compute Convolution with bounds checking on arrays */

  for ( iCount = 0; iCount  < iOutLength; iCount++ ) {
     for (iCount2 = 0; iCount2 < iResLength; iCount2++ ) {
         if (iCount - iCount2 >= 0 &&
            iCount - iCount2 < iInpLength )
            faOutput[iCount] += faImpRes[iCount2]*faInput[iCount-iCount2];
     }
  }

  /* [c] : write output */

  for ( iCount = 0; iCount < iOutLength; iCount++ )
      if ( iCount >= iInpLength || iCount >= iResLength ){
        if ( iCount >= iInpLength && iCount >= iResLength )
             printf( "NaN    NaN    %6.4f \n", faOutput[iCount] );
        if ( iCount >= iInpLength && iCount < iResLength )
             printf( "NaN    %6.4f  %6.4f \n", faImpRes[iCount], faOutput[iCount] );
        if ( iCount < iInpLength && iCount >= iResLength )
             printf( "%6.4f NaN    %6.4f \n", faInput[iCount], faOutput[iCount] );
      } else
             printf( "%6.4f %6.4f %6.4f \n", faInput[iCount], faImpRes[iCount], faOutput[iCount] );

}