C 编程：移动平均滤波器

Question

我试图在 C 中实现一个等效的 Matlab 函数 smooth(y, span)。 该函数的Matlab代码为：

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];

这是我的 C 代码：

#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;}

我获得了与 Matlab 代码相同的值，直到代码中止的信号的某个索引。有谁知道为什么以及如何解决它？

提前致谢！

Answer 1

如果您需要在cumulative_sum 中索引cum_vector[size]，那么您需要确保cum_vector 具有size+1 元素。在您的代码中它只有size。因此，您的书写越界，这可能会导致稍后崩溃。

在moving_average 中，您还可以访问cum_x_tmp[size]。

按如下方式分配您的数组：

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

在 Valgrind 或类似的内存检查器下运行你的程序会指出这个问题，以及另一个答案中指出的内存泄漏。

Answer 2

你有内存泄漏：

在每次调用moving_average 时，都会泄漏为cum_x_tmp 分配的内存：

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

当你不再需要它时，你需要释放它。 此外，您应该对从malloc 返回的值添加一个空检查，以检测分配失败的时间。否则，您将在释放空指针时调用未定义行为。

Answer 3

清理了原始代码以实际编译和正常工作。又名 Matlab 的 smooth() 函数。如果您发现有一种方法可以使此过程更快和/或更安全，请回帖。

#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;
}