Mutante
Mutante

Reputation: 278

Moving Average in C Language

I am trying to do a moving average filter in C language, I've adapted a matlab program that works correctly, the input of my filter is a .pcm archive (a sweep audio signal), the problem to me is the output archive of the moving average in C language, the output goes wrong, the signal only decreases along the time (don't filtter).

Below My C code:

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



#define m 16  // MEDIA LENGTH

int main()
{
    short *x;                              // X VECTOR = 0
    double coef[m];                           // COEF VECTOR = 0
    int i;                                  // FOR COUNT VARIABLE
    int n;                                  // FOR COUNT VARIABLE
    int itera;                              // LENGTH OF THE INPUT FILE
    double aux;                                // AUX VARIABLE TO MAKE THE SUM
    double *y;                              // AUX VECTOR TO RECEIVE X VECTOR
    double *saida;                          // OUTPUT VECTOR
    FILE *arq;                                 // POINTER TO THE FILE

    for (i = 0; i < m; i++) {
        coef[i] = 1.0/m;
    }
    arq = fopen("sweep_100_3400.pcm", "rb");
    if (arq == NULL)
    {
        printf("ERROR WHILE OPENING THE FILE\n");
        return;
    }
    // compute size of file
    fseek(arq,0,SEEK_END);
    itera = ftell(arq)/sizeof(short);
    rewind(arq);

    // alloc mem for x, read the vector from input file
    x = malloc(itera*sizeof(short));
    fread(x,sizeof(short),itera,arq);
    fclose(arq);
    // alloc mem for y

    y = malloc(itera*sizeof(double));
    saida = malloc(itera*sizeof(double));

    for (i=0; i < itera; i++) {
        y[0] = x[i];
        aux=0;
        for (n=0; n<m; n++){
            aux += coef[n] * y[n];
        }
        saida[i]=aux;
        for (n=m; n <2; n--){
            x[n] = x[n-1];
        }
    }
    arq=fopen("saida_medial_movel_c.pcm","wb");
    fwrite(saida, sizeof(double),itera,arq);
    fclose(arq);
    free(saida);
    free(y);
    free(x);
}

The image below is the output of the matlab program to moving average with length 16:

enter image description here

This image is the output in C language with moving average with length 16:

enter image description here

Someone knows what may be?

Below the code in matlab, that I've adapted:

%MOVING AVERAGE EXAMPLE SCRIPT
clear all;
close all;
clc;
% DEFINES MEDIA LENGTH
m = 16;
%VECTOS EQUAL ZERO
x = zeros (m,1); 
coef = zeros (m,1);
%INITIALIZE VECTOR 
for j = 1 : m,
        coef (j,1) = 1/m;
end
%READ INPUT FILE
fid = fopen ('sweep_100_3400.pcm','rb');
s = fread (fid, 'int16');
fclose(fid);
subplot(2,1,1);
plot(s);
grid on;
title('ENTRADA DO FILTRO');
%PROCESS
itera = length(s);
sav_y = zeros (itera, 1);
%EXECUTE PROCESS
for j = 1 : itera,
    x(1,1) = s (j,1);
    %PRODUCTS SUM
    y=0;
    for n = 1 : m,
        y = y + coef(n,1) * x(n,1);
    end
    sav_y(j,1) = y;
    %SHIFT THE VECTOR
    for n = m: -1 : 2,
        x (n,1) = x(n-1,1);
    end 
end 
%PLOT OUTPUT
subplot (2,1,2);
plot (sav_y);
grid on;
title('SAÍDA DO FLITRO');
%SAVE THE OUTPUT IN ANOTHER FILE
fid = fopen('saida_mm_manual.pcm','wb');
fwrite(fid,sav_y,'int16');
fclose(fid);

Update 1 (Using the answer above):

enter image description here

The begining of the signal still with interference, but from the middle to the end of the signal is correctly.

Upvotes: 2

Views: 7513

Answers (1)

Sultan
Sultan

Reputation: 138

There are a few issues with your main signal processing loop.

for (i=0; i < itera; i++) {
    y[0] = x[i];
    aux=0;
    for (n=0; n<m; n++){
        aux += coef[n] * y[n];
    }
    saida[i]=aux;
    for (n=m; n <2; n--){
        x[n] = x[n-1];
    }
}

Every element of your coef array is the same, so you may as well make it a single constant value of 1.0 / m. For the y array, you are not setting any of its elements apart from the first. Thus your loop adds a constant multiplied by an uninitialized value 16 times over to aux. This is what produces the garbage output. I'm not sure what your (n=m; n <2; n--) loop is supposed to do, but it will never even run since m > 2.

A (inefficient) simple moving average would look something more like this:

for (i = 0; i < itera - m; i++) {
    aux = 0;
    for (n = 0; n < m; n++){
        aux += x[i+n];
    }
    saida[i] = aux * (1.0 / m);
}

A more efficient version could avoid repeatedly processing intermediate elements, and instead just add the new number entering the sliding window and subtract the old number leaving the window each iteration. If dealing with floating point numbers, care would have to be taken with numerical stability when handling ill conditioned numbers etc, but that's a whole different problem that you don't have to concern yourself about for now.

Upvotes: 2

Related Questions