1 /*
2 * Usage: abx original_file test_file
3 *
4 * Ask you as long as the probability is below the given percentage that
5 * you recognize differences
6 *
7 * Example: abx music.wav music.mp3
8 * abx music.wav music.mp3 --help
9 *
10 * Note: several 'decoding' utilites must be on the 'right' place
11 *
12 * Bugs:
13 * fix path of decoding utilities
14 * only 16 bit support
15 * only support of the same sample frequency
16 * no exact WAV file header analysis
17 * no mouse or joystick support
18 * don't uses functionality of ath.c
19 * only 2 files are comparable
20 * worse user interface
21 * quick & dirty hack
22 * wastes memory
23 * compile time warnings
24 * buffer overruns possible
25 * no dithering if recalcs are necessary
26 * correlation only done with one channel (2 channels, sum, what is better?)
27 * lowpass+highpass filtering (300 Hz+2*5 kHz) before delay+amplitude corr
28 * cross fade at start/stop
29 * non portable keyboard
30 * fade out on quit, fade in on start
31 * level/delay ajustment should be switchable
32 * pause key missing
33 * problems with digital silence files (division by 0)
34 * Gr��e cross corr fenster 2^16...18
35 * Stellensuche, ab 0*len oder 0.1*len oder 0.25*len, nach Effektiv oder Spitzenwert
36 * Absturz bei LPAC feeding, warum?
37 * Als 'B' beim Ratespiel sollte auch '0'...'9' verwendbar sein
38 * Oder mit einem Filter 300 Hz...3 kHz vorher filtern?
39 * Multiple encoded differenziertes Signal
40 * Amplitudenanpassung schaltbar machen?
41 * Direkt auf der Kommandozeile kodieren:
42 * abx "test.wav" "!lame -b128 test.wav -"
43 */
44
45 // If the program should increase it priority while playing define USE_NICE.
46 // Program must be installed SUID root. Decompressing phase is using NORMAL priority
47 #define USE_NICE
48
49 // Not only increase priority but change to relatime scheduling. Program must be installed SUID root
50 #define USE_REALTIME
51
52 // Path of the programs: mpg123, mppdec, faad, ac3dec, ogg123, lpac, shorten, MAC, flac
53 //#define PATH_OF_EXTERNAL_TOOLS_FOR_UNCOMPRESSING "/usr/local/bin/"
54 #define PATH_OF_EXTERNAL_TOOLS_FOR_UNCOMPRESSING ""
55
56
57 #if defined HAVE_CONFIG_H
58 # include <config.h>
59 #endif
60
61 #include <assert.h>
62 #include <ctype.h>
63 #include <fcntl.h>
64 #include <limits.h>
65 #include <math.h>
66 #include <memory.h>
67 #include <signal.h>
68 #include <stdio.h>
69 #include <stdlib.h>
70 #include <string.h>
71 #include <termios.h>
72 #include <time.h>
73 #include <unistd.h>
74 #include <sys/ioctl.h>
75 #include <sys/mman.h>
76 #include <sys/stat.h>
77 #include <sys/time.h>
78 #include <sys/types.h>
79
80 #define MAX (1<<17)
81
82 #if defined HAVE_SYS_SOUNDCARD_H
83 # include <sys/soundcard.h>
84 #elif defined HAVE_LINUX_SOUNDCARD_H
85 # include <linux/soundcard.h>
86 #else
87 # include <linux/soundcard.h> /* stand alone compilable for my tests */
88 #endif
89
90 #if defined USE_NICE
91 # include <sys/resource.h>
92 #endif
93 #if defined USE_REALTIME
94 # include <sched.h>
95 #endif
96
97 #define BF ((freq)/25)
98 #define MAX_LEN (210 * 44100)
99 #define DMA_SAMPLES 512 /* My Linux driver uses a DMA buffer of 65536*16 bit, which is 32768 samples in 16 bit stereo mode */
100
Set_Realtime(void)101 void Set_Realtime ( void )
102 {
103 #if defined USE_REALTIME
104 struct sched_param sp;
105 int ret;
106
107 memset ( &sp, 0, sizeof(sp) );
108 seteuid ( 0 );
109 sp.sched_priority = sched_get_priority_min ( SCHED_FIFO );
110 ret = sched_setscheduler ( 0, SCHED_RR, &sp );
111 seteuid ( getuid() );
112 #endif
113
114 #if defined USE_NICE
115 seteuid ( 0 );
116 setpriority ( PRIO_PROCESS, getpid(), -20 );
117 seteuid ( getuid() );
118 #endif
119 }
120
121 int verbose = 0;
122
123 static struct termios stored_settings;
124
125
reset(void)126 void reset ( void )
127 {
128 tcsetattr ( 0, TCSANOW, &stored_settings );
129 }
130
131
set(void)132 void set ( void )
133 {
134 struct termios new_settings;
135
136 tcgetattr ( 0, &stored_settings );
137 new_settings = stored_settings;
138
139 new_settings.c_lflag &= ~ECHO;
140 /* Disable canonical mode, and set buffer size to 1 byte */
141 new_settings.c_lflag &= ~ICANON;
142 new_settings.c_cc[VTIME] = 0;
143 new_settings.c_cc[VMIN] = 1;
144
145 tcsetattr(0,TCSANOW,&new_settings);
146 return;
147 }
148
149
sel(void)150 int sel ( void )
151 {
152 struct timeval t;
153 fd_set fd [1];
154 int ret;
155 unsigned char c;
156
157 FD_SET (0, fd);
158 t.tv_sec = 0;
159 t.tv_usec = 0;
160
161 ret = select ( 1, fd, NULL, NULL, &t );
162
163 switch ( ret ) {
164 case 0:
165 return -1;
166 case 1:
167 ret = read (0, &c, 1);
168 return ret == 1 ? c : -1;
169 default:
170 return -2;
171 }
172 }
173
174 #define FFT_ERR_OK 0 // no error
175 #define FFT_ERR_LD 1 // len is not a power of 2
176 #define FFT_ERR_MAX 2 // len too large
177
178 typedef float f_t;
179 typedef f_t compl [2];
180 compl root [MAX >> 1]; // Sinus-/Kosinustabelle
181 size_t shuffle [MAX >> 1] [2]; // Shuffle-Tabelle
182 size_t shuffle_len;
183
184 // Bitinversion
185
swap(size_t number,int bits)186 size_t swap ( size_t number, int bits )
187 {
188 size_t ret;
189 for ( ret = 0; bits--; number >>= 1 ) {
190 ret = ret + ret + (number & 1);
191 }
192 return ret;
193 }
194
195 // Bestimmen des Logarithmus dualis
196
ld(size_t number)197 int ld ( size_t number )
198 {
199 size_t i;
200 for ( i = 0; i < sizeof(size_t)*CHAR_BIT; i++ )
201 if ( ((size_t)1 << i) == number )
202 return i;
203 return -1;
204 }
205
206 // Die eigentliche FFT
207
fft(compl * fn,const size_t newlen)208 int fft ( compl* fn, const size_t newlen )
209 {
210 static size_t len = 0;
211 static int bits = 0;
212 size_t i;
213 size_t j;
214 size_t k;
215 size_t p;
216
217 /* Tabellen initialisieren */
218
219 if ( newlen != len ) {
220 len = newlen;
221
222 if ( (bits=ld(len)) == -1 )
223 return FFT_ERR_LD;
224
225 for ( i = 0; i < len; i++ ) {
226 j = swap ( i, bits );
227 if ( i < j ) {
228 shuffle [shuffle_len] [0] = i;
229 shuffle [shuffle_len] [1] = j;
230 shuffle_len++;
231 }
232 }
233 for ( i = 0; i < (len>>1); i++ ) {
234 double x = (double) swap ( i+i, bits ) * 2*M_PI/len;
235 root [i] [0] = cos (x);
236 root [i] [1] = sin (x);
237 }
238 }
239
240 /* Eigentliche Transformation */
241
242 p = len >> 1;
243 do {
244 f_t* bp = (f_t*) root;
245 f_t* si = (f_t*) fn;
246 f_t* di = (f_t*) fn+p+p;
247
248 do {
249 k = p;
250 do {
251 f_t mulr = bp[0]*di[0] - bp[1]*di[1];
252 f_t muli = bp[1]*di[0] + bp[0]*di[1];
253
254 di[0] = si[0] - mulr;
255 di[1] = si[1] - muli;
256 si[0] += mulr;
257 si[1] += muli;
258
259 si += 2, di += 2;
260 } while ( --k );
261 si += p+p, di += p+p, bp += 2;
262 } while ( si < &fn[len][0] );
263 } while (p >>= 1);
264
265 /* Bitinversion */
266
267 for ( k = 0; k < shuffle_len; k++ ) {
268 f_t tmp;
269 i = shuffle [k] [0];
270 j = shuffle [k] [1];
271 tmp = fn [i][0]; fn [i][0] = fn [j][0]; fn [j][0] = tmp;
272 tmp = fn [i][1]; fn [i][1] = fn [j][1]; fn [j][1] = tmp;
273 }
274
275 return FFT_ERR_OK;
276 }
277
printnumber(long double x)278 void printnumber ( long double x )
279 {
280 unsigned exp = 0;
281
282 if ( x < 9.999995 ) fprintf ( stderr, "%7.5f", (double)x );
283 else if ( x < 99.99995 ) fprintf ( stderr, "%7.4f", (double)x );
284 else if ( x < 999.9995 ) fprintf ( stderr, "%7.3f", (double)x );
285 else if ( x < 9999.995 ) fprintf ( stderr, "%7.2f", (double)x );
286 else if ( x < 99999.95 ) fprintf ( stderr, "%7.1f", (double)x );
287 else if ( x < 999999.5 ) fprintf ( stderr, "%6.0f.", (double)x );
288 else if ( x < 9999999.5 ) fprintf ( stderr, "%7.0f", (double)x );
289 else if ( x < 9.9995e9 ) {
290 while ( x >= 9.9995 ) exp++ , x /= 10;
291 fprintf ( stderr, "%5.3fe%01u", (double)x, exp );
292 } else if ( x < 9.995e99 ) {
293 while ( x >= 9.5e6 ) exp+=6 , x /= 1.e6;
294 while ( x >= 9.995 ) exp++ , x /= 10;
295 fprintf ( stderr, "%4.2fe%02u", (double)x, exp );
296 } else if ( x < 9.95e999L ) {
297 while ( x >= 9.5e18 ) exp+=18, x /= 1.e18;
298 while ( x >= 9.95 ) exp++ , x /= 10;
299 fprintf ( stderr, "%3.1fe%03u", (double)x, exp );
300 } else {
301 while ( x >= 9.5e48 ) exp+=48, x /= 1.e48;
302 while ( x >= 9.5 ) exp++ , x /= 10;
303 fprintf ( stderr, "%1.0f.e%04u", (double)x, exp );
304 }
305 }
306
logdual(long double x)307 double logdual ( long double x )
308 {
309 unsigned exp = 0;
310
311 while ( x >= 18446744073709551616. )
312 x /= 18446744073709551616., exp += 64;
313 while ( x >= 256. )
314 x /= 256., exp += 8;
315 while ( x >= 2. )
316 x /= 2., exp += 1;
317 return exp + log (x)/log(2);
318 }
319
random_number(void)320 int random_number ( void )
321 {
322 struct timeval t;
323 unsigned long val;
324
325 gettimeofday ( &t, NULL );
326
327 val = t.tv_sec ^ t.tv_usec ^ rand();
328 val ^= val >> 16;
329 val ^= val >> 8;
330 val ^= val >> 4;
331 val ^= val >> 2;
332 val ^= val >> 1;
333
334 return val & 1;
335 }
336
prob(int last,int total)337 long double prob ( int last, int total )
338 {
339 long double sum = 0.;
340 long double tmp = 1.;
341 int i;
342 int j = total;
343
344 if ( 2*last == total )
345 return 1.;
346 if ( 2*last > total )
347 last = total - last;
348
349 for ( i = 0; i <= last; i++ ) {
350 sum += tmp;
351 tmp = tmp * (total-i) / (1+i);
352 while ( j > 0 && tmp > 1 )
353 j--, sum *= 0.5, tmp *= 0.5;
354 }
355 while ( j > 0 )
356 j--, sum *= 0.5;
357
358 return 2.*sum;
359 }
360
361
eval(int right)362 void eval ( int right )
363 {
364 static int count = 0;
365 static int okay = 0;
366 long double val;
367
368 count ++;
369 okay += right;
370
371 val = 1.L / prob ( okay, count );
372
373 fprintf (stderr, " %s %5u/%-5u ", right ? "OK" : "- " , okay, count );
374 printnumber (val);
375 if ( count > 1 )
376 fprintf (stderr, " %4.2f bit", 0.01 * (int)(logdual(val) / (count-1) * 100.) );
377 fprintf ( stderr, "\n" );
378 }
379
380
381 typedef signed short sample_t;
382 typedef sample_t mono_t [1];
383 typedef sample_t stereo_t [2];
384 typedef struct {
385 unsigned long n;
386 long double x;
387 long double x2;
388 long double y;
389 long double y2;
390 long double xy;
391 } korr_t;
392
393
analyze_stereo(const stereo_t * p1,const stereo_t * p2,size_t len,korr_t * const k)394 void analyze_stereo ( const stereo_t* p1, const stereo_t* p2, size_t len, korr_t* const k )
395 {
396 long double _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
397 double t1;
398 double t2;
399
400 k -> n += 2*len;
401
402 for ( ; len--; p1++, p2++ ) {
403 _x += (t1 = (*p1)[0]); _x2 += t1 * t1;
404 _y += (t2 = (*p2)[0]); _y2 += t2 * t2;
405 _xy += t1 * t2;
406 _x += (t1 = (*p1)[1]); _x2 += t1 * t1;
407 _y += (t2 = (*p2)[1]); _y2 += t2 * t2;
408 _xy += t1 * t2;
409 }
410
411 k -> x += _x ;
412 k -> x2 += _x2;
413 k -> y += _y ;
414 k -> y2 += _y2;
415 k -> xy += _xy;
416 }
417
sgn(double x)418 int sgn ( double x )
419 {
420 if ( x == 0 ) return 0;
421 if ( x < 0 ) return -1;
422 return +1;
423 }
424
report(const korr_t * const k)425 long double report ( const korr_t* const k )
426 {
427 long double r;
428 long double sx;
429 long double sy;
430 long double x;
431 long double y;
432 long double b;
433
434 r = (k->x2*k->n - k->x*k->x) * (k->y2*k->n - k->y*k->y);
435 r = r > 0.l ? (k->xy*k->n - k->x*k->y) / sqrt (r) : 1.l;
436 sx = k->n > 1 ? sqrt ( (k->x2 - k->x*k->x/k->n) / (k->n - 1) ) : 0.l;
437 sy = k->n > 1 ? sqrt ( (k->y2 - k->y*k->y/k->n) / (k->n - 1) ) : 0.l;
438 x = k->n > 0 ? k->x/k->n : 0.l;
439 y = k->n > 0 ? k->y/k->n : 0.l;
440
441 b = sx != 0 ? sy/sx * sgn(r) : 0.l;
442 if (verbose)
443 fprintf ( stderr, "r=%Lf sx=%Lf sy=%Lf x=%Lf y=%Lf b=%Lf\n", r, sx, sy, x, y, b );
444 return b;
445 }
446
447
448 /* Input: an unsigned short n.
449 * Output: the swapped bytes of n if the arch is big-endian or n itself
450 * if the arch is little-endian.
451 * Comment: should be replaced latter with a better solution than this
452 * home-brewed hack (rbrito). The name should be better also.
453 */
be16_le(unsigned short n)454 inline unsigned short be16_le(unsigned short n)
455 {
456 #ifdef _WORDS_BIGENDIAN
457 return (n << 8) | (n >> 8);
458 #else
459 return n;
460 #endif
461 }
462
463
feed(int fd,const stereo_t * p,int len)464 int feed ( int fd, const stereo_t* p, int len )
465 {
466 int i;
467 stereo_t tmp[30000]; /* An arbitrary size--to be changed latter */
468
469 if (len > sizeof(tmp)/sizeof(*tmp))
470 len = sizeof(tmp)/sizeof(*tmp);
471
472 for (i = 0; i < len; i++) {
473 tmp[i][0] = be16_le(p[i][0]);
474 tmp[i][1] = be16_le(p[i][1]);
475 }
476
477 write ( fd, tmp, sizeof(stereo_t) * len );
478 return len;
479 }
480
481
round(double f)482 short round ( double f )
483 {
484 long x = (long) floor ( f + 0.5 );
485 return x == (short)x ? (short)x : (short) ((x >> 31) ^ 0x7FFF);
486 }
487
488
feed2(int fd,const stereo_t * p1,const stereo_t * p2,int len)489 int feed2 ( int fd, const stereo_t* p1, const stereo_t* p2, int len )
490 {
491 stereo_t tmp [30000]; /* An arbitrary size, hope that no overruns occure */
492 int i;
493
494 if (len > sizeof(tmp)/sizeof(*tmp))
495 len = sizeof(tmp)/sizeof(*tmp);
496 for ( i = 0; i < len; i++ ) {
497 double f = cos ( M_PI/2*i/len );
498 f *= f;
499 tmp [i] [0] = be16_le(round ( p1 [i] [0] * f + p2 [i] [0] * (1. - f) ));
500 tmp [i] [1] = be16_le(round ( p1 [i] [1] * f + p2 [i] [1] * (1. - f) ));
501 }
502
503 write ( fd, tmp, sizeof(stereo_t) * len );
504 return len;
505 }
506
507
feedfac(int fd,const stereo_t * p1,const stereo_t * p2,int len,double fac1,double fac2)508 int feedfac ( int fd, const stereo_t* p1, const stereo_t* p2, int len, double fac1, double fac2 )
509 {
510 stereo_t tmp [30000]; /* An arbitrary size, hope that no overruns occure */
511 int i;
512
513 if (len > sizeof(tmp)/sizeof(*tmp))
514 len = sizeof(tmp)/sizeof(*tmp);
515 for ( i = 0; i < len; i++ ) {
516 tmp [i] [0] = be16_le(round ( p1 [i] [0] * fac1 + p2 [i] [0] * fac2 ));
517 tmp [i] [1] = be16_le(round ( p1 [i] [1] * fac1 + p2 [i] [1] * fac2 ));
518 }
519
520 write ( fd, tmp, sizeof(stereo_t) * len );
521 return len;
522 }
523
524
setup(int fdd,int samples,long freq)525 void setup ( int fdd, int samples, long freq )
526 {
527 int status, org, arg;
528
529 // Nach vorn verschoben
530 if ( -1 == (status = ioctl (fdd, SOUND_PCM_SYNC, 0)) )
531 perror ("SOUND_PCM_SYNC ioctl failed");
532
533 org = arg = 2;
534 if ( -1 == (status = ioctl (fdd, SOUND_PCM_WRITE_CHANNELS, &arg)) )
535 perror ("SOUND_PCM_WRITE_CHANNELS ioctl failed");
536 if (arg != org)
537 perror ("unable to set number of channels");
538 fprintf (stderr, "%1u*", arg);
539
540 org = arg = AFMT_S16_LE;
541 if ( -1 == ioctl (fdd, SNDCTL_DSP_SETFMT, &arg) )
542 perror ("SNDCTL_DSP_SETFMT ioctl failed");
543 if ((arg & org) == 0)
544 perror ("unable to set data format");
545
546 org = arg = freq;
547 if ( -1 == (status = ioctl (fdd, SNDCTL_DSP_SPEED, &arg)) )
548 perror ("SNDCTL_DSP_SPEED ioctl failed");
549 fprintf (stderr, "%5u Hz*%.3f sec\n", arg, (double)samples/arg );
550
551 }
552
553
Message(const char * s,size_t index,long freq,size_t start,size_t stop)554 void Message ( const char* s, size_t index, long freq, size_t start, size_t stop )
555 {
556 unsigned long norm_index = 100lu * index / freq;
557 unsigned long norm_start = 100lu * start / freq;
558 unsigned long norm_stop = 100lu * stop / freq;
559
560 fprintf ( stderr, "\rListening %s %2lu:%02lu.%02lu (%1lu:%02lu.%02lu...%1lu:%02lu.%02lu)%*.*s\rListening %s",
561 s,
562 norm_index / 6000, norm_index / 100 % 60, norm_index % 100,
563 norm_start / 6000, norm_start / 100 % 60, norm_start % 100,
564 norm_stop / 6000, norm_stop / 100 % 60, norm_stop % 100,
565 36 - (int)strlen(s), 36 - (int)strlen(s), "",
566 s );
567
568 fflush ( stderr );
569 }
570
571
calc_true_index(size_t index,size_t start,size_t stop)572 size_t calc_true_index ( size_t index, size_t start, size_t stop )
573 {
574 if ( start >= stop )
575 return start;
576 while ( index - start < DMA_SAMPLES )
577 index += stop - start;
578 return index - DMA_SAMPLES;
579 }
580
581
testing(const stereo_t * A,const stereo_t * B,size_t len,long freq)582 void testing ( const stereo_t* A, const stereo_t* B, size_t len, long freq )
583 {
584 int c;
585 int fd = open ( "/dev/dsp", O_WRONLY );
586 int rnd = random_number (); /* Auswahl von X */
587 int state = 0; /* derzeitiger F�ttungsmodus */
588 float fac1 = 0.5;
589 float fac2 = 0.5;
590 size_t start = 0;
591 size_t stop = len;
592 size_t index = start; /* derzeitiger Offset auf den Audiostr�men */
593 char message [80] = "A ";
594
595 setup ( fd, len, freq );
596
597 while ( 1 ) {
598 c = sel ();
599 if ( c == 27 )
600 c = sel () + 0x100;
601
602 switch ( c ) {
603 case 'A' :
604 case 'a' :
605 strcpy ( message, "A " );
606 if ( state != 0 )
607 state = 2;
608 break;
609
610 case 0x100+'0' :
611 case '0' :
612 case 'B' :
613 case 'b' :
614 strcpy ( message, " B" );
615 if ( state != 1 )
616 state = 3;
617 break;
618
619 case 'X' :
620 case 'x' :
621 strcpy ( message, " X " );
622 if ( state != rnd )
623 state = rnd + 2;
624 break;
625
626 case 'm' :
627 state = 8;
628 break;
629
630 case 'M' :
631 state = (state & 1) + 4;
632 break;
633
634 case 'x'&0x1F:
635 state = (state & 1) + 6;
636 break;
637
638 case ' ':
639 start = 0;
640 stop = len;
641 break;
642
643 case 'o' :
644 start = calc_true_index ( index, start, stop);
645 break;
646 case 'p' :
647 stop = calc_true_index ( index, start, stop);
648 break;
649 case 'h' :
650 if ( start > freq/100 )
651 start -= freq/100;
652 else
653 start = 0;
654 index = start;
655 continue;
656 case 'j' :
657 if ( start < stop-freq/100 )
658 start += freq/100;
659 else
660 start = stop;
661 index = start;
662 continue;
663 case 'k' :
664 if ( stop > start+freq/100 )
665 stop -= freq/100;
666 else
667 stop = start;
668 continue;
669 case 'l' :
670 if ( stop < len-freq/100 )
671 stop += freq/100;
672 else
673 stop = len;
674 continue;
675 case '\n':
676 index = start;
677 continue;
678
679 case 'D'+0x100:
680 strcpy ( message, "Difference (+40 dB)" );
681 state = 9;
682 fac1 = -100.;
683 fac2 = +100.;
684 break;
685
686 case 'd'+0x100:
687 strcpy ( message, "Difference (+30 dB)" );
688 state = 9;
689 fac1 = -32.;
690 fac2 = +32.;
691 break;
692
693 case 'D' & 0x1F :
694 strcpy ( message, "Difference (+20 dB)" );
695 state = 9;
696 fac1 = -10.;
697 fac2 = +10.;
698 break;
699
700 case 'D' :
701 strcpy ( message, "Difference (+10 dB)" );
702 state = 9;
703 fac1 = -3.;
704 fac2 = +3.;
705 break;
706
707 case 'd' :
708 strcpy ( message, "Difference ( 0 dB)" );
709 state = 9;
710 fac1 = -1.;
711 fac2 = +1.;
712 break;
713
714 case 0x100+'1' :
715 case 0x100+'2' :
716 case 0x100+'3' :
717 case 0x100+'4' :
718 case 0x100+'5' :
719 case 0x100+'6' :
720 case 0x100+'7' :
721 case 0x100+'8' :
722 case 0x100+'9' :
723 sprintf ( message, " B (Errors -%c dB)", (char)c );
724 state = 9;
725 fac2 = pow (10., -0.05*(c-0x100-'0') );
726 fac1 = 1. - fac2;
727 break;
728
729 case '1' :
730 case '2' :
731 case '3' :
732 case '4' :
733 case '5' :
734 case '6' :
735 case '7' :
736 case '8' :
737 case '9' :
738 sprintf ( message, " B (Errors +%c dB)", c );
739 state = 9;
740 fac2 = pow (10., 0.05*(c-'0') );
741 fac1 = 1. - fac2;
742 break;
743
744 case 'A' & 0x1F:
745 fprintf (stderr, " Vote for X:=A" );
746 eval ( rnd == 0 );
747 rnd = random_number ();
748 if ( state == 6 && state == 7 )
749 state = 6 + rnd;
750 else if ( state != rnd )
751 state = rnd + 2;
752 strcpy ( message," X " );
753 break;
754
755 case 'B' & 0x1F:
756 fprintf (stderr, " Vote for X:=B" );
757 eval ( rnd == 1 );
758 rnd = random_number ();
759 if ( state == 6 && state == 7 )
760 state = 6 + rnd;
761 else if ( state != rnd )
762 state = rnd + 2;
763 strcpy ( message," X " );
764 break;
765
766 case -1:
767 break;
768
769 default:
770 fprintf (stderr, "\a" );
771 break;
772
773 case 'Q':
774 case 'q':
775 fprintf ( stderr, "\n%-79.79s\r", "Quit program" );
776 close (fd);
777 fprintf ( stderr, "\n\n");
778 return;
779 }
780
781 switch (state) {
782 case 0: /* A */
783 if ( index + BF >= stop )
784 index += feed (fd, A+index, stop-index );
785 else
786 index += feed (fd, A+index, BF );
787 break;
788
789 case 1: /* B */
790 if ( index + BF >= stop )
791 index += feed (fd, B+index, stop-index );
792 else
793 index += feed (fd, B+index, BF );
794 break;
795
796 case 2: /* B => A */
797 if ( index + BF >= stop )
798 index += feed2 (fd, B+index, A+index, stop-index );
799 else
800 index += feed2 (fd, B+index, A+index, BF );
801 state = 0;
802 break;
803
804 case 3: /* A => B */
805 if ( index + BF >= stop )
806 index += feed2 (fd, A+index, B+index, stop-index );
807 else
808 index += feed2 (fd, A+index, B+index, BF );
809 state = 1;
810 break;
811
812 case 4: /* A */
813 strcpy ( message, "A " );
814 if ( index + BF >= stop )
815 index += feed (fd, A+index, stop-index ),
816 state++;
817 else
818 index += feed (fd, A+index, BF );
819 break;
820
821 case 5: /* B */
822 strcpy ( message, " B" );
823 if ( index + BF >= stop )
824 index += feed (fd, B+index, stop-index ),
825 state--;
826 else
827 index += feed (fd, B+index, BF );
828 break;
829
830 case 6: /* X */
831 strcpy ( message, " X " );
832 if ( index + BF >= stop )
833 index += feed (fd, (rnd ? B : A)+index, stop-index ),
834 state++;
835 else
836 index += feed (fd, (rnd ? B : A)+index, BF );
837 break;
838
839 case 7: /* !X */
840 strcpy ( message, "!X " );
841 if ( index + BF >= stop )
842 index += feed (fd, (rnd ? A : B)+index, stop-index ),
843 state--;
844 else
845 index += feed (fd, (rnd ? A : B)+index, BF );
846 break;
847
848 case 8:
849 if ( index + BF/2 >= stop )
850 index += feed2 (fd, A+index, B+index, stop-index );
851 else
852 index += feed2 (fd, A+index, B+index, BF/2 );
853 Message ( " B", index, freq, start, stop );
854 if ( index + BF >= stop )
855 index += feed (fd, B+index, stop-index );
856 else
857 index += feed (fd, B+index, BF );
858 if ( index + BF/2 >= stop )
859 index += feed2 (fd, B+index, A+index, stop-index );
860 else
861 index += feed2 (fd, B+index, A+index, BF/2 );
862 Message ( "A ", index, freq, start, stop );
863 if ( index + BF >= stop )
864 index += feed (fd, A+index, stop-index );
865 else
866 index += feed (fd, A+index, BF );
867 break;
868
869 case 9: /* Liko */
870 if ( index + BF >= stop )
871 index += feedfac (fd, A+index, B+index, stop-index, fac1, fac2 );
872 else
873 index += feedfac (fd, A+index, B+index, BF , fac1, fac2 );
874 break;
875
876 default:
877 assert (0);
878 }
879
880 if (index >= stop)
881 index = start;
882 Message ( message, calc_true_index ( index, start, stop), freq, start, stop );
883 }
884 }
885
886
has_ext(const char * name,const char * ext)887 int has_ext ( const char* name, const char* ext )
888 {
889 if ( strlen (name) < strlen (ext) )
890 return 0;
891 name += strlen (name) - strlen (ext);
892 return strcasecmp (name, ext) ? 0 : 1;
893 }
894
895
896 typedef struct {
897 const char* const extention;
898 const char* const command;
899 } decoder_t;
900
901
902 #define REDIR " 2> /dev/null"
903 #define STDOUT "/dev/fd/1"
904 #define PATH PATH_OF_EXTERNAL_TOOLS_FOR_UNCOMPRESSING
905
906 const decoder_t decoder [] = {
907 { ".mp1" , PATH"mpg123 -w - %s" REDIR }, // MPEG Layer I : www.iis.fhg.de, www.mpeg.org
908 { ".mp2" , PATH"mpg123 -w - %s" REDIR }, // MPEG Layer II : www.iis.fhg.de, www.uq.net.au/~zzmcheng, www.mpeg.org
909 { ".mp3" , PATH"mpg123 -w - %s" REDIR }, // MPEG Layer III : www.iis.fhg.de, www.mp3dev.org, www.mpeg.org
910 { ".mp3pro" , PATH"mpg123 -w - %s" REDIR }, // MPEG Layer III : www.iis.fhg.de, www.mp3dev.org, www.mpeg.org
911 { ".mpt" , PATH"mpg123 -w - %s" REDIR }, // MPEG Layer III : www.iis.fhg.de, www.mp3dev.org, www.mpeg.org
912 { ".mpp" , PATH"mppdec %s -" REDIR }, // MPEGplus : www.stud.uni-hannover.de/user/73884
913 { ".mpc" , PATH"mppdec %s -" REDIR }, // MPEGplus : www.stud.uni-hannover.de/user/73884
914 { ".mp+" , PATH"mppdec %s -" REDIR }, // MPEGplus : www.stud.uni-hannover.de/user/73884
915 { ".aac" , PATH"faad -t.wav -w %s" REDIR }, // Advanced Audio Coding: psytel.hypermart.net, www.aac-tech.com, sourceforge.net/projects/faac, www.aac-audio.com, www.mpeg.org
916 { "aac.lqt" , PATH"faad -t.wav -w %s" REDIR }, // Advanced Audio Coding: psytel.hypermart.net, www.aac-tech.com, sourceforge.net/projects/faac, www.aac-audio.com, www.mpeg.org
917 { ".ac3" , PATH"ac3dec %s" REDIR }, // Dolby AC3 : www.att.com
918 { "ac3.lqt" , PATH"ac3dec %s" REDIR }, // Dolby AC3 : www.att.com
919 { ".ogg" , PATH"ogg123 -d wav -o file:"STDOUT" %s" REDIR }, // Ogg Vorbis : www.xiph.org/ogg/vorbis/index.html
920 { ".pac" , PATH"lpac -x %s "STDOUT REDIR }, // Lossless predictive Audio Compression: www-ft.ee.tu-berlin.de/~liebchen/lpac.html (liebchen@ft.ee.tu-berlin.de)
921 { ".shn" , PATH"shorten -x < %s" REDIR }, // Shorten : shnutils.freeshell.org, www.softsound.com/Shorten.html (shnutils@freeshell.org, shorten@softsound.com)
922 { ".wav.gz" , PATH"gzip -d < %s | sox -twav - -twav -sw -"REDIR }, // gziped WAV
923 { ".wav.sz" , PATH"szip -d < %s | sox -twav - -twav -sw -"REDIR }, // sziped WAV
924 { ".wav.sz2", PATH"szip2 -d < %s | sox -twav - -twav -sw -"REDIR }, // sziped WAV
925 { ".raw" , PATH"sox -r44100 -sw -c2 -traw %s -twav -sw -"REDIR }, // raw files are treated as CD like audio
926 { ".cdr" , PATH"sox -r44100 -sw -c2 -traw %s -twav -sw -"REDIR }, // CD-DA files are treated as CD like audio, no preemphasis info available
927 { ".rm" , "echo %s '???'" REDIR }, // Real Audio : www.real.com
928 { ".epc" , "echo %s '???'" REDIR }, // ePAC : www.audioveda.com, www.lucent.com/ldr
929 { ".mov" , "echo %s '???'" REDIR }, // QDesign Music 2 : www.qdesign.com
930 { ".vqf" , "echo %s '???'" REDIR }, // TwinVQ : www.yamaha-xg.com/english/xg/SoundVQ, www.vqf.com, sound.splab.ecl.ntt.co.jp/twinvq-e
931 { ".wma" , "echo %s '???'" REDIR }, // Microsoft Media Audio: www.windowsmedia.com, www.microsoft.com/windows/windowsmedia
932 { ".flac" , PATH"flac -c -d %s" REDIR }, // Free Lossless Audio Coder: flac.sourceforge.net/
933 { ".fla" , PATH"flac -c -d %s" REDIR }, // Free Lossless Audio Coder: flac.sourceforge.net/
934 { ".ape" , "( "PATH"MAC %s _._.wav -d > /dev/null; cat _._.wav; rm _._.wav )" REDIR }, // Monkey's Audio Codec : www.monkeysaudio.com (email@monkeysaudio.com)
935 { ".rka" , "( "PATH"rkau %s _._.wav > /dev/null; cat _._.wav; rm _._.wav )" REDIR }, // RK Audio:
936 { ".rkau" , "( "PATH"rkau %s _._.wav > /dev/null; cat _._.wav; rm _._.wav )" REDIR }, // RK Audio:
937 { ".mod" , PATH"xmp -b16 -c -f44100 --stereo -o- %s | sox -r44100 -sw -c2 -traw - -twav -sw -"
938 REDIR }, // Amiga's Music on Disk:
939 { "" , PATH"sox %s -twav -sw -" REDIR }, // Rest, may be sox can handle it
940 };
941
942 #undef REDIR
943 #undef STDOUT
944 #undef PATH
945
946
readwave(stereo_t * buff,size_t maxlen,const char * name,size_t * len)947 int readwave ( stereo_t* buff, size_t maxlen, const char* name, size_t* len )
948 {
949 char* command = malloc (2*strlen(name) + 512);
950 char* name_q = malloc (2*strlen(name) + 128);
951 unsigned short header [22];
952 FILE* fp;
953 size_t i;
954 size_t j;
955
956 // The *nice* shell quoting
957 i = j = 0;
958 if ( name[i] == '-' )
959 name_q[j++] = '.',
960 name_q[j++] = '/';
961
962 while (name[i]) {
963 if ( !isalnum (name[i]) && name[i]!='-' && name[i]!='_' && name[i]!='.' )
964 name_q[j++] = '\\';
965 name_q[j++] = name[i++];
966 }
967 name_q[j] = '\0';
968
969 fprintf (stderr, "Reading %s", name );
970 for ( i = 0; i < sizeof(decoder)/sizeof(*decoder); i++ )
971 if ( has_ext (name, decoder[i].extention) ) {
972 sprintf ( command, decoder[i].command, name_q );
973 break;
974 }
975
976 free (name_q);
977 if ( (fp = popen (command, "r")) == NULL ) {
978 fprintf (stderr, "Can't exec:\n%s\n", command );
979 exit (1);
980 }
981 free (command);
982
983 fprintf (stderr, " ..." );
984 fread ( header, sizeof(*header), sizeof(header)/sizeof(*header), fp );
985
986 switch (be16_le(header[11])) {
987 case 2:
988 *len = fread ( buff, sizeof(stereo_t), maxlen, fp );
989 for (i = 0; i < *len; i ++) {
990 buff[i][0] = be16_le(buff[i][0]);
991 buff[i][1] = be16_le(buff[i][1]);
992 }
993 break;
994 case 1:
995 *len = fread ( buff, sizeof(sample_t), maxlen, fp );
996 for ( i = *len; i-- > 0; )
997 buff[i][0] = buff[i][1] = ((sample_t*)buff) [i];
998 break;
999 case 0:
1000 fprintf (stderr, "\b\b\b\b, Standard Open Source Bug detected, try murksaround ..." );
1001 *len = fread ( buff, sizeof(stereo_t), maxlen, fp );
1002 header[11] = 2;
1003 header[12] = 65534; /* use that of the other channel */
1004 break;
1005 default:
1006 fprintf (stderr, "Only 1 or 2 channels are supported, not %u\n", header[11] );
1007 pclose (fp);
1008 return -1;
1009 }
1010 pclose ( fp );
1011 fprintf (stderr, "\n" );
1012 return be16_le(header[12]) ? be16_le(header[12]) : 65534;
1013 }
1014
1015
cross_analyze(const stereo_t * p1,const stereo_t * p2,size_t len)1016 double cross_analyze ( const stereo_t* p1, const stereo_t *p2, size_t len )
1017 {
1018 float P1 [MAX] [2];
1019 float P2 [MAX] [2];
1020 int i;
1021 int maxindex;
1022 double sum1;
1023 double sum2;
1024 double max;
1025 double y1;
1026 double y2;
1027 double y3;
1028 double yo;
1029 double xo;
1030 double tmp;
1031 double tmp1;
1032 double tmp2;
1033 int ret = 0;
1034 int cnt = 5;
1035
1036 // Calculating effective voltage
1037 sum1 = sum2 = 0.;
1038 for ( i = 0; i < len; i++ ) {
1039 sum1 += (double)p1[i][0] * p1[i][0];
1040 sum2 += (double)p2[i][0] * p2[i][0];
1041 }
1042 sum1 = sqrt ( sum1/len );
1043 sum2 = sqrt ( sum2/len );
1044
1045 // Searching beginning of signal (not stable for pathological signals)
1046 for ( i = 0; i < len; i++ )
1047 if ( abs (p1[i][0]) >= sum1 && abs (p2[i][0]) >= sum2 )
1048 break;
1049 p1 += i;
1050 p2 += i;
1051 len -= i;
1052
1053 if ( len <= MAX )
1054 return 0;
1055
1056 // Filling arrays for FFT
1057 do {
1058 sum1 = sum2 = 0.;
1059 for ( i = 0; i < MAX; i++ ) {
1060 #ifdef USEDIFF
1061 tmp1 = p1 [i][0] - p1 [i+1][0];
1062 tmp2 = p2 [i+ret][0] - p2 [i+ret+1][0];
1063 #else
1064 tmp1 = p1 [i][0];
1065 tmp2 = p2 [i+ret][0];
1066 #endif
1067 sum1 += tmp1*tmp1;
1068 sum2 += tmp2*tmp2;
1069 P1 [i][0] = tmp1;
1070 P2 [i][0] = tmp2;
1071 P1 [i][1] = 0.;
1072 P2 [i][1] = 0.;
1073 }
1074
1075 fft (P1, MAX);
1076 fft (P2, MAX);
1077
1078 for ( i = 0; i < MAX; i++ ) {
1079 double a0 = P1 [i][0];
1080 double a1 = P1 [i][1];
1081 double b0 = P2 [(MAX-i)&(MAX-1)][0];
1082 double b1 = P2 [(MAX-i)&(MAX-1)][1];
1083 P1 [i][0] = a0*b0 - a1*b1;
1084 P1 [i][1] = a0*b1 + a1*b0;
1085 }
1086
1087 fft (P1, MAX);
1088
1089 max = P1 [maxindex = 0][0];
1090 for ( i = 1; i < MAX; i++ )
1091 if ( P1[i][0] > max )
1092 max = P1 [maxindex = i][0];
1093
1094 y2 = P1 [ maxindex ][0];
1095 y1 = P1 [(maxindex-1)&(MAX-1)][0] - y2;
1096 y3 = P1 [(maxindex+1)&(MAX-1)][0] - y2;
1097
1098 xo = 0.5 * (y1-y3) / (y1+y3);
1099 yo = 0.5 * ( (y1+y3)*xo + (y3-y1) ) * xo;
1100
1101 if (maxindex > MAX/2 )
1102 maxindex -= MAX;
1103
1104 ret += maxindex;
1105 tmp = 100./MAX/sqrt(sum1*sum2);
1106 if (verbose)
1107 printf ( "[%5d]%8.4f [%5d]%8.4f [%5d]%8.4f [%10.4f]%8.4f\n",
1108 ret- 1, (y1+y2)*tmp,
1109 ret , y2 *tmp,
1110 ret+ 1, (y3+y2)*tmp,
1111 ret+xo, (yo+y2)*tmp );
1112
1113 } while ( maxindex && cnt-- );
1114
1115 return ret + xo;
1116 }
1117
1118
to_short(int x)1119 short to_short ( int x )
1120 {
1121 return x == (short)x ? (short)x : (short) ((x >> 31) ^ 0x7FFF);
1122 }
1123
1124
DC_cancel(stereo_t * p,size_t len)1125 void DC_cancel ( stereo_t* p, size_t len )
1126 {
1127 double sum1 = 0;
1128 double sum2 = 0;
1129 size_t i;
1130 int diff1;
1131 int diff2;
1132
1133 for (i = 0; i < len; i++ ) {
1134 sum1 += p[i][0];
1135 sum2 += p[i][1];
1136 }
1137 if ( fabs(sum1) < len && fabs(sum2) < len )
1138 return;
1139
1140 diff1 = round ( sum1 / len );
1141 diff2 = round ( sum2 / len );
1142 if (verbose)
1143 fprintf (stderr, "Removing DC (left=%d, right=%d)\n", diff1, diff2 );
1144
1145 for (i = 0; i < len; i++ ) {
1146 p[i][0] = to_short ( p[i][0] - diff1);
1147 p[i][1] = to_short ( p[i][1] - diff2);
1148 }
1149 }
1150
multiply(char c,stereo_t * p,size_t len,double fact)1151 void multiply ( char c, stereo_t* p, size_t len, double fact )
1152 {
1153 size_t i;
1154
1155 if ( fact == 1. )
1156 return;
1157 if (verbose)
1158 fprintf (stderr, "Multiplying %c by %7.5f\n", c, fact );
1159
1160 for (i = 0; i < len; i++ ) {
1161 p[i][0] = to_short ( p[i][0] * fact );
1162 p[i][1] = to_short ( p[i][1] * fact );
1163 }
1164 }
1165
1166
maximum(stereo_t * p,size_t len)1167 int maximum ( stereo_t* p, size_t len )
1168 {
1169 int max = 0;
1170 size_t i;
1171
1172 for (i = 0; i < len; i++ ) {
1173 if (abs(p[i][0]) > max) max = abs(p[i][0]);
1174 if (abs(p[i][1]) > max) max = abs(p[i][1]);
1175 }
1176 return max;
1177 }
1178
1179
usage(void)1180 void usage ( void )
1181 {
1182 fprintf ( stderr,
1183 "usage: abx [-v] File_A File_B\n"
1184 "\n"
1185 "File_A and File_B loaded and played. File_A should be the better/reference\n"
1186 "file, File_B the other. You can press the following keys:\n"
1187 "\n"
1188 " a/A: Listen to File A\n"
1189 " b/B: Listen to File B\n"
1190 " x/X: Listen to the randomly selected File X, which is A or B\n"
1191 " Ctrl-A: You vote for X=A\n"
1192 " Ctrl-B: You vote for X=B\n"
1193 " m: Alternating playing A and B. Fast switching\n"
1194 " M: Alternating playing A and B. Slow switching\n"
1195 " d/D/Ctrl-D/Alt-d/Alt-D:\n"
1196 " Listen to the difference A-B (+0 dB...+40 dB)\n"
1197 " o/p: Chunk select\n"
1198 " hjkl: Chunk fine adjust (hj: start, kl: stop)\n"
1199 " Space: Chunk deselect\n"
1200 " 0...9: Listen to B, but difference A-B is amplified by 0-9 dB\n"
1201 " Q: Quit the program\n"
1202 "\n"
1203 );
1204 }
1205
1206
main(int argc,char ** argv)1207 int main ( int argc, char** argv )
1208 {
1209 stereo_t* _A = calloc ( MAX_LEN, sizeof(stereo_t) );
1210 stereo_t* _B = calloc ( MAX_LEN, sizeof(stereo_t) );
1211 stereo_t* A = _A;
1212 stereo_t* B = _B;
1213 size_t len_A;
1214 size_t len_B;
1215 size_t len;
1216 int max_A;
1217 int max_B;
1218 int max;
1219 long freq1;
1220 long freq2;
1221 int shift;
1222 double fshift;
1223 double ampl;
1224 int ampl_X;
1225 korr_t k;
1226
1227 if (argc > 1 && 0 == strcmp (argv[1], "-v") ) {
1228 verbose = 1;
1229 argc--;
1230 argv++;
1231 }
1232
1233 switch ( argc ) {
1234 case 0:
1235 case 1:
1236 case 2:
1237 default:
1238 usage ();
1239 return 1;
1240 case 3:
1241 usage();
1242 break;
1243 }
1244
1245 freq1 = readwave ( A, MAX_LEN, argv[1], &len_A );
1246 DC_cancel ( A, len_A );
1247 freq2 = readwave ( B, MAX_LEN, argv[2], &len_B );
1248 DC_cancel ( B, len_B );
1249
1250 if ( freq1 == 65534 && freq2 != 65534 )
1251 freq1 = freq2;
1252 else if ( freq2 == 65534 && freq1 != 65534 )
1253 freq2 = freq1;
1254 else if ( freq1 == 65534 && freq2 == 65534 )
1255 freq1 = freq2 = 44100;
1256
1257 if ( freq1 != freq2 ) {
1258 fprintf ( stderr, "Different sample frequencies currently not supported\n");
1259 fprintf ( stderr, "A: %ld, B: %ld\n", freq1, freq2 );
1260 return 2;
1261 }
1262
1263 len = len_A < len_B ? len_A : len_B;
1264 fshift = cross_analyze ( A, B, len );
1265 shift = floor ( fshift + 0.5 );
1266
1267 if ( verbose ) {
1268 fprintf ( stderr, "Delay Ch1 is %.4f samples\n", fshift );
1269 fprintf ( stderr, "Delay Ch2 is %.4f samples\n",
1270 cross_analyze ( (stereo_t*)(((sample_t*)A)+1), (stereo_t*)(((sample_t*)B)+1), len ) );
1271 }
1272
1273 if (shift > 0) {
1274 if (verbose)
1275 fprintf ( stderr, "Delaying A by %d samples\n", +shift);
1276 B += shift;
1277 len_B -= shift;
1278 }
1279 if (shift < 0) {
1280 if (verbose)
1281 fprintf ( stderr, "Delaying B by %d samples\n", -shift);
1282 A -= shift;
1283 len_A += shift;
1284 }
1285
1286 len = len_A < len_B ? len_A : len_B;
1287 memset ( &k, 0, sizeof(k) );
1288 analyze_stereo ( A, B, len, &k );
1289 ampl = report (&k);
1290 max_A = maximum ( A, len );
1291 max_B = maximum ( B, len );
1292
1293 if ( ampl <= 0.98855 ) { /* < -0.05 dB */
1294 max = max_A*ampl < max_B ? max_B : max_A*ampl;
1295 ampl_X = (int)(29203 / max);
1296 if ( ampl_X < 2 ) ampl_X = 1;
1297 multiply ( 'A', A, len, ampl*ampl_X );
1298 multiply ( 'B', B, len, ampl_X );
1299 } else if ( ampl >= 1.01158 ) { /* > +0.05 dB */
1300 max = max_A < max_B/ampl ? max_B/ampl : max_A;
1301 ampl_X = (int)(29203 / max);
1302 if ( ampl_X < 2 ) ampl_X = 1;
1303 multiply ( 'A', A, len, ampl_X );
1304 multiply ( 'B', B, len, 1./ampl*ampl_X );
1305 } else {
1306 max = max_A < max_B ? max_B : max_A;
1307 ampl_X = (int)(29203 / max);
1308 if ( ampl_X < 2 ) ampl_X = 1;
1309 multiply ( 'A', A, len, ampl_X );
1310 multiply ( 'B', B, len, ampl_X );
1311 }
1312
1313 set ();
1314 Set_Realtime ();
1315 testing ( A, B, len, freq1 );
1316 reset ();
1317
1318 free (_A);
1319 free (_B);
1320 return 0;
1321 }
1322
1323 /* end of abx.c */
1324