C programming : Moving average filter

Question

I have tried to implement an equivalent in C of the Matlab function smooth(y, span). The Matlab code of the function is :

n = length(y);
span = min(span,n);
width = span-1+mod(span,2); % force it to be odd
c = filter(ones(width,1)/width,1,y);
cbegin = cumsum(y(1:width-2));
cbegin = cbegin(1:2:end)./(1:2:(width-2))';
cend = cumsum(y(n:-1:n-width+3));
cend = cend(end:-2:1)./(width-2:-2:1)';
c = [cbegin;c(width:end);cend];

And here is my C code :

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

// Compute cumulative sum of a signal of a given size.
void cumulative_sum(float *cum_vector, float *signal, int size)
{
    int i;
    cum_vector[0]=0;
    cum_vector[1]=signal[0];
    for(i=2;i<=size;i++)
    {
        cum_vector[i] = cum_vector[i-1] + signal[i-1];
    }
 }

 // Moving average filter of a signal of a given size with a pre-defined span.
 void moving_average(float *vector, int size, int span)
{
if (span > size)
{
    span = size; //span = min(span,size);
}
int i;
int j=0;
int width = span - 1 + (span % 2); //force it to be odd
float *cum_x_tmp;
cum_x_tmp = (float *)malloc(size*sizeof(float));
cumulative_sum(cum_x_tmp,vector,size);

for (i=0;i<size;i++)
{
    if (i< (width - 2)/2 +((width-2)%2 !=0))
    {
        vector[i] = cum_x_tmp[2*i+1]/(2*i+1);
    }
    else if (i>=(width - 2)/2 +((width-2)%2 !=0) && i<= (width - 2)/2 +((width-2)%2 !=0) + size - width)
    {
        int ind1 =i + floor(width/2)+1;
        int ind2 =i - floor(width/2)-1+1;
        vector[i] = (cum_x_tmp[ind1] - cum_x_tmp[ind2])/width;
    }
    else if (i>(width - 2)/2 +((width-2)%2 !=0) + size - width)
    {
        int ind4 = size-1-width+3-1+2*j+1;
        vector[i] = (cum_x_tmp[size] - cum_x_tmp[ind4])/(width -2*(j+1));
        j++;
    }
}
}int main(int argc, const char * argv[]) {

int N=900;
float data;
FILE *fp;
fp = fopen("file.txt","r");

float *signal;
signal = (float*)malloc(N*sizeof(float));

int i;
for (i=0;i<N;i++)
{
    fscanf(fp, "%f", &data);
    signal[i]=data;
}
fclose(fp);

moving_average(signal,N,50);

free(signal);

return 0;}

I obtained the same values as the Matlab code until a certain index of the signal at which the code abords. Does anyone know why and how to fix it ?

Thanks in advance !

Answer 1

The original code cleaned up to actually compile and work correctly. aka Matlab's smooth() function. Post back if you see a way to make this faster and/or safer.

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

void get_cum_sum(double *cum_vector, double *signal, int size)
{
   int i;
   cum_vector[0] = signal[0];
   for(i = 1; i < size; i++)
   {
      cum_vector[i] = signal[i] + cum_vector[i-1];
   }
}

void smooth(double *vector, int size, int span)
{
   int i, window, half_window, tmp_index;
   double *cum_x_tmp;

   if (span > size)
   {
      span = size;
   }
   window = span - 1 + (span % 2);
   half_window = window / 2;

   cum_x_tmp = (double *)malloc((size) * sizeof(double));

   get_cum_sum(cum_x_tmp, vector, size);

   for (i = 0; i < size; i++)
   {
      if (i <= half_window)
      {
         vector[i] = (cum_x_tmp[2 * i]) / (2 * i + 1);
      }
      else if (i >= (size - half_window))
      {
         tmp_index = (size - (i + 1)) * 2 + 1;
         vector[i] = (cum_x_tmp[size - 1] - cum_x_tmp[i - (tmp_index / 2) - 1]) / tmp_index;
      }
      else
      {
         vector[i] = (cum_x_tmp[i + half_window] - cum_x_tmp[i - half_window - 1]) / window;
      }
   }
   free(cum_x_tmp);
}

int main(int argc, const char * argv[])
{
   int i;
   int N = 7;
   double data;
   FILE *fp;
   double *signal;

   signal = (double*)malloc(N*sizeof(double));
   fp = fopen("data.txt","r");

   for (i = 0; i < N; i++)
   {
      fscanf(fp, "%lf", &data);
      signal[i] = data;
   }

   fclose(fp);

   smooth(signal, N, 7);

   for (i = 0; i < N; i++)
   {
      printf("signal[%d] = %f\n", i, signal[i]);
   }

   free(signal);

   return 0;
}

Answer 2

If you need to index cum_vector[size] in cumulative_sum, then you need to make sure that cum_vector has size+1 elements. In your code it has only size. Thus, you are writing out of bounds, which is likely to cause a crash later on.

In moving_average you also access cum_x_tmp[size].

Allocate your array as follows:

cum_x_tmp = (float *)malloc((size+1)*sizeof(float));

Running your program under Valgrind or a similar memory checker will point out this issue, as well as the memory leak pointed out in the other answer.

Answer 3

You've got a memory leak:

On each invocation of moving_average, you are leaking the memory allocated for cum_x_tmp:

    float *cum_x_tmp;    
    cum_x_tmp = (float *)malloc(size*sizeof(float));

You need to free that memory when you no longer need it. Also, you should add a null-check to the value returned from malloc, to detect when allocation fails. Otherwise, you'll invoke Undefined Behavior when deallocating the null pointer.