# moving median function

• 05-04-2006
supermeew
moving median function
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