Index: /libmpc/branches/r2d/libmpcenc/bitstream.c =================================================================== --- /libmpc/branches/r2d/libmpcenc/bitstream.c (revision 220) +++ /libmpc/branches/r2d/libmpcenc/bitstream.c (revision 221) @@ -108,20 +108,4 @@ } -void writeBits (mpc_encoder_t * e, mpc_uint32_t input, unsigned int bits ) -{ - e->outputBits += bits; - - if (e->bitsCount + bits > sizeof(e->bitsBuff) * 8) { - int tmp = (sizeof(e->bitsBuff) * 8 - e->bitsCount); - bits -= tmp; - e->bitsBuff = (e->bitsBuff << tmp) | (input >> bits); - e->bitsCount = sizeof(e->bitsBuff) * 8; - emptyBits(e); - input &= (1 << bits) - 1; - } - e->bitsBuff = (e->bitsBuff << bits) | input; - e->bitsCount += bits; -} - unsigned int encodeSize(mpc_uint64_t size, char * buff, mpc_bool_t addCodeSize) { Index: /libmpc/branches/r2d/libmpcenc/encode_sv7.c =================================================================== --- /libmpc/branches/r2d/libmpcenc/encode_sv7.c (revision 220) +++ /libmpc/branches/r2d/libmpcenc/encode_sv7.c (revision 221) @@ -129,7 +129,13 @@ int n; const Huffman_t * Table, * Tables[2]; + mpc_int32_t * Res_L = e->Res_L; + mpc_int32_t * Res_R = e->Res_R; + mpc_bool_t * DSCF_Flag_L = e->DSCF_Flag_L; + mpc_bool_t * DSCF_Flag_R = e->DSCF_Flag_R; + mpc_int32_t * SCF_Last_L = e->SCF_Last_L; + mpc_int32_t * SCF_Last_R = e->SCF_Last_R; for( n = MaxBand; n >= 0; n--) - if (e->Res_L[n] != 0 || e->Res_R[n] != 0) break; + if (Res_L[n] != 0 || Res_R[n] != 0) break; n++; @@ -147,14 +153,14 @@ if (MaxBand) { - writeBits(e, HuffRes[0][e->Res_L[MaxBand - 1]].Code, HuffRes[0][e->Res_L[MaxBand - 1]].Length); - writeBits(e, HuffRes[0][e->Res_R[MaxBand - 1]].Code, HuffRes[0][e->Res_R[MaxBand - 1]].Length); + writeBits(e, HuffRes[0][Res_L[MaxBand - 1]].Code, HuffRes[0][Res_L[MaxBand - 1]].Length); + writeBits(e, HuffRes[0][Res_R[MaxBand - 1]].Code, HuffRes[0][Res_R[MaxBand - 1]].Length); for ( n = MaxBand - 2; n >= 0; n--) { - int tmp = e->Res_L[n] - e->Res_L[n + 1]; + int tmp = Res_L[n] - Res_L[n + 1]; if (tmp < 0) tmp += 17; - writeBits(e, HuffRes[e->Res_L[n + 1] > 2][tmp].Code, HuffRes[e->Res_L[n + 1] > 2][tmp].Length); - - tmp = e->Res_R[n] - e->Res_R[n + 1]; + writeBits(e, HuffRes[Res_L[n + 1] > 2][tmp].Code, HuffRes[Res_L[n + 1] > 2][tmp].Length); + + tmp = Res_R[n] - Res_R[n + 1]; if (tmp < 0) tmp += 17; - writeBits(e, HuffRes[e->Res_R[n + 1] > 2][tmp].Code, HuffRes[e->Res_R[n + 1] > 2][tmp].Length); + writeBits(e, HuffRes[Res_R[n + 1] > 2][tmp].Code, HuffRes[Res_R[n + 1] > 2][tmp].Length); } @@ -162,8 +168,9 @@ mpc_uint32_t tmp = 0; int cnt = 0, tot = 0; + mpc_bool_t * MS_Flag = e->MS_Flag; for( n = 0; n < MaxBand; n++) { - if ( e->Res_L[n] != 0 || e->Res_R[n] != 0 ) { - tmp = (tmp << 1) | e->MS_Flag[n]; - cnt += e->MS_Flag[n]; + if ( Res_L[n] != 0 || Res_R[n] != 0 ) { + tmp = (tmp << 1) | MS_Flag[n]; + cnt += MS_Flag[n]; tot++; } @@ -179,5 +186,5 @@ if (e->framesInBlock == 0){ for( n = 0; n < 32; n++) - e->DSCF_Flag_L[n] = e->DSCF_Flag_R[n] = 1; // new block -> force key frame + DSCF_Flag_L[n] = DSCF_Flag_R[n] = 1; // new block -> force key frame } @@ -186,9 +193,9 @@ for ( n = 0; n < MaxBand; n++ ) { int tmp = 0, cnt = -1; - if (e->Res_L[n]) { + if (Res_L[n]) { tmp = (e->SCF_Index_L[n][1] == e->SCF_Index_L[n][0]) * 2 + (e->SCF_Index_L[n][2] == e->SCF_Index_L[n][1]); cnt++; } - if (e->Res_R[n]) { + if (Res_R[n]) { tmp = (tmp << 2) | ((e->SCF_Index_R[n][1] == e->SCF_Index_R[n][0]) * 2 + (e->SCF_Index_R[n][2] == e->SCF_Index_R[n][1])); cnt++; @@ -201,11 +208,12 @@ for ( n = 0; n < MaxBand; n++ ) { - if ( e->Res_L[n] ) { + if ( Res_L[n] ) { int m; - if (e->DSCF_Flag_L[n] == 1) { - writeBits(e, e->SCF_Index_L[n][0] + 6, 7); - e->DSCF_Flag_L[n] = 0; + mpc_int32_t * SCFI_L_n = e->SCF_Index_L[n]; + if (DSCF_Flag_L[n] == 1) { + writeBits(e, SCFI_L_n[0] + 6, 7); + DSCF_Flag_L[n] = 0; } else { - unsigned int tmp = (e->SCF_Index_L[n][0] - e->SCF_Last_L[n] + 31) & 127; + unsigned int tmp = (SCFI_L_n[0] - SCF_Last_L[n] + 31) & 127; if (tmp < 64) writeBits(e, HuffDSCF_2[tmp].Code, HuffDSCF_2[tmp].Length); @@ -216,6 +224,6 @@ } for( m = 0; m < 2; m++){ - if (e->SCF_Index_L[n][m+1] != e->SCF_Index_L[n][m]) { - unsigned int tmp = (e->SCF_Index_L[n][m+1] - e->SCF_Index_L[n][m] + 31) & 127; + if (SCFI_L_n[m+1] != SCFI_L_n[m]) { + unsigned int tmp = (SCFI_L_n[m+1] - SCFI_L_n[m] + 31) & 127; if (tmp < 64) writeBits(e, HuffDSCF_1[tmp].Code, HuffDSCF_1[tmp].Length); @@ -226,13 +234,14 @@ } } - e->SCF_Last_L[n] = e->SCF_Index_L[n][2]; - } - if ( e->Res_R[n] ) { + SCF_Last_L[n] = SCFI_L_n[2]; + } + if ( Res_R[n] ) { int m; - if (e->DSCF_Flag_R[n] == 1) { - writeBits(e, e->SCF_Index_R[n][0] + 6, 7); - e->DSCF_Flag_R[n] = 0; + mpc_int32_t * SCFI_R_n = e->SCF_Index_R[n]; + if (DSCF_Flag_R[n] == 1) { + writeBits(e, SCFI_R_n[0] + 6, 7); + DSCF_Flag_R[n] = 0; } else { - unsigned int tmp = (e->SCF_Index_R[n][0] - e->SCF_Last_R[n] + 31) & 127; + unsigned int tmp = (SCFI_R_n[0] - SCF_Last_R[n] + 31) & 127; if (tmp < 64) writeBits(e, HuffDSCF_2[tmp].Code, HuffDSCF_2[tmp].Length); @@ -243,6 +252,6 @@ } for( m = 0; m < 2; m++){ - if (e->SCF_Index_R[n][m+1] != e->SCF_Index_R[n][m]) { - unsigned int tmp = (e->SCF_Index_R[n][m+1] - e->SCF_Index_R[n][m] + 31) & 127; + if (SCFI_R_n[m+1] != SCFI_R_n[m]) { + unsigned int tmp = (SCFI_R_n[m+1] - SCFI_R_n[m] + 31) & 127; if (tmp < 64) writeBits(e, HuffDSCF_1[tmp].Code, HuffDSCF_1[tmp].Length); @@ -253,5 +262,5 @@ } } - e->SCF_Last_R[n] = e->SCF_Index_R[n][2]; + SCF_Last_R[n] = SCFI_R_n[2]; } } @@ -259,5 +268,5 @@ /*********************************** Samples *********************************/ for ( n = 0; n < MaxBand; n++ ) { - int Res = e->Res_L[n]; + int Res = Res_L[n]; const int * q = e->Q[n].L; static const unsigned int thres[] = {0, 0, 3, 7, 9, 1, 3, 4, 8}; @@ -336,5 +345,5 @@ } - Res = e->Res_R[n]; + Res = Res_R[n]; } while (q == e->Q[n].L && (q = e->Q[n].R)); } Index: /libmpc/branches/r2d/libmpcenc/libmpcenc.h =================================================================== --- /libmpc/branches/r2d/libmpcenc/libmpcenc.h (revision 220) +++ /libmpc/branches/r2d/libmpcenc/libmpcenc.h (revision 221) @@ -94,5 +94,20 @@ mpc_uint32_t writeBlock ( mpc_encoder_t *, const char *, const mpc_bool_t, mpc_uint32_t); void writeMagic (mpc_encoder_t * e); - void writeBits (mpc_encoder_t * e, mpc_uint32_t input, unsigned int bits ); + void emptyBits(mpc_encoder_t * e); + static mpc_inline void writeBits (mpc_encoder_t * e, mpc_uint32_t input, unsigned int bits ) + { + e->outputBits += bits; + + if (e->bitsCount + bits > sizeof(e->bitsBuff) * 8) { + int tmp = (sizeof(e->bitsBuff) * 8 - e->bitsCount); + bits -= tmp; + e->bitsBuff = (e->bitsBuff << tmp) | (input >> bits); + e->bitsCount = sizeof(e->bitsBuff) * 8; + emptyBits(e); + input &= (1 << bits) - 1; + } + e->bitsBuff = (e->bitsBuff << bits) | input; + e->bitsCount += bits; + } void writeSeekTable (mpc_encoder_t * e); void writeBitstream_SV8 ( mpc_encoder_t*, int); Index: /libmpc/branches/r2d/libmpcenc/quant.c =================================================================== --- /libmpc/branches/r2d/libmpcenc/quant.c (revision 220) +++ /libmpc/branches/r2d/libmpcenc/quant.c (revision 221) @@ -158,19 +158,35 @@ ISNR_Schaetzer ( const float* input, const float SNRcomp, const int res ) { - int k; - const float fac = A [res]; - const float invfac = C [res]; + const float fac = A [res] * NoiseInjectionCompensation1D [res]; + const float invfac = C [res] / NoiseInjectionCompensation1D [res]; float signal = 1.e-30f; float error = 1.e-30f; - float err, sig; + const float * in_end = input + 36; // Summation of the absolute power and the quadratic error - for ( k = 0; k < 36; k++ ) { - sig = input[k] * NoiseInjectionCompensation1D [res]; - err = mpc_nearbyintf(sig * fac) * invfac - sig; - + do { + float err; + err = mpc_nearbyintf(input[0] * fac) * invfac - input[0]; error += err * err; - signal += sig * sig; - } + signal += input[0] * input[0]; + + err = mpc_nearbyintf(input[1] * fac) * invfac - input[1]; + error += err * err; + signal += input[1] * input[1]; + + err = mpc_nearbyintf(input[2] * fac) * invfac - input[2]; + error += err * err; + signal += input[2] * input[2]; + + err = mpc_nearbyintf(input[3] * fac) * invfac - input[3]; + error += err * err; + signal += input[3] * input[3]; + + input += 4; + + } while (input < in_end); + + error *= NoiseInjectionCompensation1D [res] * NoiseInjectionCompensation1D [res]; + signal *= NoiseInjectionCompensation1D [res] * NoiseInjectionCompensation1D [res]; // Utilization of SNRcomp only if SNR > 1 !!! Index: /libmpc/branches/r2d/libmpcpsy/fft4g.c =================================================================== --- /libmpc/branches/r2d/libmpcpsy/fft4g.c (revision 220) +++ /libmpc/branches/r2d/libmpcpsy/fft4g.c (revision 221) @@ -26,16 +26,433 @@ static void makewt ( const int nw, int* ip, float* w ); static void makect ( const int nc, int* ip, float* c ); -static mpc_inline void bitrv2 ( const int n, int* ip, float* a ); // -static mpc_inline void cftfsub ( const int n, float* a, float* w ); // -static mpc_inline void rftfsub ( const int n, float* a, int nc, float* c ); // -static mpc_inline void cft1st ( const int n, float* a, float* w ); // -static mpc_inline void cftmdl_i386 ( const int n, const int l, float* a, float* w ); // 5648 -// static mpc_inline void cftmdl_3DNow ( const int n, const int l, float* a, float* w ); // 4954 - -#if 0 -# define cftmdl(n,l,a,w) cftmdl_3DNow ( n, l, a, w ) -#else + # define cftmdl(n,l,a,w) cftmdl_i386 ( n, l, a, w ) -#endif + + +/* -------- child routines -------- */ +static mpc_inline void +bitrv2 ( const int n, int* ip, float* a ) +{ + int j, j1, k, k1, l, m, m2; + float xr, xi, yr, yi; + + ip[0] = 0; + l = n; + m = 1; + while ( (m << 3) < l ) { + l >>= 1; + for ( j = 0; j < m; j++ ) { + ip[m + j] = ip[j] + l; + } + m <<= 1; + } + m2 = 2 * m; + if ( (m << 3) == l ) { + for ( k = 0; k < m; k++ ) { + for ( j = 0; j < k; j++ ) { + j1 = 2 * j + ip[k]; + k1 = 2 * k + ip[j]; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += 2 * m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 -= m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += 2 * m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + j1 = 2 * k + m2 + ip[k]; + k1 = j1 + m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + } else { + for ( k = 1; k < m; k++ ) { + for ( j = 0; j < k; j++ ) { + j1 = 2 * j + ip[k]; + k1 = 2 * k + ip[j]; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + j1 += m2; + k1 += m2; + xr = a[j1]; + xi = a[j1 + 1]; + yr = a[k1]; + yi = a[k1 + 1]; + a[j1] = yr; + a[j1 + 1] = yi; + a[k1] = xr; + a[k1 + 1] = xi; + } + } + } + return; +} + + +static mpc_inline void +cft1st ( const int n, float* a, float* w ) +{ + int j, k1; + float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; + float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; + + x0r = a[ 0] + a[ 2]; + x0i = a[ 1] + a[ 3]; + x1r = a[ 0] - a[ 2]; + x1i = a[ 1] - a[ 3]; + x2r = a[ 4] + a[ 6]; + x2i = a[ 5] + a[ 7]; + x3r = a[ 4] - a[ 6]; + x3i = a[ 5] - a[ 7]; + a[ 0] = x0r + x2r; + a[ 1] = x0i + x2i; + a[ 4] = x0r - x2r; + a[ 5] = x0i - x2i; + a[ 2] = x1r - x3i; + a[ 3] = x1i + x3r; + a[ 6] = x1r + x3i; + a[ 7] = x1i - x3r; + wk1r = w[ 2]; + x0r = a[ 8] + a[10]; + x0i = a[ 9] + a[11]; + x1r = a[ 8] - a[10]; + x1i = a[ 9] - a[11]; + x2r = a[12] + a[14]; + x2i = a[13] + a[15]; + x3r = a[12] - a[14]; + x3i = a[13] - a[15]; + a[ 8] = x0r + x2r; + a[ 9] = x0i + x2i; + a[12] = x2i - x0i; + a[13] = x0r - x2r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[10] = wk1r * (x0r - x0i); + a[11] = wk1r * (x0r + x0i); + x0r = x3i + x1r; + x0i = x3r - x1i; + a[14] = wk1r * (x0i - x0r); + a[15] = wk1r * (x0i + x0r); + + k1 = 0; + j = 16; + do { + k1 += 2; + wk2r = w[k1]; + wk2i = w[k1 + 1]; + wk1r = w[2*k1]; + wk1i = w[2*k1 + 1]; + wk3r = wk1r - 2 * wk2i * wk1i; + wk3i = 2 * wk2i * wk1r - wk1i; + x0r = a[j] + a[j + 2]; + x0i = a[j + 1] + a[j + 3]; + x1r = a[j] - a[j + 2]; + x1i = a[j + 1] - a[j + 3]; + x2r = a[j + 4] + a[j + 6]; + x2i = a[j + 5] + a[j + 7]; + x3r = a[j + 4] - a[j + 6]; + x3i = a[j + 5] - a[j + 7]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + x0r -= x2r; + x0i -= x2i; + a[j + 4] = wk2r * x0r - wk2i * x0i; + a[j + 5] = wk2r * x0i + wk2i * x0r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[j + 2] = wk1r * x0r - wk1i * x0i; + a[j + 3] = wk1r * x0i + wk1i * x0r; + x0r = x1r + x3i; + x0i = x1i - x3r; + a[j + 6] = wk3r * x0r - wk3i * x0i; + a[j + 7] = wk3r * x0i + wk3i * x0r; + wk1r = w[2*k1 + 2]; + wk1i = w[2*k1 + 3]; + wk3r = wk1r - 2 * wk2r * wk1i; + wk3i = 2 * wk2r * wk1r - wk1i; + x0r = a[j + 8] + a[j + 10]; + x0i = a[j + 9] + a[j + 11]; + x1r = a[j + 8] - a[j + 10]; + x1i = a[j + 9] - a[j + 11]; + x2r = a[j + 12] + a[j + 14]; + x2i = a[j + 13] + a[j + 15]; + x3r = a[j + 12] - a[j + 14]; + x3i = a[j + 13] - a[j + 15]; + a[j + 8] = x0r + x2r; + a[j + 9] = x0i + x2i; + x0r -= x2r; + x0i -= x2i; + a[j + 12] = -wk2i * x0r - wk2r * x0i; + a[j + 13] = -wk2i * x0i + wk2r * x0r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[j + 10] = wk1r * x0r - wk1i * x0i; + a[j + 11] = wk1r * x0i + wk1i * x0r; + x0r = x1r + x3i; + x0i = x1i - x3r; + a[j + 14] = wk3r * x0r - wk3i * x0i; + a[j + 15] = wk3r * x0i + wk3i * x0r; + } while ( j += 16, j < n ); + return; +} + +static mpc_inline void +cftmdl_i386 ( const int n, const int l, float* a, float* w ) +{ + int j, j1, j2, j3, k, k1, m, m2; + float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; + float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; + + m = l << 2; + + for ( j = 0; j < l; j += 2 ) { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + a[j2] = x0r - x2r; + a[j2 + 1] = x0i - x2i; + a[j1] = x1r - x3i; + a[j1 + 1] = x1i + x3r; + a[j3] = x1r + x3i; + a[j3 + 1] = x1i - x3r; + } + + wk1r = w[2]; + for ( j = m; j < l + m; j += 2 ) { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + a[j2] = x2i - x0i; + a[j2 + 1] = x0r - x2r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[j1] = wk1r * (x0r - x0i); + a[j1 + 1] = wk1r * (x0r + x0i); + x0r = x3i + x1r; + x0i = x3r - x1i; + a[j3] = wk1r * (x0i - x0r); + a[j3 + 1] = wk1r * (x0i + x0r); + } + + k1 = 0; + m2 = 2 * m; + for ( k = m2; k < n; k += m2 ) { + k1 += 2; + wk2r = w[k1]; + wk2i = w[k1 + 1]; + wk1r = w[2*k1]; + wk1i = w[2*k1 + 1]; + wk3r = wk1r - 2 * wk2i * wk1i; + wk3i = 2 * wk2i * wk1r - wk1i; + j = k; + do { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + x0r -= x2r; + x0i -= x2i; + a[j2] = wk2r * x0r - wk2i * x0i; + a[j2 + 1] = wk2r * x0i + wk2i * x0r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[j1] = wk1r * x0r - wk1i * x0i; + a[j1 + 1] = wk1r * x0i + wk1i * x0r; + x0r = x1r + x3i; + x0i = x1i - x3r; + a[j3] = wk3r * x0r - wk3i * x0i; + a[j3 + 1] = wk3r * x0i + wk3i * x0r; + } while ( j += 2, j < l + k ); + + wk1r = w[2*k1 + 2]; + wk1i = w[2*k1 + 3]; + wk3r = wk1r - 2 * wk2r * wk1i; + wk3i = 2 * wk2r * wk1r - wk1i; + j = k + m; + do { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + x0r -= x2r; + x0i -= x2i; + a[j2] = -wk2i * x0r - wk2r * x0i; + a[j2 + 1] = -wk2i * x0i + wk2r * x0r; + x0r = x1r - x3i; + x0i = x1i + x3r; + a[j1] = wk1r * x0r - wk1i * x0i; + a[j1 + 1] = wk1r * x0i + wk1i * x0r; + x0r = x1r + x3i; + x0i = x1i - x3r; + a[j3] = wk3r * x0r - wk3i * x0i; + a[j3 + 1] = wk3r * x0i + wk3i * x0r; + } while ( j += 2, j < l+k+m ); + } + return; +} + +static mpc_inline void +cftfsub ( const int n, float* a, float* w ) +{ + int j, j1, j2, j3, l; + float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; + + l = 2; + if ( n > 8 ) { + cft1st ( n, a, w ); + l = 8; + while ( (l << 2) < n ) { + cftmdl ( n, l, a, w ); + l <<= 2; + } + } + if ( (l << 2) == n ) { + j = 0; + do { + j1 = j + l; + j2 = j1 + l; + j3 = j2 + l; + x0r = a[j] + a[j1]; + x0i = a[j + 1] + a[j1 + 1]; + x1r = a[j] - a[j1]; + x1i = a[j + 1] - a[j1 + 1]; + x2r = a[j2] + a[j3]; + x2i = a[j2 + 1] + a[j3 + 1]; + x3r = a[j2] - a[j3]; + x3i = a[j2 + 1] - a[j3 + 1]; + a[j] = x0r + x2r; + a[j + 1] = x0i + x2i; + a[j2] = x0r - x2r; + a[j2 + 1] = x0i - x2i; + a[j1] = x1r - x3i; + a[j1 + 1] = x1i + x3r; + a[j3] = x1r + x3i; + a[j3 + 1] = x1i - x3r; + } while ( j += 2, j < l ); + } else { + j = 0; + do { + j1 = j + l; + x0r = a[j] - a[j1]; + x0i = a[j + 1] - a[j1 + 1]; + a[j] += a[j1]; + a[j + 1] += a[j1 + 1]; + a[j1] = x0r; + a[j1 + 1] = x0i; + } while ( j += 2, j < l ); + } + return; +} + + +static mpc_inline void +rftfsub ( const int n, float* a, int nc, float* c ) +{ + int j, k, kk, ks, m; + float wkr, wki, xr, xi, yr, yi; + + m = n >> 1; + ks = 2 * nc / m; + kk = ks; + j = 2; + k = n; + do { + k -= 2; + nc -= ks; + wkr = 0.5f - c[nc]; + wki = c[kk]; + xr = a[j] - a[k]; + xi = a[j + 1] + a[k + 1]; + yr = wkr * xr - wki * xi; + yi = wkr * xi + wki * xr; + a[j] -= yr; + a[j + 1] -= yi; + a[k] += yr; + a[k + 1] -= yi; + kk += ks; + } while ( j += 2, j < m ); + return; +} // generates lookup-tables @@ -131,521 +548,3 @@ } - -/* -------- child routines -------- */ -static void -bitrv2 ( const int n, int* ip, float* a ) -{ - int j, j1, k, k1, l, m, m2; - float xr, xi, yr, yi; - - ip[0] = 0; - l = n; - m = 1; - while ( (m << 3) < l ) { - l >>= 1; - for ( j = 0; j < m; j++ ) { - ip[m + j] = ip[j] + l; - } - m <<= 1; - } - m2 = 2 * m; - if ( (m << 3) == l ) { - for ( k = 0; k < m; k++ ) { - for ( j = 0; j < k; j++ ) { - j1 = 2 * j + ip[k]; - k1 = 2 * k + ip[j]; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += m2; - k1 += 2 * m2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += m2; - k1 -= m2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += m2; - k1 += 2 * m2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - j1 = 2 * k + m2 + ip[k]; - k1 = j1 + m2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - } else { - for ( k = 1; k < m; k++ ) { - for ( j = 0; j < k; j++ ) { - j1 = 2 * j + ip[k]; - k1 = 2 * k + ip[j]; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - j1 += m2; - k1 += m2; - xr = a[j1]; - xi = a[j1 + 1]; - yr = a[k1]; - yi = a[k1 + 1]; - a[j1] = yr; - a[j1 + 1] = yi; - a[k1] = xr; - a[k1 + 1] = xi; - } - } - } - return; -} - - -static void -cftfsub ( const int n, float* a, float* w ) -{ - int j, j1, j2, j3, l; - float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - l = 2; - if ( n > 8 ) { - cft1st ( n, a, w ); - l = 8; - while ( (l << 2) < n ) { - cftmdl ( n, l, a, w ); - l <<= 2; - } - } - if ( (l << 2) == n ) { - j = 0; - do { - j1 = j + l; - j2 = j1 + l; - j3 = j2 + l; - x0r = a[j] + a[j1]; - x0i = a[j + 1] + a[j1 + 1]; - x1r = a[j] - a[j1]; - x1i = a[j + 1] - a[j1 + 1]; - x2r = a[j2] + a[j3]; - x2i = a[j2 + 1] + a[j3 + 1]; - x3r = a[j2] - a[j3]; - x3i = a[j2 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - a[j2] = x0r - x2r; - a[j2 + 1] = x0i - x2i; - a[j1] = x1r - x3i; - a[j1 + 1] = x1i + x3r; - a[j3] = x1r + x3i; - a[j3 + 1] = x1i - x3r; - } while ( j += 2, j < l ); - } else { - j = 0; - do { - j1 = j + l; - x0r = a[j] - a[j1]; - x0i = a[j + 1] - a[j1 + 1]; - a[j] += a[j1]; - a[j + 1] += a[j1 + 1]; - a[j1] = x0r; - a[j1 + 1] = x0i; - } while ( j += 2, j < l ); - } - return; -} - - -static void -cft1st ( const int n, float* a, float* w ) -{ - int j, k1; - float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; - float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - x0r = a[ 0] + a[ 2]; - x0i = a[ 1] + a[ 3]; - x1r = a[ 0] - a[ 2]; - x1i = a[ 1] - a[ 3]; - x2r = a[ 4] + a[ 6]; - x2i = a[ 5] + a[ 7]; - x3r = a[ 4] - a[ 6]; - x3i = a[ 5] - a[ 7]; - a[ 0] = x0r + x2r; - a[ 1] = x0i + x2i; - a[ 4] = x0r - x2r; - a[ 5] = x0i - x2i; - a[ 2] = x1r - x3i; - a[ 3] = x1i + x3r; - a[ 6] = x1r + x3i; - a[ 7] = x1i - x3r; - wk1r = w[ 2]; - x0r = a[ 8] + a[10]; - x0i = a[ 9] + a[11]; - x1r = a[ 8] - a[10]; - x1i = a[ 9] - a[11]; - x2r = a[12] + a[14]; - x2i = a[13] + a[15]; - x3r = a[12] - a[14]; - x3i = a[13] - a[15]; - a[ 8] = x0r + x2r; - a[ 9] = x0i + x2i; - a[12] = x2i - x0i; - a[13] = x0r - x2r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[10] = wk1r * (x0r - x0i); - a[11] = wk1r * (x0r + x0i); - x0r = x3i + x1r; - x0i = x3r - x1i; - a[14] = wk1r * (x0i - x0r); - a[15] = wk1r * (x0i + x0r); - - k1 = 0; - j = 16; - do { - k1 += 2; - wk2r = w[k1]; - wk2i = w[k1 + 1]; - wk1r = w[2*k1]; - wk1i = w[2*k1 + 1]; - wk3r = wk1r - 2 * wk2i * wk1i; - wk3i = 2 * wk2i * wk1r - wk1i; - x0r = a[j] + a[j + 2]; - x0i = a[j + 1] + a[j + 3]; - x1r = a[j] - a[j + 2]; - x1i = a[j + 1] - a[j + 3]; - x2r = a[j + 4] + a[j + 6]; - x2i = a[j + 5] + a[j + 7]; - x3r = a[j + 4] - a[j + 6]; - x3i = a[j + 5] - a[j + 7]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - x0r -= x2r; - x0i -= x2i; - a[j + 4] = wk2r * x0r - wk2i * x0i; - a[j + 5] = wk2r * x0i + wk2i * x0r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j + 2] = wk1r * x0r - wk1i * x0i; - a[j + 3] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j + 6] = wk3r * x0r - wk3i * x0i; - a[j + 7] = wk3r * x0i + wk3i * x0r; - wk1r = w[2*k1 + 2]; - wk1i = w[2*k1 + 3]; - wk3r = wk1r - 2 * wk2r * wk1i; - wk3i = 2 * wk2r * wk1r - wk1i; - x0r = a[j + 8] + a[j + 10]; - x0i = a[j + 9] + a[j + 11]; - x1r = a[j + 8] - a[j + 10]; - x1i = a[j + 9] - a[j + 11]; - x2r = a[j + 12] + a[j + 14]; - x2i = a[j + 13] + a[j + 15]; - x3r = a[j + 12] - a[j + 14]; - x3i = a[j + 13] - a[j + 15]; - a[j + 8] = x0r + x2r; - a[j + 9] = x0i + x2i; - x0r -= x2r; - x0i -= x2i; - a[j + 12] = -wk2i * x0r - wk2r * x0i; - a[j + 13] = -wk2i * x0i + wk2r * x0r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j + 10] = wk1r * x0r - wk1i * x0i; - a[j + 11] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j + 14] = wk3r * x0r - wk3i * x0i; - a[j + 15] = wk3r * x0i + wk3i * x0r; - } while ( j += 16, j < n ); - return; -} - -// extern void mpc_cdecl cftmdl_3DNow_1 ( const int n, const int l, float* a, float* w ); -// extern void mpc_cdecl cftmdl_3DNow_2 ( const int n, const int l, float* a, float* w ); - - -static void -cftmdl_i386 ( const int n, const int l, float* a, float* w ) -{ - int j, j1, j2, j3, k, k1, m, m2; - float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; - float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; - - m = l << 2; - - for ( j = 0; j < l; j += 2 ) { - j1 = j + l; - j2 = j1 + l; - j3 = j2 + l; - x0r = a[j] + a[j1]; - x0i = a[j + 1] + a[j1 + 1]; - x1r = a[j] - a[j1]; - x1i = a[j + 1] - a[j1 + 1]; - x2r = a[j2] + a[j3]; - x2i = a[j2 + 1] + a[j3 + 1]; - x3r = a[j2] - a[j3]; - x3i = a[j2 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - a[j2] = x0r - x2r; - a[j2 + 1] = x0i - x2i; - a[j1] = x1r - x3i; - a[j1 + 1] = x1i + x3r; - a[j3] = x1r + x3i; - a[j3 + 1] = x1i - x3r; - } - - wk1r = w[2]; - for ( j = m; j < l + m; j += 2 ) { - j1 = j + l; - j2 = j1 + l; - j3 = j2 + l; - x0r = a[j] + a[j1]; - x0i = a[j + 1] + a[j1 + 1]; - x1r = a[j] - a[j1]; - x1i = a[j + 1] - a[j1 + 1]; - x2r = a[j2] + a[j3]; - x2i = a[j2 + 1] + a[j3 + 1]; - x3r = a[j2] - a[j3]; - x3i = a[j2 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - a[j2] = x2i - x0i; - a[j2 + 1] = x0r - x2r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j1] = wk1r * (x0r - x0i); - a[j1 + 1] = wk1r * (x0r + x0i); - x0r = x3i + x1r; - x0i = x3r - x1i; - a[j3] = wk1r * (x0i - x0r); - a[j3 + 1] = wk1r * (x0i + x0r); - } - - k1 = 0; - m2 = 2 * m; - for ( k = m2; k < n; k += m2 ) { - k1 += 2; - wk2r = w[k1]; - wk2i = w[k1 + 1]; - wk1r = w[2*k1]; - wk1i = w[2*k1 + 1]; - wk3r = wk1r - 2 * wk2i * wk1i; - wk3i = 2 * wk2i * wk1r - wk1i; - j = k; - do { - j1 = j + l; - j2 = j1 + l; - j3 = j2 + l; - x0r = a[j] + a[j1]; - x0i = a[j + 1] + a[j1 + 1]; - x1r = a[j] - a[j1]; - x1i = a[j + 1] - a[j1 + 1]; - x2r = a[j2] + a[j3]; - x2i = a[j2 + 1] + a[j3 + 1]; - x3r = a[j2] - a[j3]; - x3i = a[j2 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - x0r -= x2r; - x0i -= x2i; - a[j2] = wk2r * x0r - wk2i * x0i; - a[j2 + 1] = wk2r * x0i + wk2i * x0r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j1] = wk1r * x0r - wk1i * x0i; - a[j1 + 1] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3r * x0r - wk3i * x0i; - a[j3 + 1] = wk3r * x0i + wk3i * x0r; - } while ( j += 2, j < l + k ); - - wk1r = w[2*k1 + 2]; - wk1i = w[2*k1 + 3]; - wk3r = wk1r - 2 * wk2r * wk1i; - wk3i = 2 * wk2r * wk1r - wk1i; - j = k + m; - do { - j1 = j + l; - j2 = j1 + l; - j3 = j2 + l; - x0r = a[j] + a[j1]; - x0i = a[j + 1] + a[j1 + 1]; - x1r = a[j] - a[j1]; - x1i = a[j + 1] - a[j1 + 1]; - x2r = a[j2] + a[j3]; - x2i = a[j2 + 1] + a[j3 + 1]; - x3r = a[j2] - a[j3]; - x3i = a[j2 + 1] - a[j3 + 1]; - a[j] = x0r + x2r; - a[j + 1] = x0i + x2i; - x0r -= x2r; - x0i -= x2i; - a[j2] = -wk2i * x0r - wk2r * x0i; - a[j2 + 1] = -wk2i * x0i + wk2r * x0r; - x0r = x1r - x3i; - x0i = x1i + x3r; - a[j1] = wk1r * x0r - wk1i * x0i; - a[j1 + 1] = wk1r * x0i + wk1i * x0r; - x0r = x1r + x3i; - x0i = x1i - x3r; - a[j3] = wk3r * x0r - wk3i * x0i; - a[j3 + 1] = wk3r * x0i + wk3i * x0r; - } while ( j += 2, j < l+k+m ); - } - return; -} - - -// static void -// cftmdl_3DNow ( const int n, const int l, float* a, float* w ) -// { -// int j, j1, j2, j3, k, k1, m, m2; -// float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; -// float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; -// -// cftmdl_3DNow_1 (n,l,a,w); -// -// m = l << 2; -// k1 = 0; -// m2 = 2 * m; -// for ( k = m2; k < n; k += m2 ) { -// k1 += 2; -// wk2r = w[k1]; -// wk2i = w[k1 + 1]; -// wk1r = w[2*k1]; -// wk1i = w[2*k1 + 1]; -// wk3r = wk1r - 2 * wk2i * wk1i; -// wk3i = 2 * wk2i * wk1r - wk1i; -// j = k; -// do { -// j1 = j + l; -// j2 = j1 + l; -// j3 = j2 + l; -// x0r = a[j] + a[j1]; -// x0i = a[j + 1] + a[j1 + 1]; -// x1r = a[j] - a[j1]; -// x1i = a[j + 1] - a[j1 + 1]; -// x2r = a[j2] + a[j3]; -// x2i = a[j2 + 1] + a[j3 + 1]; -// x3r = a[j2] - a[j3]; -// x3i = a[j2 + 1] - a[j3 + 1]; -// a[j] = x0r + x2r; -// a[j + 1] = x0i + x2i; -// x0r -= x2r; -// x0i -= x2i; -// a[j2] = wk2r * x0r - wk2i * x0i; -// a[j2 + 1] = wk2r * x0i + wk2i * x0r; -// x0r = x1r - x3i; -// x0i = x1i + x3r; -// a[j1] = wk1r * x0r - wk1i * x0i; -// a[j1 + 1] = wk1r * x0i + wk1i * x0r; -// x0r = x1r + x3i; -// x0i = x1i - x3r; -// a[j3] = wk3r * x0r - wk3i * x0i; -// a[j3 + 1] = wk3r * x0i + wk3i * x0r; -// } while ( j += 2, j < l + k ); -// -// wk1r = w[2*k1 + 2]; -// wk1i = w[2*k1 + 3]; -// wk3r = wk1r - 2 * wk2r * wk1i; -// wk3i = 2 * wk2r * wk1r - wk1i; -// j = k + m; -// do { -// j1 = j + l; -// j2 = j1 + l; -// j3 = j2 + l; -// x0r = a[j] + a[j1]; -// x0i = a[j + 1] + a[j1 + 1]; -// x1r = a[j] - a[j1]; -// x1i = a[j + 1] - a[j1 + 1]; -// x2r = a[j2] + a[j3]; -// x2i = a[j2 + 1] + a[j3 + 1]; -// x3r = a[j2] - a[j3]; -// x3i = a[j2 + 1] - a[j3 + 1]; -// a[j] = x0r + x2r; -// a[j + 1] = x0i + x2i; -// x0r -= x2r; -// x0i -= x2i; -// a[j2] = -wk2i * x0r - wk2r * x0i; -// a[j2 + 1] = -wk2i * x0i + wk2r * x0r; -// x0r = x1r - x3i; -// x0i = x1i + x3r; -// a[j1] = wk1r * x0r - wk1i * x0i; -// a[j1 + 1] = wk1r * x0i + wk1i * x0r; -// x0r = x1r + x3i; -// x0i = x1i - x3r; -// a[j3] = wk3r * x0r - wk3i * x0i; -// a[j3 + 1] = wk3r * x0i + wk3i * x0r; -// } while ( j += 2, j < l+k+m ); -// } -// return; -// } - - -static void -rftfsub ( const int n, float* a, int nc, float* c ) -{ - int j, k, kk, ks, m; - float wkr, wki, xr, xi, yr, yi; - - m = n >> 1; - ks = 2 * nc / m; - kk = ks; - j = 2; - k = n; - do { - k -= 2; - nc -= ks; - wkr = 0.5f - c[nc]; - wki = c[kk]; - xr = a[j] - a[k]; - xi = a[j + 1] + a[k + 1]; - yr = wkr * xr - wki * xi; - yi = wkr * xi + wki * xr; - a[j] -= yr; - a[j + 1] -= yi; - a[k] += yr; - a[k + 1] -= yi; - kk += ks; - } while ( j += 2, j < m ); - return; -} - /* end of fft4g.c */ Index: /libmpc/branches/r2d/libmpcpsy/psy.c =================================================================== --- /libmpc/branches/r2d/libmpcpsy/psy.c (revision 220) +++ /libmpc/branches/r2d/libmpcpsy/psy.c (revision 221) @@ -858,5 +858,5 @@ int l; float new_erg; - float th; + float th, TransDetect = m->TransDetect; const float* ep; @@ -874,6 +874,6 @@ if ( new_erg > old_erg [0][k] ) { // bigger than the old? - if ( new_erg > old_erg [0][k] * m->TransDetect || - new_erg > old_erg [1][k] * m->TransDetect*2 ) // is signal transient? + if ( new_erg > old_erg [0][k] * TransDetect || + new_erg > old_erg [1][k] * TransDetect*2 ) // is signal transient? transient [k] = 1; }