Double rounding when writing to file leads to accumulating errors

[zsoltlazar]

After a read-write-read cycle of an EDF file the sampel values get changed by a significant margin.
Digging up the responsible lines in the code (4604 and 4606 in _extensions/c/edflib.c -> edfwrite_physical_samples()) suggest that the two subsequent roundings might be the cause.
As a proof I included a simple code that compares the present implementation with the suggested (correct) one.

`#include <stdio.h>

int f1(double x, double p, double P, int d, int D) {
/*
converting double to int
implementation according to edflib.c
*/
int value;
double bitvalue, offset;
bitvalue = (P - p)/(D - d);
offset = P/bitvalue - D;
value = x/bitvalue; // rounding 1
value -= offset; // rounding 2
return value;
}

int f2(double x, double p, double P, int d, int D) {
/*
converting double to int
corrected implementation
*/
int value;
double bitvalue, offset;
bitvalue = (P - p)/(D - d);
offset = P/bitvalue - D;
value = x/bitvalue - offset; // rounding
return value;
}

double inv(int v, double p, double P, int d, int D) {
/*
converting int to double
implementation according to edflib.c
/
double bufvalue;
double bitvalue, offset;
bitvalue = (P - p)/(D - d);
offset = P/bitvalue - D;
bufvalue = bitvalue(offset + v);
return bufvalue;
}

int main() {
int D = 2047; // dig_max
int d = -2048; // dim_min
double P = 249; // phys_max
double p = -251; // phys_min
double x1, x2; // signal values
int v1, v2; // stored values

    printf("Show that the two double->int conversions are not equivalent\n"$
    for(x1 = -1; x1 < 1; x1 += 0.13) {
            printf("x: %5f\t", x1);
            printf("f1: %d ", f1(x1, p, P, d, D));
            printf("f2: %d\n", f2(x1, p, P, d, D));
    }

    printf("Show that the edflib.c implementation (with two roundings)\n");
    printf("can lead to the accumulation of errors upon subsequent\n");
    printf("back and forth conversions\n");
    x1 = x2 = 0.3;
    for (int i = 0; i < 10; i++) {
            v1 = f1(x1, p, P, d, D); //  double->int the wrong way
            x1 = inv(v1, p, P, d, D);//  int->double
            v2 = f2(x2, p, P, d, D); //  double->int the right way
            x2 = inv(v2, p, P, d, D);//  int->double
            printf("x1: %4f\tv1: %d\t x2: %4f\tv2: %d\n", x1, v1, x2, v2); $
    }
    return 0;

}

`

[holgern]

Thank you for reporting this
I added a unit test to comfirm that it is working now