convsinc.c (corrected)

/*
 convsinc.c

 Simple demo of a sinc-based convolution used to downsample a data vector to
 demonstrate the absence of "ringing" when using sinc filtering to downsample data.
 See e.g. http://en.wikipedia.org/wiki/Sinc_filter for some info on sinc convolution.

 Written, and put in the public domain by Antisthenes@[DPreview forums] 2014/09/21
 Related link: http://www.dpreview.com/forums/post/54425513

 On a Unix-like system, compile with:
   cc -lm -o convsinc convsinc.c

 Run:
   ./convsinc

 Output:
   The program will then blurp out the numerical values of the downsampled data.
   I'm lazy, so there are no fancy graphs or anything.
 */


#define ORIGINAL_WIDTH      55
#define DOWNSAMPLED_WIDTH   29

// Number of "half-periods" of the sin(x)/x function we'll use in the convolution.
// Increasing this number increases the computation time.  7 already provides a
// good accuracy.
#define CONVOLUTION_SPAN    7

// Define this symbol as "step", "pulse", "square" or "sine" to specify the type
// of function with which the original data should be initialized with.
#define INITIALIZE_WITH     step


// In principle, nothing needs to be changed / configured in the source code
// from below this comment to run this trivial demo...

// step:
//       _____
//      |
// _____|

// pulse:
//
//      |
//______|______

// square:
//       _____       _____       _____
//      |     |     |     |     |     |
// _____|     |_____|     |_____|     |

// sine:
//   _       _       _
//  / \     / \     / \
// /   \   /   \   /   \   /
//      \_/     \_/     \_/

enum F_TYPE {
    step,
    pulse,
    square,
    sine
};

#define _USE_MATH_DEFINES

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, const char * argv[])
{
    int
        i, j, s, originalWidth, downsampledWidth;
    double
        *originalData, *downsampledData, *ptr1, *ptr2, accumulator, scaleFactor, x, d;
    enum F_TYPE
        initializationType = INITIALIZE_WITH;

    originalWidth = ORIGINAL_WIDTH;
    downsampledWidth = DOWNSAMPLED_WIDTH;
    if( argc > 1 ) {
        originalWidth = atoi( argv[1] );
        if( argc > 2 )
            downsampledWidth = atoi( argv[2] );
    }

    downsampledData = (double *) malloc( downsampledWidth * sizeof(double) );

    /*
     In addition to space to hold the original data to be downsampled, we also allocate
     some buffer space "left" and "right" of the original data and fill it with dummy data.
     This buffer space / dummy data on the side makes the convolution algorithm cleaner as
     we don't have to worry about the "edges" of the convolution window.
    */
    originalData = (double *) malloc( (CONVOLUTION_SPAN * 2 + originalWidth) * sizeof(double) );

    // Initialize the original data
    ptr1 = originalData;
    switch ( initializationType ) {
        case step:
            // Fill everything, including the left and right buffer zones, with initial data
            for ( i = 0; i < (CONVOLUTION_SPAN * 2 + originalWidth) / 2; i++ )
                *ptr1++ = 1;
            for (  ; i < (CONVOLUTION_SPAN * 2 + originalWidth); i++ )
                *ptr1++ = 2;
            break;

        case pulse:
            for ( i = 0; i < (CONVOLUTION_SPAN * 2 + originalWidth); i++ )
                *ptr1++ = 1;
            originalData[ (CONVOLUTION_SPAN * 2 + originalWidth) / 2 ] = 2;
            break;

        case square:
            s = originalWidth / 6;
            for ( i = 0; i < (CONVOLUTION_SPAN * 2 + originalWidth); i++ )
                *ptr1++ = 1;
            ptr1 = originalData + CONVOLUTION_SPAN + s;
            for ( i = 0; i < 3; i++ ) {
                for ( j = 0; j < s; j++ )
                    *ptr1++ = 2;
                ptr1 += s;
            }
            break;

        case sine:
            scaleFactor = 6 * M_PI / originalWidth;
            for ( i = -CONVOLUTION_SPAN; i < (CONVOLUTION_SPAN + originalWidth); i++ )
                *ptr1++ = sin( i * scaleFactor );
            break;

        default:
            fprintf(stderr, "Unsupported initialization function type" );
            exit( -1 );
    }

    // Convolve and downsample the original data into a new, smaller vector
    scaleFactor = (M_PI * downsampledWidth ) / originalWidth;
    ptr2 = downsampledData;
    for ( d = 0; d < downsampledWidth; d++) {
        accumulator = 0;
        i = (.5 + originalWidth * d / downsampledWidth );
        ptr1 = originalData + i;
        for ( i = -CONVOLUTION_SPAN; i <= CONVOLUTION_SPAN; i++) {
            x = i * scaleFactor;
            if ( fabs(x) < 1e-6 )
                accumulator += (*ptr1++ * downsampledWidth) / originalWidth;
            else
                accumulator += *ptr1++ * sin( x ) / (i * M_PI);
        }
        *ptr2++ = accumulator;
    }

    // Output the results
    ptr2 = downsampledData;
    for ( i = 0; i < downsampledWidth; i++)
        printf( "%.9lf\n", *ptr2++ );

    return 0;
}