abx.c

音频编码

/*
 *  Usage: abx original_file test_file
 *
 *  Ask you as long as the probability is below the given percentage that
 *  you recognize differences
 *
 *  Example: abx music.wav music.mp3
 *           abx music.wav music.mp3 --help
 *
 *  Note: several 'decoding' utilites must be on the 'right' place
 *
 *  Bugs: 
 *      fix path of decoding utilities
 *      only 16 bit support
 *      only support of the same sample frequency
 *      no exact WAV file header analysis
 *      no mouse or joystick support
 *      don't uses functionality of ath.c
 *      only 2 files are comparable
 *      worse user interface
 *      quick & dirty hack
 *      wastes memory
 *      compile time warnings
 *      buffer overruns possible
 *      no dithering if recalcs are necessary
 *      correlation only done with one channel (2 channels, sum, what is better?)
 *      lowpass+highpass filtering (300 Hz+2*5 kHz) before delay+amplitude corr
 *      cross fade at start/stop
 *      non portable keyboard
 *      fade out on quit, fade in on start
 *      level/delay ajustment should be switchable
 *      pause key missing
 *      problems with digital silence files (division by 0)
 *      Gr鲞e cross corr fenster 2^16...18
 *      Stellensuche, ab 0*len oder 0.1*len oder 0.25*len, nach Effektiv oder Spitzenwert
 *      Absturz bei LPAC feeding, warum?
 *      Als 'B' beim Ratespiel sollte auch '0'...'9' verwendbar sein 
 *      Oder mit einem Filter 300 Hz...3 kHz vorher filtern?
 *      Multiple encoded differenziertes Signal
 *      Amplitudenanpassung schaltbar machen?
 *      Direkt auf der Kommandozeile kodieren:
 *      abx "test.wav" "!lame -b128 test.wav -"
 */

// If the program should increase it priority while playing define USE_NICE. 
// Program must be installed SUID root. Decompressing phase is using NORMAL priority
#define USE_NICE

// Not only increase priority but change to relatime scheduling. Program must be installed SUID root
#define USE_REALTIME

// Path of the programs: mpg123, mppdec, faad, ac3dec, ogg123, lpac, shorten, MAC, flac
//#define PATH_OF_EXTERNAL_TOOLS_FOR_UNCOMPRESSING   "/usr/local/bin/"
#define PATH_OF_EXTERNAL_TOOLS_FOR_UNCOMPRESSING   ""


#if defined HAVE_CONFIG_H
# include <config.h>
#endif

#include <assert.h>
#include <ctype.h>
#include <fcntl.h>
#include <limits.h>
#include <math.h>
#include <memory.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <termios.h>
#include <time.h>
#include <unistd.h>
#include <sys/ioctl.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <sys/types.h>

#define  MAX  (1<<17)

#if   defined HAVE_SYS_SOUNDCARD_H
# include <sys/soundcard.h>
#elif defined HAVE_LINUX_SOUNDCARD_H
# include <linux/soundcard.h>
#else
# include <linux/soundcard.h>         /* stand alone compilable for my tests */
#endif

#if defined USE_NICE
# include <sys/resource.h>
#endif
#if defined USE_REALTIME
# include <sched.h>
#endif

#define  BF           ((freq)/25)
#define  MAX_LEN      (210 * 44100)
#define  DMA_SAMPLES  512	    	/* My Linux driver uses a DMA buffer of 65536*16 bit, which is 32768 samples in 16 bit stereo mode */

void  Set_Realtime ( void )
{
#if defined USE_REALTIME
    struct sched_param  sp;
    int                 ret;

    memset      ( &sp, 0, sizeof(sp) );
    seteuid     ( 0 );
    sp.sched_priority = sched_get_priority_min ( SCHED_FIFO );
    ret               = sched_setscheduler ( 0, SCHED_RR, &sp );
    seteuid     ( getuid() );
#endif

#if defined USE_NICE
    seteuid     ( 0 );
    setpriority ( PRIO_PROCESS, getpid(), -20 );
    seteuid     ( getuid() );
#endif
}

int verbose = 0;

static struct termios stored_settings;


void reset ( void )
{
    tcsetattr ( 0, TCSANOW, &stored_settings );
}


void set ( void )
{
    struct termios new_settings;
    
    tcgetattr ( 0, &stored_settings );
    new_settings = stored_settings;
    
    new_settings.c_lflag    &= ~ECHO;
    /* Disable canonical mode, and set buffer size to 1 byte */
    new_settings.c_lflag    &= ~ICANON;
    new_settings.c_cc[VTIME] = 0;
    new_settings.c_cc[VMIN]  = 1;
    
    tcsetattr(0,TCSANOW,&new_settings);
    return;
}


int sel ( void )
{
    struct timeval  t;
    fd_set          fd [1];
    int             ret;
    unsigned char   c;

    FD_SET (0, fd);    
    t.tv_sec  = 0;
    t.tv_usec = 0;

    ret = select ( 1, fd, NULL, NULL, &t );
    
    switch ( ret ) {
    case  0: 
        return -1;
    case  1: 
        ret = read (0, &c, 1);
        return ret == 1  ?  c  :  -1;
    default: 
        return -2;
    }
}

#define FFT_ERR_OK      0               // no error
#define FFT_ERR_LD      1               // len is not a power of 2
#define FFT_ERR_MAX     2               // len too large

typedef float   f_t;
typedef f_t     compl [2];
compl           root    [MAX >> 1];             // Sinus-/Kosinustabelle
size_t          shuffle [MAX >> 1] [2];         // Shuffle-Tabelle
size_t          shuffle_len;

// Bitinversion

size_t  swap ( size_t number, int bits )
{
    size_t  ret;
    for ( ret = 0; bits--; number >>= 1 ) {
        ret = ret + ret + (number & 1);
    }
    return ret;
}

// Bestimmen des Logarithmus dualis

int  ld ( size_t number )
{
    size_t i;
    for ( i = 0; i < sizeof(size_t)*CHAR_BIT; i++ )
        if ( ((size_t)1 << i) == number )
            return i;
    return -1;
}

// Die eigentliche FFT

int  fft ( compl* fn, const size_t newlen )
{
    static size_t  len  = 0;
    static int     bits = 0;
    size_t         i;
    size_t         j;
    size_t         k;
    size_t         p;

    /* Tabellen initialisieren */

    if ( newlen != len ) {
        len  = newlen;

        if ( (bits=ld(len)) == -1 )
            return FFT_ERR_LD;

        for ( i = 0; i < len; i++ ) {
            j = swap ( i, bits );
            if ( i < j ) {
                shuffle [shuffle_len] [0] = i;
                shuffle [shuffle_len] [1] = j;
                shuffle_len++;
            }
        }
        for ( i = 0; i < (len>>1); i++ ) {
            double x = (double) swap ( i+i, bits ) * 2*M_PI/len;
            root [i] [0] = cos (x);
            root [i] [1] = sin (x);
        }
    }

    /* Eigentliche Transformation */

    p = len >> 1;
    do {
        f_t*  bp = (f_t*) root;
        f_t*  si = (f_t*) fn;
        f_t*  di = (f_t*) fn+p+p;

        do {
            k = p;
            do {
                f_t  mulr = bp[0]*di[0] - bp[1]*di[1];
                f_t  muli = bp[1]*di[0] + bp[0]*di[1];

                di[0]  = si[0] - mulr;
                di[1]  = si[1] - muli;
                si[0] += mulr;
                si[1] += muli;

                si += 2, di += 2;
            } while ( --k );
            si += p+p, di += p+p, bp += 2;
        } while ( si < &fn[len][0] );
    } while (p >>= 1);

    /* Bitinversion */

    for ( k = 0; k < shuffle_len; k++ ) {
        f_t  tmp;
        i   = shuffle [k] [0];
        j   = shuffle [k] [1];
        tmp = fn [i][0]; fn [i][0] = fn [j][0]; fn [j][0] = tmp;
        tmp = fn [i][1]; fn [i][1] = fn [j][1]; fn [j][1] = tmp;
    }

    return FFT_ERR_OK;
}

void  printnumber ( long double x )
{
    unsigned  exp = 0;
    
    if      ( x < 9.999995  ) fprintf ( stderr, "%7.5f",  (double)x );
    else if ( x < 99.99995  ) fprintf ( stderr, "%7.4f",  (double)x );
    else if ( x < 999.9995  ) fprintf ( stderr, "%7.3f",  (double)x );
    else if ( x < 9999.995  ) fprintf ( stderr, "%7.2f",  (double)x );
    else if ( x < 99999.95  ) fprintf ( stderr, "%7.1f",  (double)x );
    else if ( x < 999999.5  ) fprintf ( stderr, "%6.0f.", (double)x );
    else if ( x < 9999999.5 ) fprintf ( stderr, "%7.0f",  (double)x );
    else if ( x < 9.9995e9  ) {
        while ( x >= 9.9995 ) exp++  , x /= 10;
        fprintf ( stderr, "%5.3fe%01u", (double)x, exp );
    } else if ( x < 9.995e99 ) {
        while ( x >= 9.5e6  ) exp+=6 , x /= 1.e6;
        while ( x >= 9.995  ) exp++  , x /= 10;
        fprintf ( stderr, "%4.2fe%02u", (double)x, exp );
    } else if ( x < 9.95e999L ) {
        while ( x >= 9.5e18 ) exp+=18, x /= 1.e18;
        while ( x >= 9.95   ) exp++  , x /= 10;
        fprintf ( stderr, "%3.1fe%03u", (double)x, exp );
    } else {
        while ( x >= 9.5e48 ) exp+=48, x /= 1.e48;
        while ( x >= 9.5    ) exp++  , x /= 10;
        fprintf ( stderr, "%1.0f.e%04u", (double)x, exp );
    }
}

double logdual ( long double x )
{
    unsigned exp = 0;
    
    while ( x >= 18446744073709551616. )
        x /= 18446744073709551616., exp += 64;
    while ( x >= 256. )
        x /= 256., exp += 8;
    while ( x >= 2. )
        x /= 2., exp += 1;
    return exp + log (x)/log(2);
}

int  random_number ( void )
{
    struct timeval  t;
    unsigned long   val;

    gettimeofday ( &t, NULL );

    val  = t.tv_sec ^ t.tv_usec ^ rand();
    val ^= val >> 16;
    val ^= val >>  8;
    val ^= val >>  4;
    val ^= val >>  2;
    val ^= val >>  1;

    return val & 1;
}

long double  prob ( int last, int total )
{
    long double  sum = 0.;
    long double  tmp = 1.;
    int          i;
    int          j   = total;
    
    if ( 2*last == total )
        return 1.;
    if ( 2*last > total )
        last = total - last;
        
    for ( i = 0; i <= last; i++ ) {
        sum += tmp;
        tmp  = tmp * (total-i) / (1+i);
        while ( j > 0  &&  tmp > 1 )
            j--, sum *= 0.5, tmp *= 0.5;
    }
    while ( j > 0 )
        j--, sum *= 0.5;

    return 2.*sum;
}


void  eval ( int right )
{
    static int   count = 0;
    static int   okay  = 0;
    long double  val;
    
    count ++;
    okay  += right;
    
    val    = 1.L / prob ( okay, count );
    
    fprintf (stderr, "   %s %5u/%-5u ", right ? "OK" : "- " , okay, count );
    printnumber (val);
    if ( count > 1 )
        fprintf (stderr, "   %4.2f bit", 0.01 * (int)(logdual(val) / (count-1) * 100.) );
    fprintf ( stderr, "\n" );
}


typedef signed short  sample_t;
typedef sample_t      mono_t   [1];
typedef sample_t      stereo_t [2];
typedef struct {
    unsigned long       n;
    long double         x;
    long double         x2;
    long double         y;
    long double         y2;
    long double         xy;
} korr_t;


void  analyze_stereo ( const stereo_t* p1, const stereo_t* p2, size_t len, korr_t* const k )
{
    long double  _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
    double       t1;
    double       t2;
    
    k -> n  += 2*len;
    
    for ( ; len--; p1++, p2++ ) {
        _x  += (t1 = (*p1)[0]); _x2 += t1 * t1;
        _y  += (t2 = (*p2)[0]); _y2 += t2 * t2;
                                _xy += t1 * t2;
        _x  += (t1 = (*p1)[1]); _x2 += t1 * t1;
        _y  += (t2 = (*p2)[1]); _y2 += t2 * t2;
                                _xy += t1 * t2;
    }
    
    k -> x  += _x ;
    k -> x2 += _x2;
    k -> y  += _y ;
    k -> y2 += _y2;
    k -> xy += _xy;
}

int  sgn ( double x )
{
    if ( x == 0 ) return 0;
    if ( x <  0 ) return -1;
    return +1;
}    

long double  report ( const korr_t* const k )
{
    long double  r;
    long double  sx;
    long double  sy;
    long double  x;
    long double  y;
    long double  b;
    
    r  = (k->x2*k->n - k->x*k->x) * (k->y2*k->n - k->y*k->y);
    r  = r  > 0.l  ?  (k->xy*k->n - k->x*k->y) / sqrt (r)  :  1.l;