void CWaveBuf::BMedian35(short* lpSource, short* lpDest, short cx)
{	
	short offsetnextpoint = m_waveFormat.scan_count;// number of chans within source buffer
	short nbspan = TransformBufferSpan[13];	// nb of points within moving windows /2
	// usually we use = 30; ie 61 points total. Our spikes are defined over 60 points.
	short offsetspan = offsetnextpoint * nbspan;	// offset between center of window & end
	short NN = nbspan*2 +1;			// nb of points within window
	short *parray = new short[NN];		// temporary array 
	
	// init parray = consecutive points (image of the data points)	
	int i=0;					// index variable
	short* lp = lpSource;			// pointer to origin of source buffer
	lp -= offsetspan;				// first point of the moving window
	for (i = 0; i < NN; lp += offsetnextpoint, i++)	// load parray
		*(parray+i) = *lp;			// with buffer values

	// sort parray into ascending order using heapsort algorithm
	// 	// "l"= index which will be decremented from initial value down to 0 during
	// 'hiring' (heap creation) phase. Once it reaches 0, the index "ir" will be 
	// decremented from its initial value down to 0 during the 'retirement-and-
	// promotion' (heap selection) phase.
	int l = nbspan + 1;				// temp index
	int ir = NN - 1;				// temp index
	int j;					// temp index
	short val;					// temp storage	

	for (;;)					// pseudo-loop over parray
	{
		// -------------------------
		if (l > 0)				// still in hiring phase
		{
			l--;
			val = *(parray + l);
		}
		else				// in retirement-and-promotion phase
		{
			val = *(parray+ir);		// clear a space at the end of the array
			*(parray+ir) = *(parray);	// retire the top of the heap into it
			ir--;			// decrease the size of the corporation
			if (ir == 0)		// done with the last promotion
			{
				*(parray) = val;	// the least competent worker of all
				break;		// exit the sorting algorithm
			}
		}
		// -------------------------
		i = l + 1;				// wether we are in the hiring or promotion, we
		j = i + i;				// here set up to sift down element raa to its 
						// proper level
		while (j-1 <= ir)
		{			
			if (j-1 < ir)				
				if (*(parray+j-1) < *(parray+j))
					j++;	// compare to the better underlining
			if (val < *(parray+j-1))	// demote val
			{
				*(parray+i-1) = *(parray+j-1);
				i = j;
				j = j + j;
			}			// this is val's level. Set j to terminate the
			else			// sift-down
				j = ir + 2;
		}
		*(parray+i-1) = val;			// put val into its slot
	}

	// end of initial sort

	short newvalue;				// temp variable used in the loop
	short oldval;				// temp variable used in the loop
	short offset2span = 2*offsetspan+offsetnextpoint;// offset first to last pt/window
	lp= lpSource + offsetspan;


	for (cx; cx>0; cx--, lp += offsetnextpoint, lpDest++)
	{
		newvalue = *lp;			// new value to insert into array
		oldval = *(lp-offset2span);		// value to discard from array

		// locate position of old value to discard
		// use bisection - cf Numerical Recipes pp 90
		// on exit, j= index of oldval

		// 
		int jhigh = NN-1;			// upper index
		int jlow = 0;			// mid point index
		while (jlow <= jhigh)
		{
			j = (jlow+jhigh)/2;			
			if (oldval > *(parray+j))
				jlow = j +1;
			else if (oldval < *(parray+j))
				jhigh = j-1;
			else
				break;
		}

		// insert new value in the correct position

		// case 1: search (and replace) towards higher values
		if (newvalue > *(parray + j))
		{	
			for (j; newvalue > *(parray + j); j++)
			{								
				if (j == NN)
					break;
				*(parray+ j) = *(parray + j + 1);
			} 
			*(parray + j-1) = newvalue;
		}

		// case 2: search (and replace) towards lower values
		else if (newvalue < *(parray + j))
		{
			for (j; newvalue < *(parray + j); j--)
			{
				if (j == 0)
				{
					if (newvalue < *parray)
						j--;
					break;
				}
				*(parray + j)  = *(parray + j -1);
			}
			*(parray + j+1) = newvalue;
		}

		// case 3: already found!
		else
			*(parray + j) = newvalue;


		// save median value in the output array
		// m_waveFormat.binzero : binary value of the zero volt. For ex, for 12 bits binary encoded data, zero is
		// encoded as 2048.
		*lpDest = *(lp-offsetspan) - *(parray+nbspan) + (short) m_waveFormat.binzero;
	}

	delete parray;
}

/*


#include	<stdio.h>

/* Include this header file to use functions from libsndfile. */
#include	<sndfile.h>

/*    This will be the length of the buffer used to hold.frames while
**    we process them.
*/
#define		BUFFER_LEN	1024

/* libsndfile can handle more than 6 channels but we'll restrict it to 6. */
#define		MAX_CHANNELS	6

/* Function prototype. */
static void process_data (double *data, int count, int channels) ;


int
main (void)
{   /* This is a buffer of double precision floating point values
    ** which will hold our data while we process it.
    */
    static double data [BUFFER_LEN] ;

    /* A SNDFILE is very much like a FILE in the Standard C library. The
    ** sf_open function return an SNDFILE* pointer when they sucessfully
	** open the specified file.
    */
    SNDFILE      *infile, *outfile ;

    /* A pointer to an SF_INFO stutct is passed to sf_open.
    ** On read, the library fills this struct with information about the file.
    ** On write, the struct must be filled in before calling sf_open.
    */
    SF_INFO		sfinfo ;
    int			readcount ;
    const char	*infilename = "input.wav" ;
    const char	*outfilename = "output.wav" ;

    /* Here's where we open the input file. We pass sf_open the file name and
    ** a pointer to an SF_INFO struct.
    ** On successful open, sf_open returns a SNDFILE* pointer which is used
    ** for all subsequent operations on that file.
    ** If an error occurs during sf_open, the function returns a NULL pointer.
	**
	** If you are trying to open a raw headerless file you will need to set the
	** format and channels fields of sfinfo before calling sf_open(). For
	** instance to open a raw 16 bit stereo PCM file you would need the following
	** two lines:
	**
	**		sfinfo.format   = SF_FORMAT_RAW | SF_FORMAT_PCM_16 ;
	**		sfinfo.channels = 2 ;
    */
    if (! (infile = sf_open (infilename, SFM_READ, &sfinfo)))
    {   /* Open failed so print an error message. */
        printf ("Not able to open input file %s.\n", infilename) ;
        /* Print the error message from libsndfile. */
        puts (sf_strerror (NULL)) ;
        return  1 ;
        } ;

    if (sfinfo.channels > MAX_CHANNELS)
    {   printf ("Not able to process more than %d channels\n", MAX_CHANNELS) ;
        return  1 ;
        } ;
    /* Open the output file. */
    if (! (outfile = sf_open (outfilename, SFM_WRITE, &sfinfo)))
    {   printf ("Not able to open output file %s.\n", outfilename) ;
        puts (sf_strerror (NULL)) ;
        return  1 ;
        } ;

    /* While there are.frames in the input file, read them, process
    ** them and write them to the output file.
    */
    while ((readcount = sf_read_double (infile, data, BUFFER_LEN)))
    {   process_data (data, readcount, sfinfo.channels) ;
        sf_write_double (outfile, data, readcount) ;
        } ;

    /* Close input and output files. */
    sf_close (infile) ;
    sf_close (outfile) ;

    return 0 ;
} /* main */

static void
process_data (double *data, int count, int channels)
{	double channel_gain [MAX_CHANNELS] = { 0.5, 0.8, 0.1, 0.4, 0.4, 0.9 } ;
    int k, chan ;

    /* Process the data here.
    ** If the soundfile contains more then 1 channel you need to take care of
    ** the data interleaving youself.
    ** Current we just apply a channel dependant gain.
    */

    for (chan = 0 ; chan < channels ; chan ++)
        for (k = chan ; k < count ; k+= channels)
            data [k] *= channel_gain [chan] ;

    return ;
} /* process_data */