Index: /mppenc/trunk/CMakeLists.txt =================================================================== --- /mppenc/trunk/CMakeLists.txt (revision 108) +++ /mppenc/trunk/CMakeLists.txt (revision 109) @@ -4,5 +4,5 @@ if(NOT MSVC) -set(CMAKE_C_FLAGS "-fno-strict-aliasing -Os -fomit-frame-pointer -pipe") +set(CMAKE_C_FLAGS "-Wstrict-aliasing -fstrict-aliasing -ffast-math -O3 -fomit-frame-pointer -pipe") endif(NOT MSVC) Index: /mppenc/trunk/src/cvd.c =================================================================== --- /mppenc/trunk/src/cvd.c (revision 108) +++ /mppenc/trunk/src/cvd.c (revision 109) @@ -219,8 +219,13 @@ logfast ( float x ) { - double y = x * x; - y *= y; - y *= y; - return (((int*)(&y))[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) ); + union { double d; Int32_t n[2]; } tmp; + tmp.d = x * x; + tmp.d *= tmp.d; + tmp.d *= tmp.d; +#if ENDIAN == HAVE_LITTLE_ENDIAN + return (tmp.n[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) ); +#else + return (tmp.n[0] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) ); +#endif } Index: /mppenc/trunk/src/fastmath.c =================================================================== --- /mppenc/trunk/src/fastmath.c (revision 108) +++ /mppenc/trunk/src/fastmath.c (revision 109) @@ -19,4 +19,5 @@ #include "mppenc.h" +#include "fastmath.h" #ifdef FAST_MATH @@ -30,13 +31,5 @@ void Init_FastMath ( void ) { - int i; - float X; - float Y; - double xm; - double x0; - double xp; - double x; - double y; - float* p; + int i; mpc_floatint X, Y; double xm, x0, xp, x, y; float* p; p = (float*) tabatan2; @@ -64,10 +57,10 @@ p = (float*) tabsqrt_ex; for ( i = 0; i < 255; i++ ) { - *(int*)&X = (i << 23); - *(int*)&Y = (i << 23) + (1<<23) - 1; - *p++ = sqrt(X); + X.n = (i << 23); + Y.n = (i << 23) + (1<<23) - 1; + *p++ = sqrt(X.f); } - *(int*)&X = (255 << 23) - 1; - *p++ = sqrt(X); + X.n = (255 << 23) - 1; + *p++ = sqrt(X.f); p = (float*) tabsqrt_m; Index: /mppenc/trunk/src/fastmath.h =================================================================== --- /mppenc/trunk/src/fastmath.h (revision 108) +++ /mppenc/trunk/src/fastmath.h (revision 109) @@ -18,9 +18,25 @@ */ -#if 1 -# define ROUND32(x) ( floattmp = (x) + (int)0x00FD8000L, *(int*)(&floattmp) - (int)0x4B7D8000L ) -#else -# define ROUND32(x) ( (int) floor ((x) + 0.5) ) -#endif +typedef union mpc_floatint +{ + float f; + Int32_t n; +} mpc_floatint; + +static __inline Int32_t mpc_nearbyintf(float fVal) +{ + mpc_floatint tmp; + tmp.f = fVal + 0x00FF8000; + tmp.n = tmp.n - 0x4B7F8000; + return tmp.n; +} +#define mpc_lrintf mpc_nearbyintf +static __inline Int32_t mpc_round32(float fVal) +{ + mpc_floatint tmp; + tmp.f = fVal + 0x00FD8000; + tmp.n = tmp.n - 0x4B7D8000; + return tmp.n; +} #ifdef FAST_MATH @@ -29,21 +45,20 @@ my_atan2 ( float x, float y ) { - float t; - int i; - float ret; - float floattmp; + float t, ret; int i; mpc_floatint mx, my; - if ( (*(int*)&x & 0x7FFFFFFF) < (*(int*)&y & 0x7FFFFFFF) ) { - i = ROUND32 (t = TABSTEP * (x / y)); + mx.f = x; + my.f = y; + if ( (mx.n & 0x7FFFFFFF) < (my.n & 0x7FFFFFFF) ) { + i = mpc_round32 (t = TABSTEP * (mx.f / my.f)); ret = tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (t-i); - if ( *(int*)&y < 0 ) + if ( my.n < 0 ) ret = (float)(ret - M_PI); } - else if ( *(int*)&x < 0) { - i = ROUND32 (t = TABSTEP * (y / x)); + else if ( mx.n < 0 ) { + i = mpc_round32 (t = TABSTEP * (my.f / mx.f)); ret = - M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t); } - else if ( *(int*)&x > 0) { - i = ROUND32 (t = TABSTEP * (y / x)); + else if ( mx.n > 0 ) { + i = mpc_round32 (t = TABSTEP * (my.f / mx.f)); ret = + M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t); } @@ -58,10 +73,6 @@ my_cos ( float x ) { - float t; - int i; - float ret; - float floattmp; - - i = ROUND32 (t = TABSTEP * x); + float t, ret; int i; + i = mpc_round32 (t = TABSTEP * x); ret = tabcos [13*TABSTEP+i][0] + tabcos [13*TABSTEP+i][1] * (t-i); return ret; @@ -72,6 +83,7 @@ my_ifloor ( float x ) { - x = x + (0x0C00000L + 0.500000001); - return *(int*)&x - 1262485505; + mpc_floatint mx; + mx.f = (float) (x + (0x0C00000L + 0.500000001)); + return mx.n - 1262485505; } @@ -80,12 +92,10 @@ my_sqrt ( float x ) { - float ret; - int i; - int ex = *(int*)&x >> 23; // get the exponent - float floattmp; - - *(int*)&x = (*(int*)&x & 0x7FFFFF) | 0x42800000; // delete the exponent - i = ROUND32 (x); // Integer-part of the mantissa (round ????????????) - ret = tabsqrt_m [i-TABSTEP][0] + tabsqrt_m [i-TABSTEP][1] * (x-i); // calculate value + float ret; int i, ex; mpc_floatint mx; + mx.f = x; + ex = mx.n >> 23; // get the exponent + mx.n = (mx.n & 0x7FFFFF) | 0x42800000; // delete the exponent + i = mpc_round32 (mx.f); // Integer-part of the mantissa (round ????????????) + ret = tabsqrt_m [i-TABSTEP][0] + tabsqrt_m [i-TABSTEP][1] * (mx.f-i); // calculate value ret *= tabsqrt_ex [ex]; return ret; Index: /mppenc/trunk/src/quant.c =================================================================== --- /mppenc/trunk/src/quant.c (revision 108) +++ /mppenc/trunk/src/quant.c (revision 109) @@ -19,4 +19,5 @@ #include "mppenc.h" +#include "fastmath.h" /* V A R I A B L E S */ @@ -149,27 +150,21 @@ ISNR_Schaetzer ( const float* input, const float SNRcomp, const int res ) { - int k; - float fac = A [res]; - float invfac = C [res]; - float Signal = 1.e-30f; - float Fehler = 1.e-30f; - float tmp ; - float tmp2; - float tmp3; + int k; float signal, error, sig, err; + const float fac = A [res]; + const float invfac = C [res]; + const float noiseComp = NoiseInjectionCompensation1D [res]; // Summation of the absolute power and the quadratic error - for ( k = 0; k < 36; k++ ) { - tmp2 = input[k] * NoiseInjectionCompensation1D [res]; - // q = ftol(in), correct rounding - tmp = tmp2 * fac + 0xFF8000; - tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac; - tmp = tmp3 - tmp2; - - Fehler += tmp * tmp; - Signal += tmp2 * tmp2; + signal = error = 1.e-30f; + for ( k = 0; k < 36; k++ ) + { + sig = input[k] * noiseComp; + err = mpc_nearbyintf(sig * fac) * invfac - sig; + signal += sig * sig; + error += err * err; } // Utilization of SNRcomp only if SNR > 1 !!! - return Signal > Fehler ? Fehler / (SNRcomp * Signal) : Fehler / Signal; + return signal > error ? error / (signal * SNRcomp) : error / signal; } @@ -178,60 +173,43 @@ ISNR_Schaetzer_Trans ( const float* input, const float SNRcomp, const int res ) { - int k; - float fac = A [res]; - float invfac = C [res]; - float Signal; - float Fehler; - float ret ; - float tmp ; - float tmp2; - float tmp3; - + int k; float signal, error, result, sig, err; + const float fac = A [res]; + const float invfac = C [res]; + const float noiseComp = NoiseInjectionCompensation1D [res]; + // Summation of the absolute power and the quadratic error - k = 0; - Signal = Fehler = 1.e-30f; - for ( ; k < 12; k++ ) { - tmp2 = input[k] * NoiseInjectionCompensation1D [res]; - // q = ftol(in), correct rounding - tmp = tmp2 * fac + 0xFF8000; - tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac; - tmp = tmp3 - tmp2; - - Fehler += tmp * tmp; - Signal += tmp2 * tmp2; - } - tmp = Signal > Fehler ? Fehler / (SNRcomp * Signal) : Fehler / Signal; - ret = tmp; - Signal = Fehler = 1.e-30f; - for ( ; k < 24; k++ ) { - tmp2 = input[k] * NoiseInjectionCompensation1D [res]; - // q = ftol(in), correct rounding - tmp = tmp2 * fac + 0xFF8000; - tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac; - tmp = tmp3 - tmp2; - - Fehler += tmp * tmp; - Signal += tmp2 * tmp2; - } - tmp = Signal > Fehler ? Fehler / (SNRcomp * Signal) : Fehler / Signal; - if ( tmp > ret ) ret = tmp; - //ret += tmp; - Signal = Fehler = 1.e-30f; - for ( ; k < 36; k++ ) { - tmp2 = input[k] * NoiseInjectionCompensation1D [res]; - // q = ftol(in), correct rounding - tmp = tmp2 * fac + 0xFF8000; - tmp3 = (*(int*) & tmp - 0x4B7F8000) * invfac; - tmp = tmp3 - tmp2; - - Fehler += tmp * tmp; - Signal += tmp2 * tmp2; - } - tmp = Signal > Fehler ? Fehler / (SNRcomp * Signal) : Fehler / Signal; - if ( tmp > ret ) ret = tmp; - //ret += tmp; - //ret *= 0.33333333333f; - - return ret; + signal = error = 1.e-30f; + for ( k = 0 ; k < 12; k++ ) + { + sig = input[k] * noiseComp; + err = mpc_nearbyintf(sig * fac) * invfac - sig; + signal += sig * sig; + error += err * err; + } + result = signal > error ? error / (signal * SNRcomp) : error / signal; + + signal = error = 1.e-30f; + for ( ; k < 24; k++ ) + { + sig = input[k] * noiseComp; + err = mpc_nearbyintf(sig * fac) * invfac - sig; + signal += sig * sig; + error += err * err; + } + sig = signal > error ? error / (signal * SNRcomp) : error / signal; + if ( sig > result ) result = sig; + + signal = error = 1.e-30f; + for ( ; k < 36; k++ ) + { + sig = input[k] * noiseComp; + err = mpc_nearbyintf(sig * fac) * invfac - sig; + signal += sig * sig; + error += err * err; + } + sig = signal > error ? error / (signal * SNRcomp) : error / signal; + if ( sig > result ) result = sig; + + return result; } @@ -241,40 +219,37 @@ QuantizeSubband ( unsigned int* qu_output, const float* input, const int res, float* errors ) { - int n; - int offset = D [res]; - float mult = A [res] * NoiseInjectionCompensation1D [res]; - float invmult = C [res]; - float tmp; - int quant; - float signal; - - for ( n = 0; n < 36 - MAX_NS_ORDER; n++, input++, qu_output++ ) { - // q = ftol(in), correct rounding - tmp = *input * mult + 0xFF8000; - quant = (unsigned int)(*(int*) & tmp - 0x4B7F8000 + offset); - + int n, quant; float signal; + const int offset = D [res]; + const float noiseComp = NoiseInjectionCompensation1D [res]; + const float mult = A [res] * noiseComp; + const float invmult = C [res]; + + for ( n = 0; n < 36 - MAX_NS_ORDER; n++) + { + quant = (unsigned int) ( mpc_lrintf( input[n] * mult ) + offset ); // limitation to 0...2D - if ((unsigned int)quant > (unsigned int)2*offset ) { - quant = mini ( quant, 2*offset ); - quant = maxi ( quant, 0 ); + if ( quant > offset * 2 ) + { + quant = mini ( quant, offset * 2 ); + quant = maxi ( quant, 0 ); } - *qu_output = quant; - } - - for ( ; n < 36; n++, input++, qu_output++ ) { - // q = ftol(in), correct rounding - signal = *input * mult; - tmp = signal + 0xFF8000; - quant = (unsigned int)(*(int*) & tmp - 0x4B7F8000 + offset); + qu_output[n] = quant; + } + + for ( ; n < 36; n++ ) + { + signal = input[n] * mult; + quant = (unsigned int) ( mpc_lrintf( signal ) + offset ); // calculate the current error and save it for error refeeding - errors [n + 6] = invmult * (quant - offset) - signal * NoiseInjectionCompensation1D [res]; + errors [n + 6] = invmult * (quant - offset) - signal * noiseComp; // limitation to 0...2D - if ((unsigned int)quant > (unsigned int)2*offset ) { - quant = mini ( quant, 2*offset ); - quant = maxi ( quant, 0 ); + if ( quant > offset * 2 ) + { + quant = mini ( quant, offset * 2); + quant = maxi ( quant, 0 ); } - *qu_output = quant; + qu_output[n] = quant; } } @@ -285,25 +260,21 @@ QuantizeSubbandWithNoiseShaping ( unsigned int* qu_output, const float* input, const int res, float* errors, const float* FIR ) { -#define E(x) *((int*)errors+(x)) - - float signal; - float tmp; - float mult = A [res]; - float invmult = C [res]; - int offset = D [res]; - int n; - int quant; - - E(0) = E(1) = E(2) = E(3) = E(4) = E(5) = 0; // arghh, it produces pops on each frame boundary! - - for ( n = 0; n < 36; n++, input++, qu_output++ ) { - signal = *input * NoiseInjectionCompensation1D [res] - (FIR[5]*errors[n+0] + FIR[4]*errors[n+1] + FIR[3]*errors[n+2] + FIR[2]*errors[n+3] + FIR[1]*errors[n+4] + FIR[0]*errors[n+5]); - - // quant = ftol(signal), correct rounding - tmp = signal * mult + 0xFF8000; - quant = *(int*) & tmp - 0x4B7F8000; + int n, quant; float signal; + const float mult = A [res]; + const float invmult = C [res]; + const int offset = D [res]; + const float noiseComp = NoiseInjectionCompensation1D [res]; + + memset(errors, 0, 6 * sizeof *errors); // arghh, it produces pops on each frame boundary! + + for ( n = 0; n < 36; n++ ) + { + signal = input[n] * noiseComp - + (FIR[5]*errors[n+0] + FIR[4]*errors[n+1] + FIR[3]*errors[n+2] + + FIR[2]*errors[n+3] + FIR[1]*errors[n+4] + FIR[0]*errors[n+5]); + quant = mpc_lrintf(signal * mult); // calculate the current error and save it for error refeeding - errors [n + 6] = invmult * quant - signal * NoiseInjectionCompensation1D [res]; + errors [n + 6] = invmult * quant - signal * noiseComp; // limitation to +/-D @@ -311,5 +282,5 @@ quant = maxf ( quant, -offset ); - *qu_output = (unsigned int)(quant + offset); + qu_output[n] = (unsigned int) ( quant + offset ); } }