Could anyone help on how to call this function which will compute a moving median from an input array and return an output array with replaced median values
Code:
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;
}
into this simple program which uses sndfile.h
Code:
/*
#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 */
All i'm trying to do is replace the process _data function call with the WaveBuf::BMedian35 function, it's just im not sure what parameters to use, i realise some code will have to be changed in the function such as the m_waveFormat.scan_count, but its just the correct use of the function that im worried about,
Anyones help is greatly appreciated