Moteurs audio, 32 bits / 64 bits, distinguer le vrai du faux .

#! /usr/bin/env python                                                                   

# Last Change: Mon Oct 30 01:00 PM 2006 J                                                



import numpy as N

    

n   = int(1e5)

dbd = -80

    

# A is double, between -1 and 1 (0dB)                                                    

A   = 2*N.random.rand((n))-1                                                             

# B is double, between -1e-4 and 1e4

B   = (2*N.random.rand((n))-1) * 10 ** (dbd/20)                                          

    

S   = A + B

        

powerA  = N.sqrt(N.sum(A**2)) * 1./n                                                     

powerB  = N.sqrt(N.sum(B**2)) * 1./n                                                     

print "Real db Diff (power) is %f" % (20*N.log10(powerA/powerB))                         

    

A32 = N.float32(A)

B32 = N.float32(B)                                                                       

    

S32     = A32 + B32

    

# numpy broadcast to 64 bits automatically, so the difference is done

# in 64 bits

diff    = S32 - S



powerdiff   = N.sqrt(N.sum(diff ** 2)) * 1./n

maxdiff     = N.max(N.abs(diff))

                                                                                         

print "average diff of dynamic compared to bigger signal (0 dB) is %f" % (20*N.log10(powerdiff))

print "worse dynamic compared to bigger signal (0 dB) is %f" % (20*N.log10(maxdiff))

Citation :

/*

 * Last Change: Mon Oct 30 01:00 PM 2006 J

 *

 * Small test to see pratical difference beween summing in 32 and summing in 64 bits

 * on the platform where this program is run

 */



/*

 * For ISO C, all floating functions are for double by default.

 * If they have the suffix f, they are float

 */



/* For FLT_ constants */

#include <float.h>



/* for drand48 (SVID standart) and exit */

#include <stdlib.h> 

/* for log10 */

#include <math.h>



#include <stdio.h>



int do_compare(size_t nelements, double dbd);

int generate_signals(double *in, double *in2, size_t nel, double dbdiff);



int lev2dB(double ref, double val, double* ret);

double db2coeff(double dblev);

double avdbdiff(const double* ref, const double* val, size_t nel);

double avdbref(const double* in, size_t nel, double ref);

double mindbref(const double* in, size_t nel, double ref);



long double power(const double*in, size_t nel);



int sum32(const float* in1, const float* in2, float* out, const size_t n);

int sum64(const double* in1, const double* in2, double* out, const size_t n);

int diff64(double* inout, const double *in2, size_t nel);



int double2float(const double* in, float *out, size_t nel);

int float2double(const float* in, double* out, size_t nel);



/* alloc or die */

void*    xmalloc(size_t n);



/* Not sure about precision, but should not matter, since 

 * it is used only for db computation*/

double exp10(double in)

{

    return exp(in * log(10));

}



/* Entry point */

int main(void)

{

    int     status;

    size_t  nel = 10000000;

    double  dbdiff  = -80;



    status  = do_compare(nel, dbdiff);

    if (status < -1) {

        fprintf(stderr, "Error while comparing...n");

        return -1;

    }

    return 0;

}



int do_compare(size_t nel, double dbd)

{

    double  *a, *b, *s, *s64;

    float   *a32, *b32, *s32;

    int     status;

    double  sumdbdiff;



    /*

     * Allocations

     */

    a   = xmalloc(sizeof(*a) * nel);

    b   = xmalloc(sizeof(*b) * nel);

    s   = xmalloc(sizeof(*s) * nel);

    a32 = xmalloc(sizeof(*a32) * nel);

    b32 = xmalloc(sizeof(*b32) * nel);

    s32 = xmalloc(sizeof(*s32) * nel);

    s64 = xmalloc(sizeof(*s64) * nel);



    srand48(time(0));

    /*

     * Generate test signals

     */

    status  = generate_signals(a, b, nel, dbd);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }



    /*

     * Make a copy of test signals into 32 bits buffers

     */

    status  = double2float(a, a32, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }

    status  = double2float(b, b32, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }



    /*

     * Do the summation in 32 bits

     */

    status  = sum32(a32, b32, s32, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }

    status  = float2double(s32, s64, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }



    /*

     * Do the summation in 64 bits

     */

    status  = sum64(a, b, s, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }



    /*

     * Compute the difference (in 64 bits) (in place into s)

     */

    status  = diff64(s, s64, nel);

    if (status) {

        fprintf(stderr, "%s:%s, %d errorn", __FILE__, __func__, __LINE__);

        exit(EXIT_FAILURE);

    }



    /*

     * Compare the results

     */

    sumdbdiff = avdbref(s, nel, 1.0);

    printf("average db diff between 32 bits summing and 64 summing is %fn", sumdbdiff);

    sumdbdiff = mindbref(s, nel, 1.0);

    printf("worse db diff between 32 bits summing and 64 summing is %fn", sumdbdiff);

#if 0

#endif



    /* Those free are just for checking under valgrind */

    free(a);

    free(b);

    free(s);

    free(a32);

    free(b32);

    free(s32);

    free(s64);



    return 0;

}



/* compute the power of in: 1/nel * sqrt(sum(in ** 2)) */

long double power(const double*in, size_t nel)

{

    long double acc = 0;

    size_t      i;



    for(i = 0; i < nel; ++i) {

        acc += in[i] * in[i];

    }



    return sqrtl(acc) / nel;

}



int diff64(double* inout, const double *in2, size_t nel)

{

    size_t      i;



    for(i = 0; i < nel; ++i) {

        inout[i] -= in2[i];

    }



    return 0;

}



/* Generate (uniform) random signals in and in2. in2 has dbdiff dB difference on average 

 * compare to in1, and in1 is normalized in the range [-1, 1] */

int generate_signals(double *in1, double *in2, size_t nel, double dbdiff)

{

    size_t  i;

    double  amp = db2coeff(dbdiff);



    for(i = 0; i < nel; ++i) {

        in1[i]  = 2 * drand48() - 1;

        in2[i]  = amp * (2 * drand48() - 1);

    }



    printf("Asked level is %f, real level on average is %fn", dbdiff, avdbdiff(in1, in2, nel));



    return 0;

}



/* Convert a and b content to 32 bits and put the results into a32 and b 32.

 * I am not sure about the correct way to do the convertion... */

int double2float(const double* in, float *out, size_t nel)

{

    size_t  i;

    for(i = 0; i < nel; ++i) {

        out[i]  = (float)in[i];

    }



    return 0;

}



int float2double(const float* in, double* out, size_t nel)

{

    size_t  i;



    for(i = 0; i < nel; ++i) {

        out[i]  = (double)in[i];

    }



    return 0;

}



/* give the *average* db level compared to a reference 0 dB */

double avdbref(const double* in, size_t nel, double ref)

{

    long double pow;

    double      dbl;

    

    pow = power(in, nel);



    lev2dB(ref, (double)pow, &dbl);



    printf("ref is %10e, pow is %10e, level is %fn", 

            ref, (double)pow, dbl);



    return dbl;

}



/* give the min db level compared to a reference 0 dB */

double mindbref(const double* in, size_t nel, double ref)

{

    double  maxval;

    double  tmp;

    double  dbl;

    size_t  i;



    maxval  = 0;

    for(i = 0; i < nel; ++i) {

        tmp = fabs(in[i]);

        if ( tmp > maxval) {

            maxval  = tmp;

        }

    }

    

    lev2dB(ref, maxval, &dbl);



    printf("ref is %10e, max is %10e, level is %fn", 

            ref, maxval, dbl);



    return dbl;

}



/* give the *average* db difference between two signals */

double avdbdiff(const double* ref, const double* val, size_t nel)

{

    long double pow1;

    long double pow2;

    double      dbl;

    

    pow1    = power(ref, nel);

    pow2    = power(val, nel);



    lev2dB((double)pow1, (double)pow2, &dbl);



    printf("pow1 is %10e, pow2 is %10e, dbdiff is %fn", 

            (double)pow1, (double)pow2, dbl);



    return (dbl);

}



/* put the db value with reference level ref into ret

 * return 0 on success, -1 otherwise */

int lev2dB(double ref, double val, double* ret)

{

    if (fabs(val) < DBL_EPSILON || val < DBL_MIN) {

        fprintf(stderr, "%s:%s, line %d, val (%f) is below epsilon or negativen", 

                __FILE__, __func__, __LINE__, val);

        return -1;

    }

    *ret    = 20*log10(ref / val);



    return 0;

}



/* given a dB level difference, returns the amplitude coeff */

double db2coeff(double dblev)

{

    return exp10(dblev / 20.);

}



/* sum in 32 bits */

int sum32(const float* in1, const float* in2, float* out, const size_t n)

{

    size_t  i;

    for(i = 0; i < n; ++i) {

        out[i]  = in1[i] + in2[i]; 

    }



    return 0;

}



/* sum in 64 bits */

int sum64(const double* in1, const double* in2, double* out, const size_t n)

{

    size_t  i;

    for(i = 0; i < n; ++i) {

        out[i]  = in1[i] + in2[i]; 

    }



    return 0;

}



void* xmalloc(size_t n)

{

    void*   tmp;

    tmp     = malloc(n);

    if (tmp == NULL) {

        fprintf(stderr, "Error while allocatingn");

        exit(EXIT_FAILURE);

    }

    return tmp;

}