// dts, NPR // dts, PR #include #include #include #ifndef M_PI # define M_PI 3.1415926535897932384626433832795029 #endif #define MAX 8192 #define TYPE double ////////////////////////////// // // BesselI0 -- Regular Modified Cylindrical Bessel Function (Bessel I). // static double Bessel_I_0 ( double x ) { double denominator; double numerator; double z; if (x == 0.) return 1.; z = x * x; numerator = z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* 0.210580722890567e-22 + 0.380715242345326e-19 ) + 0.479440257548300e-16) + 0.435125971262668e-13 ) + 0.300931127112960e-10) + 0.160224679395361e-07 ) + 0.654858370096785e-05) + 0.202591084143397e-02 ) + 0.463076284721000e+00) + 0.754337328948189e+02 ) + 0.830792541809429e+04) + 0.571661130563785e+06 ) + 0.216415572361227e+08) + 0.356644482244025e+09 ) + 0.144048298227235e+10; denominator = z* (z* (z - 0.307646912682801e+04) + 0.347626332405882e+07) - 0.144048298227235e+10; return - numerator / denominator; } static double residual ( double x ) { return sqrt ( 1. - x*x ); } ////////////////////////////// // // KBDWindow -- Kaiser Bessel Derived Window // fills the input window array with size samples of the // KBD window with the given tuning parameter alpha. // static void KBDWindow ( TYPE* window, unsigned int size, TYPE alpha ) { double sumvalue = 0.; unsigned int i; for ( i = 0; i < size/2; i++ ) window [i] = sumvalue += Bessel_I_0 ( M_PI * alpha * residual (4.*i/size - 1.) ); // need to add one more value to the nomalization factor at size/2: sumvalue += Bessel_I_0 ( M_PI * alpha * residual (4.*(size/2)/size-1.) ); // normalize the window and fill in the righthand side of the window: for ( i = 0; i < size/2; i++ ) window [size-1-i] = window [i] = sqrt ( window [i] / sumvalue ); } static void CosWindow ( TYPE* window, unsigned int size ) { double x; unsigned int i; for ( i = 0; i < size/2; i++ ) { x = cos ( (i+0.5) * (M_PI / size) ); window [size/2-1-i] = window [size/2+i] = x; } } static void CosSinWindow ( TYPE* window, unsigned int size ) { double x; unsigned int i; for ( i = 0; i < size/2; i++ ) { x = sin ( (i+0.5) * (M_PI / size) ); x = cos (x * x * M_PI/2); window [size/2-1-i] = window [size/2+i] = x; } } static void SincWindow ( TYPE* window, unsigned int size ) { static long W [] = { 0, -1, -1, -1, -1, -1, -1, -2, -2, -2, -2, -3, -3, -4, -4, -5, -5, -6, -7, -7, -8, -9, -10, -11, -13, -14, -16, -17, -19, -21, -24, -26, -29, -31, -35, -38, -41, -45, -49, -53, -58, -63, -68, -73, -79, -85, -91, -97, -104, -111, -117, -125, -132, -139, -147, -154, -161, -169, -176, -183, -190, -196, -202, -208, -213, -218, -222, -225, -227, -228, -228, -227, -224, -221, -215, -208, -200, -189, -177, -163, -146, -127, -106, -83, -57, -29, 2, 36, 72, 111, 153, 197, 244, 294, 347, 401, 459, 519, 581, 645, 711, 779, 848, 919, 991, 1064, 1137, 1210, 1283, 1356, 1428, 1498, 1567, 1634, 1698, 1759, 1817, 1870, 1919, 1962, 2001, 2032, 2057, 2075, 2085, 2087, 2080, 2063, 2037, 2000, 1952, 1893, 1822, 1739, 1644, 1535, 1414, 1280, 1131, 970, 794, 605, 402, 185, -45, -288, -545, -814,-1095,-1388,-1692,-2006,-2330,-2663,-3004,-3351,-3705,-4063,-4425,-4788, -5153,-5517,-5879,-6237,-6589,-6935,-7271,-7597,-7910,-8209,-8491,-8755,-8998,-9219,-9416,-9585, -9727,-9838,-9916,-9959,-9966,-9935,-9863,-9750,-9592,-9389,-9139,-8840,-8492,-8092,-7640,-7134, -6574,-5959,-5288,-4561,-3776,-2935,-2037,-1082, -70, 998, 2122, 3300, 4533, 5818, 7154, 8540, 9975,11455,12980,14548,16155,17799,19478,21189,22929,24694,26482,28289,30112,31947,33791,35640, 37489,39336,41176,43006,44821,46617,48390,50137,51853,53534,55178,56778,58333,59838,61289,62684, 64019,65290,66494,67629,68692,69679,70590,71420,72169,72835,73415,73908,74313,74630,74856,74992, 75038, 0 }; double x; unsigned int i; unsigned int j; for ( i = 0; i <= size/2; i++ ) { x = i * 512. / size; j = (int) x; x = W [j] * (1 - x + j) + W [j+1] * (x - j); x /= 75038.; window [size-1-i] = window [i] = x; } } #pragma warning ( disable: 4035 ) static void Multiply ( TYPE* z, const TYPE* x, const TYPE* y ) { double X [MAX]; int i; int j; int maxx = 0; int maxy = 0; for ( i = 0; i < MAX; i++ ) { if ( x[i] != 0. ) maxx = i+1; if ( y[i] != 0. ) maxy = i+1; } memset ( X, 0, sizeof X ); for ( i = 0; i < maxx; i++ ) for ( j = 0; j < maxy; j++ ) X [i+j] += (double)x[i] * y[j]; for ( i = 0; i < MAX; i++ ) z [i] = X [i]; } double dB ( double x ) { return 10 * log10 (1.e-99 + x*x); } const char* dir; int flag; float xmin1; float xmax1; float ymin1; float ymax1; float xmin2; float xmax2; float ymin2; float ymax2; float xmin3; float xmax3; float ymin3; float ymax3; void Setup ( FILE* fp, float xmin, float xmax, float ymin, float ymax ) { fprintf ( fp, "@ world xmin %f\n" "@ world xmax %f\n" "@ world ymin %f\n" "@ world ymax %f\n" "@ view xmin 0.08\n" "@ view xmax 0.96\n" "@ view ymin 0.08\n" "@ view ymax 0.96\n", xmin, xmax, ymin, ymax ); } static void Message ( const char* name, TYPE* x, const TYPE sf, int fftsize ) { char filename [128]; TYPE y [MAX/2 + 1]; TYPE S [MAX + 1]; TYPE C [MAX + 1]; double tmpc; double tmps; double tmp; FILE* fp; TYPE maxx = 0.; int maxi = 0; int maxu = 0; int i; int j; int w; if ( fftsize == 0 ) fftsize = MAX; fprintf ( stderr, "Impulse response: %s\n", name ); for ( i = 0; i < MAX; i++ ) { if ( fabs (x[i]) > maxx ) { maxi = i; maxx = fabs (x[i]); } if ( x[i] != 0. ) maxu = i + 1; } mkdir (dir, 0777); sprintf ( filename, "%s/%s time.txt", dir, name ); if ( flag & 1 ) { fp = fopen ( filename, "w" ); Setup ( fp, xmin1, xmax1, ymin1, ymax1 ); for ( i = 0; i < maxu; i++ ) fprintf ( fp, "%8.4f\t%12.9f\n", (i - maxi) * 1000. / sf , x[i] / maxx ); fclose ( fp ); } fprintf ( stderr, "Frequency response: %s\n", name ); for ( j = 0; j <= fftsize/4; j++ ) { tmp = 2. * M_PI / fftsize * j; C[fftsize/2-j] = C[fftsize/2+j] = - ( C[fftsize -j] = C[j] = cos (tmp) ); S[fftsize -j] = S[fftsize/2+j] = - ( S[fftsize/2-j] = S[j] = sin (tmp) ); } for ( i = 0; i <= fftsize/2; i++ ) { tmpc = 0.; tmps = 0.; for ( j = 0; j < maxu; j++ ) { w = i*j & (fftsize - 1); tmpc += x[j] * C[w]; tmps += x[j] * S[w]; } tmp = tmpc*tmpc + tmps*tmps; y [i] = sqrt ( tmp ); } sprintf ( filename, "%s/%s freq.txt", dir, name ); if ( flag & 2 ) { fp = fopen ( filename, "w" ); Setup ( fp, xmin2, xmax2, ymin2, ymax2 ); for ( i = 0; i <= fftsize/2; i++ ) fprintf ( fp, "%8.2f\t%12.9f\n", i * sf / fftsize, y[i] / y[0] ); fclose ( fp ); } sprintf ( filename, "%s/%s frq log.txt", dir, name ); if ( flag & 4 ) { fp = fopen ( filename, "w" ); Setup ( fp, xmin3, xmax3, ymin3, ymax3 ); for ( i = 0; i <= fftsize/2; i++ ) fprintf ( fp, "%8.2f\t%12.9f\n", i * sf / fftsize, dB (y[i] / y[0]) ); fclose ( fp ); } } int #ifdef _WIN32 _cdecl #endif main ( int argc, char** argv ) { TYPE A [MAX]; TYPE B [MAX]; TYPE sf = argc == 1 ? 44100. : atof (argv[1]); if ( sf <= 192. ) sf *= 1000.; if ( sf < 3000. ) sf = 44100.; dir = "Overall impulse response"; flag = 7; memset ( A, 0, sizeof A ); CosWindow ( A, 512 ); Multiply ( A, A, A ); Message ( "Klemm1", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 2 ); Multiply ( A, A, A ); Message ( "Klemm2", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 2.5 ); Multiply ( A, A, A ); Message ( "Klemm2.5", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 3 ); Multiply ( A, A, A ); Message ( "Klemm3", A, sf, 0 ); memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Multiply ( A, A, A ); Message ( "Layer 1, Layer 2, MPC, dts (NPR)", A, sf, 0 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 2048 ); Multiply ( A, A, A ); Message ( "Ogg Vorbis (Long)", A, sf, 0 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 256 ); Multiply ( A, A, A ); Message ( "Ogg Vorbis (Short)", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 5.0 ); Multiply ( A, A, A ); Message ( "AC-3 (Long)", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 256, 5.0 ); Multiply ( A, A, A ); Message ( "AC-3 (Short)", A, sf, 0 ); memset ( A, 0, sizeof A ); CosWindow ( A, 2048 ); Multiply ( A, A, A ); Message ( "AAC (Long, Cos)", A, sf, 0 ); memset ( A, 0, sizeof A ); CosWindow ( A, 256 ); Multiply ( A, A, A ); Message ( "AAC (Short, Cos)", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 2048, 4.0 ); Multiply ( A, A, A ); Message ( "AAC (Long, KBD)", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 256, 6.0 ); Multiply ( A, A, A ); Message ( "AAC (Short, KBD)", A, sf, 0 ); memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Multiply ( A, A, A ); memset ( B, 0, sizeof B ); CosWindow ( B, 1152 ); Multiply ( B, B, B ); Multiply ( A, A, B ); Message ( "MP3 (Long)", A, sf, 0 ); memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Multiply ( A, A, A ); memset ( B, 0, sizeof B ); CosWindow ( B, 384 ); Multiply ( B, B, B ); Multiply ( A, A, B ); Message ( "MP3 (Short)", A, sf, 0 ); memset ( A, 0, sizeof A ); SincWindow ( A, 1024 ); Multiply ( A, A, A ); memset ( B, 0, sizeof B ); CosWindow ( B, 2304 ); Multiply ( B, B, B ); Multiply ( A, A, B ); Message ( "MP3Pro (Long)", A, sf, 0 ); memset ( A, 0, sizeof A ); SincWindow ( A, 1024 ); Multiply ( A, A, A ); memset ( B, 0, sizeof B ); CosWindow ( B, 768 ); Multiply ( B, B, B ); Multiply ( A, A, B ); Message ( "MP3Pro (Short)", A, sf, 0 ); dir = "Frequency resolution (long)"; flag = 6; memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Message ( "Layer 1 + 2", A, sf, 0 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 2048 ); Message ( "Ogg Vorbis", A, sf, 0 ); memset ( A, 0, sizeof A ); CosWindow ( A, 2048 ); Message ( "AAC cos" , A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 2048, 4.0 ); Message ( "AAC KBD", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 2048, 5.0 ); Message ( "AC-3", A, sf, 0 ); dir = "Frequency resolution (short)"; flag = 6; memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Message ( "Layer 1 + 2", A, sf, 0 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 256 ); Message ( "Ogg Vorbis", A, sf, 0 ); memset ( A, 0, sizeof A ); CosWindow ( A, 256 ); Message ( "AAC cos" , A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 256, 6.0 ); Message ( "AAC KBD", A, sf, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 256, 5.0 ); Message ( "AC-3", A, sf, 0 ); dir = "Generic Shape"; flag = 1; xmin2 = -500; xmax2 = 500; ymin2 = 0; ymax2 = 1; memset ( A, 0, sizeof A ); SincWindow ( A, 4096 ); Message ( "Sinc", A, 4096, 0 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 4096 ); Message ( "CosSin", A, 4096, 0 ); memset ( A, 0, sizeof A ); CosWindow ( A, 4096 ); Message ( "Cos" , A, 4096, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 4096, 4.0 ); Message ( "KBD-4", A, 4096, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 4096, 5.0 ); Message ( "KBD-5", A, 4096, 0 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 4096, 6.0 ); Message ( "KBD-6", A, 4096, 0 ); dir = "Generic Frequency Resolution"; flag = 6; xmin2 = 0; xmax2 = 500; ymin2 = 0; ymax2 = 1; xmin3 = 0; xmax3 = 1000; ymin3 = -100; ymax3 = 0; memset ( A, 0, sizeof A ); SincWindow ( A, 512 ); Message ( "Sinc", A, 48000, 512 ); memset ( A, 0, sizeof A ); CosSinWindow ( A, 512 ); Message ( "CosSin", A, 48000, 512 ); memset ( A, 0, sizeof A ); CosWindow ( A, 512 ); Message ( "Cos" , A, 48000, 512 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 4.0 ); Message ( "KBD-4", A, 48000, 512 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 5.0 ); Message ( "KBD-5", A, 48000, 512 ); memset ( A, 0, sizeof A ); KBDWindow ( A, 512, 6.0 ); Message ( "KBD-6", A, 48000, 512 ); return 0; }