mlame_corr.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220
  1. /*
  2. Bug:
  3. - runs only on little endian machines for WAV files
  4. - Not all WAV files are recognized
  5. */
  6. #include <stdio.h>
  7. #include <unistd.h>
  8. #include <math.h>
  9. #include <sys/types.h>
  10. #include <sys/stat.h>
  11. #include <fcntl.h>
  12. #include <memory.h>
  13. typedef signed short stereo [2];
  14. typedef signed short mono;
  15. typedef struct {
  16. unsigned long long n;
  17. long double x;
  18. long double x2;
  19. long double y;
  20. long double y2;
  21. long double xy;
  22. } korr_t;
  23. void analyze_stereo ( const stereo* p, size_t len, korr_t* k )
  24. {
  25. long double _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
  26. double t1;
  27. double t2;
  28. k -> n += len;
  29. for ( ; len--; p++ ) {
  30. _x += (t1 = (*p)[0]); _x2 += t1 * t1;
  31. _y += (t2 = (*p)[1]); _y2 += t2 * t2;
  32. _xy += t1 * t2;
  33. }
  34. k -> x += _x ;
  35. k -> x2 += _x2;
  36. k -> y += _y ;
  37. k -> y2 += _y2;
  38. k -> xy += _xy;
  39. }
  40. void analyze_dstereo ( const stereo* p, size_t len, korr_t* k )
  41. {
  42. static double l0 = 0;
  43. static double l1 = 0;
  44. long double _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
  45. double t1;
  46. double t2;
  47. k -> n += len;
  48. for ( ; len--; p++ ) {
  49. _x += (t1 = (*p)[0] - l0); _x2 += t1 * t1;
  50. _y += (t2 = (*p)[1] - l1); _y2 += t2 * t2;
  51. _xy += t1 * t2;
  52. l0 = (*p)[0];
  53. l1 = (*p)[1];
  54. }
  55. k -> x += _x ;
  56. k -> x2 += _x2;
  57. k -> y += _y ;
  58. k -> y2 += _y2;
  59. k -> xy += _xy;
  60. }
  61. void analyze_mono ( const mono* p, size_t len, korr_t* k )
  62. {
  63. long double _x = 0, _x2 = 0;
  64. double t1;
  65. k -> n += len;
  66. for ( ; len--; p++ ) {
  67. _x += (t1 = (*p)); _x2 += t1 * t1;
  68. }
  69. k -> x += _x ;
  70. k -> x2 += _x2;
  71. k -> y += _x ;
  72. k -> y2 += _x2;
  73. k -> xy += _x2;
  74. }
  75. void analyze_dmono ( const mono* p, size_t len, korr_t* k )
  76. {
  77. static double l0 = 0;
  78. long double _x = 0, _x2 = 0;
  79. double t1;
  80. k -> n += len;
  81. for ( ; len--; p++ ) {
  82. _x += (t1 = (*p) - l0); _x2 += t1 * t1;
  83. l0 = *p;
  84. }
  85. k -> x += _x ;
  86. k -> x2 += _x2;
  87. k -> y += _x ;
  88. k -> y2 += _x2;
  89. k -> xy += _x2;
  90. }
  91. int sgn ( long double x )
  92. {
  93. if ( x > 0 ) return +1;
  94. if ( x < 0 ) return -1;
  95. return 0;
  96. }
  97. int report ( const korr_t* k )
  98. {
  99. long double scale = sqrt ( 1.e5 / (1<<29) ); // Sine Full Scale is +10 dB, 7327 = 100%
  100. long double r;
  101. long double rd;
  102. long double sx;
  103. long double sy;
  104. long double x;
  105. long double y;
  106. long double b;
  107. r = (k->x2*k->n - k->x*k->x) * (k->y2*k->n - k->y*k->y);
  108. r = r > 0.l ? (k->xy*k->n - k->x*k->y) / sqrt (r) : 1.l;
  109. sx = k->n > 1 ? sqrt ( (k->x2 - k->x*k->x/k->n) / (k->n - 1) ) : 0.l;
  110. sy = k->n > 1 ? sqrt ( (k->y2 - k->y*k->y/k->n) / (k->n - 1) ) : 0.l;
  111. x = k->n > 0 ? k->x/k->n : 0.l;
  112. y = k->n > 0 ? k->y/k->n : 0.l;
  113. b = atan2 ( sy, sx * sgn (r) ) * ( 8 / M_PI);
  114. // 6 5 4 3 2
  115. // _______________________________
  116. // |\ | /|
  117. // 7 | \ | / | 1
  118. // | \ | / |
  119. // | \ | / |
  120. // | \ | / |
  121. // 8 |--------------+--------------| 0
  122. // | / | \ |
  123. // | / | \ |
  124. // -7 | / | \ | -1
  125. // | / | \ |
  126. // |/_____________|_____________\|
  127. //
  128. // -6 -5 -4 -3 -2
  129. if ( r > 0.98 ) {
  130. printf ("-mm"); // disguised mono file
  131. return;
  132. }
  133. if ( fabs (b-2) > 0.666 ) {
  134. printf ("-ms"); // low profit for joint stereo
  135. return;
  136. }
  137. if ( r < 0.333 ) {
  138. printf ("-ms"); // low profit for joint stereo
  139. return;
  140. }
  141. }
  142. void readfile ( const char* name, int fd )
  143. {
  144. unsigned short header [22];
  145. stereo s [4096];
  146. mono m [8192];
  147. size_t samples;
  148. korr_t k0;
  149. korr_t k1;
  150. memset ( &k0, 0, sizeof(k0) );
  151. memset ( &k1, 0, sizeof(k1) );
  152. read ( fd, header, sizeof(header) );
  153. switch ( header[11] ) {
  154. case 1:
  155. printf ("-mm\n");
  156. break;
  157. case 2:
  158. while ( ( samples = read (fd, s, sizeof(s)) ) > 0 ) {
  159. analyze_stereo ( s, samples / sizeof (*s), &k0 );
  160. analyze_dstereo ( s, samples / sizeof (*s), &k1 );
  161. }
  162. report (&k0);
  163. report (&k1);
  164. break;
  165. default:
  166. fprintf ( stderr, "%u Channels not supported: %s\n", header[11], name );
  167. break;
  168. }
  169. }
  170. int main ( int argc, char** argv )
  171. {
  172. char* name;
  173. int fd;
  174. if (argc < 2)
  175. readfile ( "<stdin>", 0 );
  176. else
  177. while ( (name = *++argv) != NULL ) {
  178. if ( (fd = open ( name, O_RDONLY )) >= 0 ) {
  179. readfile ( name, fd );
  180. close ( fd );
  181. } else {
  182. fprintf ( stderr, "Can't open: %s\n", name );
  183. }
  184. }
  185. return 0;
  186. }