Index: mppenc/branches/r2d/libmpcpsy/ans.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/ans.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/ans.c (revision 59) @@ -0,0 +1,305 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +/* + * Depending on how transient it is, it can be further reduced (up to 0=No ANS). + * Estimate coefficient for feedback at Order=1 over Mask_fu - Mask_fo. + * 3 quantization routines: Order=0, Order=1, Order=2...6 + * Order doesn't specify the power of the noise shaping, but only the flexibility of the form. + * Don't reset utilization of the "remains" at the frame borders; + * "remains"-utilization as scalefactor-independent values, + * so that a utilization beyond Subframe/Frame Borders is even possible. + */ + +#include "mppenc.h" + + +static float InvFourier [MAX_NS_ORDER + 1] [16]; +static float Cos_Tab [16] [MAX_NS_ORDER + 1]; +static float Sin_Tab [16] [MAX_NS_ORDER + 1]; +unsigned int NS_Order; // Maximum order for ANS +unsigned int NS_Order_L [32]; +unsigned int NS_Order_R [32]; // frame-wise order of the Noiseshaping (0: off, 1...5: on) +float FIR_L [32] [MAX_NS_ORDER]; +float FIR_R [32] [MAX_NS_ORDER]; // contains FIR-Filter for NoiseShaping +float ANSspec_L [MAX_ANS_LINES]; +float ANSspec_R [MAX_ANS_LINES]; // L/R-masking thresholds for ANS +float ANSspec_M [MAX_ANS_LINES]; +float ANSspec_S [MAX_ANS_LINES]; // M/S-masking thresholds for ANS + + +void +Init_ANS ( void ) +{ + int n; + int k; + + // calculate Fourier tables + for ( k = 0; k <= MAX_NS_ORDER; k++ ) { + for ( n = 0; n < 16; n++ ) { + InvFourier [k] [n] = (float) cos ( +2*M_PI/64 * (2*n) * k ) / 16.; + Cos_Tab [n] [k] = (float) cos ( -2*M_PI/64 * (2*n+1) * (k+1) ); + Sin_Tab [n] [k] = (float) sin ( -2*M_PI/64 * (2*n+1) * (k+1) ); + } + } +} + + +// calculates optimal reflection coefficients and time response of a prediction filter in LPC analysis +static __inline void +durbin_akf_to_kh1( float* k, // out: reflection coefficients + float* h, // out: time response + const float* akf ) // in : autocorrelation function (0..1 used) +{ + h[0] = k[0] = akf [1] / akf [0]; +} + +static __inline void +durbin_akf_to_kh2( float* k, // out: reflection coefficients + float* h, // out: time response + const float* akf ) // in : autocorrelation function (0..2 used) +{ + float tk,e; + + tk = akf [1] / akf[0]; + e = akf[0] * (1. - tk*tk); + h[0] = k[0] = tk; + h[0] *= 1. - (h[1] = k[1] = tk = (akf[2] - h[0] * akf[1]) / e); +} + +static __inline void +durbin_akf_to_kh3( float* k, // out: reflection coefficients + float* h, // out: time response + const float* akf ) // in : autocorrelation function (0..3 used) +{ + float a,b,tk,e; + + tk = akf[1] / akf[0]; + e = akf[0] * (1. - tk*tk); + h[0] = k[0] = tk; + + tk = (akf[2] - h[0] * akf[1]) / e; + e *= 1. - tk*tk; + h[0] *= 1. - (h[1] = k[1] = tk); + h[2] = k[2] = tk = (akf[3] - h[0] * akf[2] - h[1] * akf[1]) / e; + + h[0] = (a=h[0]) - (b=h[1])*tk; + h[1] = b - a*tk; +} + + +static __inline void +durbin_akf_to_kh ( float* k, // out: reflection coefficients + float* h, // out: time response + float* akf, // in : autocorrelation function (0..n used) + const int n ) // in : number of parameters to calculate +{ + int i,j; + float s,a,b,tk,e; + float* p; + float* q; + + e = akf [0]; + for ( i = 0; i < n; i++ ) { + s = 0.f; + p = h; + q = akf+i; + j = i; + while ( j-- ) + s += *p++ * *q--; + + tk = (akf[i+1] - s) / e; + e *= 1. - tk*tk; + h[i] = k[i] = tk; + p = h; + q = h + i - 1; + + for ( ; p < q; p++, q-- ) { + a = *p; + b = *q; + *p = a - b*tk; + *q = b - a*tk; + } + if ( p == q ) + *p *= 1. - tk; + } +} + +static const unsigned char maxANSOrder [32] = { + 6, 5, 4, 3, 2, 2, 2, 2, + 2, 2, 2, 2, 1, 1, 1, 1, + 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, +}; + +static void +FindOptimalANS ( const int MaxBand, + const unsigned char* ms, + const float* spec0, + const float* spec1, + unsigned int* NS, + float* snr_comp, + float fir [] [MAX_NS_ORDER], + const float* smr0, + const float* smr1, + const int scf [] [3], + const int Transient [32] ) +{ + int Band; + int n; + int k; + int order; + float akf [MAX_NS_ORDER + 1]; + float h [MAX_NS_ORDER]; + float reflex [MAX_NS_ORDER]; + float spec [16]; + float invspec [16]; + float norm; + float ns_loss; + float min_spec; + float min_diff; + float re; + float im; + float ns_energy; + float gain; + float NS_Gain; + float actSMR; + int max; + const float* tmp; + + ENTER(235); + for ( Band = 0; Band <= MaxBand && maxANSOrder[Band]; Band++ ) { + + if ( scf[Band][0] != scf[Band][1] || scf[Band][1] != scf[Band][2] ) + continue; + + if ( Transient[Band] ) + continue; + + max = maxANSOrder [Band]; + + if ( ms[Band] ) { // setting pointer and SMR in relation to the M/S-flag + tmp = &spec1 [Band<<4]; // pointer to MS-data + actSMR = smr1 [Band]; // selecting SMR + } + else { + tmp = &spec0 [Band<<4]; // pointer to LR-data + actSMR = smr0 [Band]; // selecting SMR + } + + if ( actSMR >= 1. ) { + NS_Gain = 1.f; // reset gain + norm = 1.e-30f; + + // Selection of the masking threshold of the current subband, also considering frequency inversion in every 2nd subband + if ( Band & 1 ) + for ( n = 0, tmp += 15; n < 16; n++ ) + norm += spec[n] = *tmp--; + else + for ( n = 0; n < 16; n++ ) + norm += spec[n] = *tmp++; + + // Preprocessing: normalization of the the power of spec[] to 1, and search for minimum of masking threshold + norm = 16.f / norm; + min_spec = 1.e+12f; + for ( n = 0; n < 16; n++ ) { + invspec[n] = 1.f / (spec[n] *= norm); + if ( spec[n] < min_spec ) // normalize spec[] + min_spec = spec[n]; + } + + // Calculation of the auto-correlation function + tmp = InvFourier [0]; + for ( k = 0; k <= max; k++, tmp += 16 ) { + akf[k] = tmp[ 0]*invspec[ 0] + tmp[ 1]*invspec[ 1] + tmp[ 2]*invspec[ 2] + tmp[ 3]*invspec[ 3] + + tmp[ 4]*invspec[ 4] + tmp[ 5]*invspec[ 5] + tmp[ 6]*invspec[ 6] + tmp[ 7]*invspec[ 7] + + tmp[ 8]*invspec[ 8] + tmp[ 9]*invspec[ 9] + tmp[10]*invspec[10] + tmp[11]*invspec[11] + + tmp[12]*invspec[12] + tmp[13]*invspec[13] + tmp[14]*invspec[14] + tmp[15]*invspec[15]; + } + + // Searching for the noise-shaper with maximum gain + for ( order = 1; order <= max; order++ ) { + switch ( order ) { // calculating best FIR-Filter for the return + case 1: durbin_akf_to_kh1 (reflex, h, akf); break; + case 2: durbin_akf_to_kh2 (reflex, h, akf); break; + case 3: durbin_akf_to_kh3 (reflex, h, akf); break; + default: durbin_akf_to_kh (reflex, h, akf, order); break; + } + + ns_loss = 1.e-30f; // estimating the gain + min_diff = 1.e+12f; + for ( n = 0; n < 16; n++ ) { + re = 1.f; // calculating the obtained noise shaping + im = 0.f; + for ( k = 0; k < order; k++ ) { + re -= h[k] * Cos_Tab[n][k]; + im += h[k] * Sin_Tab[n][k]; + } + + ns_energy = re*re + im*im; // calculated spectral shaped noise + ns_loss += ns_energy; // noise energy increases with shaping + + if ( spec[n] < min_diff * ns_energy ) // Searching for minimum distance between the shaped noise and the masking threshold + min_diff = spec[n] / ns_energy; + } + + // Updating the Filter if new gain is bigger than old gain and if the extra noise power through shaping is smaller than the SMR of this band + gain = 16. * min_diff / (min_spec * ns_loss); + if ( gain > NS_Gain && ns_loss < actSMR ) { + NS [Band] = order; + NS_Gain = gain; + memcpy ( fir [Band], h, order * sizeof(*h) ); + } + } + + if ( NS_Gain > 1.f ) { // Activation of ANS if there is gain + snr_comp[Band] *= NS_Gain; + } + } + } + + LEAVE(235); + return; +} + + +// perform ANS-analysis (calculation of FIR-filter and gain) +void +NS_Analyse ( const int MaxBand, + const unsigned char* MSflag, + const SMRTyp smr, + const int* Transient ) +{ + ENTER(10); + + // for L or M, respectively + memset ( FIR_L, 0, sizeof FIR_L ); // reset FIR + memset ( NS_Order_L, 0, sizeof NS_Order_L ); // reset Flags + FindOptimalANS ( MaxBand, MSflag, ANSspec_L, ANSspec_M, NS_Order_L, SNR_comp_L, FIR_L, smr.L, smr.M, SCF_Index_L, Transient ); + + // for R or S, respectively + memset ( FIR_R, 0, sizeof FIR_R ); // reset FIR + memset ( NS_Order_R, 0, sizeof NS_Order_R ); // reset Flags + FindOptimalANS ( MaxBand, MSflag, ANSspec_R, ANSspec_S, NS_Order_R, SNR_comp_R, FIR_R, smr.R, smr.S, SCF_Index_R, Transient ); + + LEAVE(10); + return; +} + +/* end of ans.c */ Index: mppenc/branches/r2d/libmpcpsy/cvd.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/cvd.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/cvd.c (revision 59) @@ -0,0 +1,267 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include "mppenc.h" + +/* C O N S T A N T S */ +// from MatLab-Simulation (Fourier-transforms of the Cos-Rolloff) +#if 0 +static const float Puls [11] = { + -0.02724753942504f, -0.10670808991329f, -0.06198987803623f, 0.18006206051664f, + 0.49549552704050f, 0.64201253447071f, 0.49549552704050f, 0.18006206051664f, + -0.06198987803623f, -0.10670808991329f, -0.02724753942504f +}; +#endif + +static const float Puls [ 9] = { + -0.10670808991329f, -0.06198987803623f, 0.18006206051664f, 0.49549552704050f, + 0.64201253447071f, 0.49549552704050f, 0.18006206051664f, -0.06198987803623f, + -0.10670808991329f +}; + +/* +// Generating the Cos-Rolloff of the Cepstral-analysis, Cos-Rolloff from 5512,5 Hz to 11025 Hz +// for ( k = 0; k <= 1024; k++ ) { +// if (k < 256) CosWin [k-256] = 1; +// else if (k < 512) CosWin [k-256] = 0.5 + 0.5*cos (M_PI*(k-256)/256); +// else CosWin [k-256] = 0; +// } +*/ +static const float CosWin [256] = { + 1.0000000000000000f, 0.9999623298645020f, 0.9998494386672974f, 0.9996612071990967f, 0.9993977546691895f, 0.9990590810775757f, 0.9986452460289002f, 0.9981563091278076f, 0.9975923895835877f, 0.9969534873962402f, 0.9962397813796997f, 0.9954513311386108f, 0.9945882558822632f, 0.9936507344245911f, 0.9926388263702393f, 0.9915527701377869f, 0.9903926253318787f, 0.9891586899757385f, 0.9878510832786560f, 0.9864699840545654f, 0.9850156307220459f, 0.9834882616996765f, 0.9818880558013916f, 0.9802152514457703f, 0.9784701466560364f, 0.9766530394554138f, 0.9747641086578369f, 0.9728036522865295f, 0.9707720279693604f, 0.9686695337295532f, 0.9664964079856873f, 0.9642530679702759f, 0.9619397521018982f, 0.9595569372177124f, 0.9571048617362976f, 0.9545840024948120f, 0.9519946575164795f, 0.9493372440338135f, 0.9466121792793274f, 0.9438198208808899f, 0.9409606456756592f, 0.9380350708961487f, 0.9350435137748718f, 0.9319864511489868f, 0.9288643002510071f, 0.9256775975227356f, 0.9224267601966858f, 0.9191123247146606f, 0.9157348275184631f, 0.9122946262359619f, 0.9087924361228943f, 0.9052286148071289f, 0.9016037583351135f, 0.8979184627532959f, 0.8941732048988342f, 0.8903686404228210f, 0.8865052461624146f, 0.8825836181640625f, 0.8786044120788574f, 0.8745682239532471f, 0.8704755902290344f, 0.8663271069526672f, 0.8621235489845276f, 0.8578653931617737f, + 0.8535534143447876f, 0.8491881489753723f, 0.8447702527046204f, 0.8403005003929138f, 0.8357794880867004f, 0.8312078714370728f, 0.8265864253044128f, 0.8219157457351685f, 0.8171966671943665f, 0.8124297261238098f, 0.8076158165931702f, 0.8027555346488953f, 0.7978496551513672f, 0.7928989529609680f, 0.7879040837287903f, 0.7828658819198608f, 0.7777851223945618f, 0.7726625204086304f, 0.7674987912178040f, 0.7622948288917542f, 0.7570513486862183f, 0.7517691850662231f, 0.7464491128921509f, 0.7410919070243835f, 0.7356983423233032f, 0.7302693724632263f, 0.7248056530952454f, 0.7193081378936768f, 0.7137775421142578f, 0.7082147598266602f, 0.7026206851005554f, 0.6969960331916809f, 0.6913416981697083f, 0.6856585741043091f, 0.6799474954605103f, 0.6742093563079834f, 0.6684449315071106f, 0.6626551747322083f, 0.6568408608436585f, 0.6510030031204224f, 0.6451423168182373f, 0.6392598152160645f, 0.6333563923835754f, 0.6274328231811523f, 0.6214900612831116f, 0.6155290603637695f, 0.6095505952835083f, 0.6035556793212891f, 0.5975451469421387f, 0.5915199518203735f, 0.5854809284210205f, 0.5794290900230408f, 0.5733652114868164f, 0.5672903656959534f, 0.5612053275108337f, 0.5551111102104187f, 0.5490085482597351f, 0.5428986549377441f, 0.5367822647094727f, 0.5306603908538818f, 0.5245338082313538f, 0.5184035897254944f, 0.5122706294059753f, 0.5061357617378235f, + 0.5000000000000000f, 0.4938642382621765f, 0.4877294003963471f, 0.4815963804721832f, 0.4754661619663239f, 0.4693396389484406f, 0.4632177054882050f, 0.4571013450622559f, 0.4509914219379425f, 0.4448888897895813f, 0.4387946724891663f, 0.4327096343040466f, 0.4266347587108612f, 0.4205709397792816f, 0.4145190417766571f, 0.4084800481796265f, 0.4024548530578613f, 0.3964443206787109f, 0.3904493749141693f, 0.3844709396362305f, 0.3785099089145660f, 0.3725671768188477f, 0.3666436076164246f, 0.3607401549816132f, 0.3548576533794403f, 0.3489970266819000f, 0.3431591391563416f, 0.3373448550701141f, 0.3315550684928894f, 0.3257906734943390f, 0.3200524747371674f, 0.3143413960933685f, 0.3086582720279694f, 0.3030039668083191f, 0.2973793447017670f, 0.2917852103710175f, 0.2862224578857422f, 0.2806918919086456f, 0.2751943469047546f, 0.2697306573390961f, 0.2643016278743744f, 0.2589081227779388f, 0.2535509169101715f, 0.2482308149337769f, 0.2429486215114594f, 0.2377051562070847f, 0.2325011938810349f, 0.2273375093936920f, 0.2222148776054382f, 0.2171340882778168f, 0.2120959013700485f, 0.2071010768413544f, 0.2021503448486328f, 0.1972444802522659f, 0.1923841983079910f, 0.1875702589750290f, 0.1828033626079559f, 0.1780842244625092f, 0.1734135746955872f, 0.1687921136617661f, 0.1642205268144608f, 0.1596994996070862f, 0.1552297323942184f, 0.1508118808269501f, + 0.1464466154575348f, 0.1421345919370651f, 0.1378764659166336f, 0.1336728632450104f, 0.1295244395732880f, 0.1254318058490753f, 0.1213955804705620f, 0.1174163669347763f, 0.1134947761893272f, 0.1096313893795013f, 0.1058267876505852f, 0.1020815446972847f, 0.0983962342143059f, 0.0947714000940323f, 0.0912075936794281f, 0.0877053514122963f, 0.0842651948332787f, 0.0808876454830170f, 0.0775732174515724f, 0.0743224024772644f, 0.0711356922984123f, 0.0680135712027550f, 0.0649565011262894f, 0.0619649514555931f, 0.0590393692255020f, 0.0561801902949810f, 0.0533878505229950f, 0.0506627671420574f, 0.0480053536593914f, 0.0454160086810589f, 0.0428951233625412f, 0.0404430739581585f, 0.0380602329969406f, 0.0357469581067562f, 0.0335035994648933f, 0.0313304923474789f, 0.0292279683053494f, 0.0271963365375996f, 0.0252359099686146f, 0.0233469791710377f, 0.0215298328548670f, 0.0197847411036491f, 0.0181119665503502f, 0.0165117643773556f, 0.0149843730032444f, 0.0135300243273377f, 0.0121489353477955f, 0.0108413146808743f, 0.0096073597669601f, 0.0084472559392452f, 0.0073611787520349f, 0.0063492907211185f, 0.0054117450490594f, 0.0045486823655665f, 0.0037602325901389f, 0.0030465149320662f, 0.0024076367262751f, 0.0018436938989908f, 0.0013547716662288f, 0.0009409435442649f, 0.0006022718735039f, 0.0003388077020645f, 0.0001505906548118f, 0.0000376490788767f, +}; + + +/* F U N C T I O N S */ +// sets all the harmonics +static void +SetVoiceLines ( int* VoiceLine, const float base, int val ) +{ + int n; + int max = (int) (MAX_CVD_LINE * base / 1024.f); // harmonics up to Index MAX_CVD_LINE (spectral lines outside of that don't make sense) + int line; + float frq = 1024.f / base; // frq = 1024./i is the Index of the basic harmonic + + // go through all harmonics + for ( n = 1; n <= max; n++ ) { + line = (int) (n * frq); + VoiceLine [line] = VoiceLine [line+1] = val; + } +} + + +// Analyze the Cepstrum, search for the basic harmonic +static void +CEP_Analyse2048 ( float* res1, + float* res2, + float* qual1, + float* qual2, + float* cep ) +{ + int n; + int line; + float cc [MAX_ANALYZED_IDX + 3]; // cross correlation + float ref; + float line_sum; + float sum; + float kkf; + float norm; + const float* x; + + // cross-correlation with pulse shape + // Calculate idx = MIN_ANALYZED_IDX-2 to MAX_ANALYZED_IDX+2, + // because they are read during search for maximum + // 50 -> 882 Hz, 700 -> 63 Hz base frequency + + *res1 = *res2 = 0. ; + memset ( cc, 0, sizeof cc ); + + for ( n = MIN_ANALYZED_IDX - 2; n <= MAX_ANALYZED_IDX + 2; n++ ) { + x = cep + n; + if ( x[0] > 0 ) { + norm = x[-4] * x[-4] + + x[-3] * x[-3] + + x[-2] * x[-2] + + x[-1] * x[-1] + + x[ 0] * x[ 0] + + x[ 1] * x[ 1] + + x[ 2] * x[ 2] + + x[ 3] * x[ 3] + + x[ 4] * x[ 4]; + kkf = x[-4] * Puls [0] + + x[-3] * Puls [1] + + x[-2] * Puls [2] + + x[-1] * Puls [3] + + x[ 0] * Puls [4] + + x[ 1] * Puls [5] + + x[ 2] * Puls [6] + + x[ 3] * Puls [7] + + x[ 4] * Puls [8]; + cc [n] = kkf * kkf / norm; // calculate the square of ncc to avoid sqrt() + } + } + + // search for the (relative) maximum + ref = 0.f; + line = MED_ANALYZED_IDX; + for ( n = MAX_ANALYZED_IDX; n >= MED_ANALYZED_IDX; n-- ) { + if ( + cc[n] * cep[n] * cep[n] > ref && + cc[n] > 0.40f && // e33 (02) 0.85 + cep[n] > 0.00f && // e33 (02) + cc[n ] >= cc[n+1] && + cc[n ] >= cc[n-1] && + cc[n+1] >= cc[n+2] && + cc[n-1] >= cc[n-2] + ) + { + ref = cc[n] * cep[n] * cep[n]; + line = n; + } + } + + // Calculating the center of the maximum (Interpolation) + x = cep + line; + sum = x[-3] + x[-2] + x[-1] + x[0] + x[1] + x[2] + x[3] + 1.e-30f; + line_sum = (x[1]-x[-1]) + 2 * (x[2]-x[-2]) + 3 * (x[3]-x[-3]) + sum * line + 1.e-30f; + + /* e33 (04) */ + ref = cc[line ] * cep[line ] * cep[line ] + + cc[line-1] * cep[line-1] * cep[line-1] + + cc[line+1] * cep[line+1] * cep[line+1]; + + //{ + // static unsigned int x = 0; + // + // printf ("%7.3f s ", (x/2)*1152./44100 ); + // x++; + //} + + //printf ("ref=%5.3f *res1=%7.3f f=%8.3f ", ref, line_sum / sum, 44100. / (line_sum / sum) ); + + *qual1 = ref; + if ( ref > 0.015f ) + *res1 = line_sum / sum; + + if ( CVD_used < 2 ) + return; + + // search for the (relative) maximum + ref = 0.f; + line = MIN_ANALYZED_IDX; + + for ( n = MED_ANALYZED_IDX + 1; n >= MIN_ANALYZED_IDX - 1; n-- ) { + cc [2*n ] += 0.5 * cc [n]; + cc [2*n+1] += 0.5 * (cc [n] + cc[n+1]); + cep [2*n ] += 0.5 * cep [n]; + cep [2*n+1] += 0.5 * (cep [n] + cep[n+1]); + } + + for ( n = 2*MED_ANALYZED_IDX; n >= 2*MIN_ANALYZED_IDX; n-- ) { + if ( + cc[n] * cep[n] * cep[n] > ref && + cc[n] > 0.85f && /* e33 (02) */ + cep[n] > 0.00f && /* e33 (02) */ + cc[n ] >= cc[n+1] && + cc[n ] >= cc[n-1] && + cc[n+1] >= cc[n+2] && + cc[n-1] >= cc[n-2] + ) + { + ref = cc[n] * cep[n] * cep[n]; + line = n; + } + } + + // Calculating the center of the maximum (Interpolation) + x = cep + line; + sum = x[-3] + x[-2] + x[-1] + x[0] + x[1] + x[2] + x[3] + 1.e-30f; + line_sum = (x[1]-x[-1]) + 2 * (x[2]-x[-2]) + 3 * (x[3]-x[-3]) + sum * line + 1.e-30f; + + /* e33 (04) */ + ref = cc[line ] * cep[line ] * cep[line ] + + cc[line-1] * cep[line-1] * cep[line-1] + + cc[line+1] * cep[line+1] * cep[line+1]; + + //printf ("ref=%5.3f *res2=%8.3f f=%8.3f\n", ref, 0.5 * line_sum / sum, 44100. / (0.5 * line_sum / sum) ); + + *qual2 = ref; + if ( ref >= 0.1f ) + *res2 = 0.5 * line_sum / sum; + + return; +} + +#ifndef CVD_FASTLOG +# define logfast(x) ((float) log (x)) +#else + +static __inline float /* This is a rough estimation with an accuracy of |x|<0.0037 */ +logfast ( float x ) +{ + double y = x * x; + y *= y; + y *= y; + return (((int*)(&y))[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) ); +} + +#endif + +// ClearVoiceDetection for spectrum *spec +// input : Spectrum *spec +// output: Array *vocal contains information if the FFT-Line is a harmonic component +int +CVD2048 ( const float* spec, int* vocal ) +{ + static float cep [4096]; // cep[4096] -- array, which is also used for the 2048 FFT + const float* win = CosWin; // pointer to cos-roll-off + float res1; + float res2; + float qual1; + float qual2; + int n; + + ENTER(20); + // Calculating logarithmated, windowed spectrum cep[] + // cep[512...1024] = 0 -- cep[1025...2047] doesn't matter, because the first have to be filled by fft + for ( n = 0; n < 256; n++ ) + cep[n] = logfast (*spec++); + for ( n = 256; n < 512; n++ ) + cep[n] = logfast (*spec++) * *win++; + + memset ( cep+512, 0, 513*sizeof(*cep) ); + + // Calculating cepstrum of cep[] (the function Cepstrum() outputs the cepstrum in-place) + Cepstrum2048 ( cep, MAX_ANALYZED_IDX ); + + // search the harmonic + CEP_Analyse2048 ( &res1, &res2, &qual1, &qual2, cep ); +//#include "cvd.h" + if ( res1 > 0.f || res2 > 0.f ) { + if ( res1 > 0. ) SetVoiceLines ( vocal, res1, 100 ); + if ( res2 > 0. ) SetVoiceLines ( vocal, res2, 20 ); + LEAVE(20); + return 1; + } + LEAVE(20); + return 0; +} Index: mppenc/branches/r2d/libmpcpsy/fastmath.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/fastmath.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/fastmath.c (revision 59) @@ -0,0 +1,85 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include "mppenc.h" + +#ifdef FAST_MATH + +const float tabatan2 [ 2*TABSTEP+1] [2]; +const float tabcos [26*TABSTEP+1] [2]; +const float tabsqrt_ex [256]; +const float tabsqrt_m [ TABSTEP+1] [2]; + + +void Init_FastMath ( void ) +{ + int i; + float X; + float Y; + double xm; + double x0; + double xp; + double x; + double y; + float* p; + + p = (float*) tabatan2; + for ( i = -TABSTEP; i <= TABSTEP; i++ ) { + xm = atan ((i-0.5)/TABSTEP); + x0 = atan ((i+0.0)/TABSTEP); + xp = atan ((i+0.5)/TABSTEP); + x = x0/2 + (xm + xp)/4; + y = xp - xm; + *p++ = x; + *p++ = y; + } + + p = (float*) tabcos; + for ( i = -13*TABSTEP; i <= 13*TABSTEP; i++ ) { + xm = cos ((i-0.5)/TABSTEP); + x0 = cos ((i+0.0)/TABSTEP); + xp = cos ((i+0.5)/TABSTEP); + x = x0/2 + (xm + xp)/4; + y = xp - xm; + *p++ = x; + *p++ = y; + } + + p = (float*) tabsqrt_ex; + for ( i = 0; i < 255; i++ ) { + *(int*)&X = (i << 23); + *(int*)&Y = (i << 23) + (1<<23) - 1; + *p++ = sqrt(X); + } + *(int*)&X = (255 << 23) - 1; + *p++ = sqrt(X); + + p = (float*) tabsqrt_m; + for ( i = 1*TABSTEP; i <= 2*TABSTEP; i++ ) { + xm = sqrt ((i-0.5)/TABSTEP); + x0 = sqrt ((i+0.0)/TABSTEP); + xp = sqrt ((i+0.5)/TABSTEP); + x = x0/2 + (xm + xp)/4; + y = xp - xm; + *p++ = x; + *p++ = y; + } +} + +#endif Index: mppenc/branches/r2d/libmpcpsy/fastmath.h =================================================================== --- mppenc/branches/r2d/libmpcpsy/fastmath.h (revision 59) +++ mppenc/branches/r2d/libmpcpsy/fastmath.h (revision 59) @@ -0,0 +1,94 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#if 1 +# define ROUND32(x) ( floattmp = (x) + (int)0x00FD8000L, *(int*)(&floattmp) - (int)0x4B7D8000L ) +#else +# define ROUND32(x) ( (int) floor ((x) + 0.5) ) +#endif + +#ifdef FAST_MATH + +static __inline float +my_atan2 ( float x, float y ) +{ + float t; + int i; + float ret; + float floattmp; + + if ( (*(int*)&x & 0x7FFFFFFF) < (*(int*)&y & 0x7FFFFFFF) ) { + i = ROUND32 (t = TABSTEP * (x / y)); + ret = tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (t-i); + if ( *(int*)&y < 0 ) + ret = (float)(ret - M_PI); + } + else if ( *(int*)&x < 0) { + i = ROUND32 (t = TABSTEP * (y / x)); + 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)); + ret = + M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t); + } + else { + ret = 0.; + } + return ret; +} + + +static __inline float +my_cos ( float x ) +{ + float t; + int i; + float ret; + float floattmp; + + i = ROUND32 (t = TABSTEP * x); + ret = tabcos [13*TABSTEP+i][0] + tabcos [13*TABSTEP+i][1] * (t-i); + return ret; +} + + +static __inline int +my_ifloor ( float x ) +{ + x = x + (0x0C00000L + 0.500000001); + return *(int*)&x - 1262485505; +} + + +static __inline float +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 + ret *= tabsqrt_ex [ex]; + return ret; +} + +#endif Index: mppenc/branches/r2d/libmpcpsy/fft4g.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/fft4g.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/fft4g.c (revision 59) @@ -0,0 +1,670 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include "mppenc.h" + +/* F U N C T I O N S */ +static void makewt ( const int nw, int* ip, float* w ); +static void makect ( const int nc, int* ip, float* c ); +static __inline void bitrv2 ( const int n, int* ip, float* a ); // +static __inline void cftfsub ( const int n, float* a, float* w ); // +static __inline void rftfsub ( const int n, float* a, int nc, float* c ); // +static __inline void cft1st ( const int n, float* a, float* w ); // +static __inline void cftmdl_i386 ( const int n, const int l, float* a, float* w ); // 5648 +static __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 + +// generates lookup-tables +void +Generate_FFT_Tables ( const int n, int* ip, float* w ) +{ + int nw; + int nc; + + nw = n >> 2; + makewt ( nw, ip, w ); + + nc = n >> 2; + makect ( nc, ip, w + nw ); +} + + +// patched to only-forward +void +rdft ( const int n, float* a, int* ip, float* w ) +{ + float xi; + + ENTER(30); + if ( n > 4) { + bitrv2 ( n, ip + 2, a ); + cftfsub ( n, a, w ); + rftfsub ( n, a, ip[1], w + ip[0] ); + } + else if ( n == 4 ) { + cftfsub ( n, a, w ); + } + xi = a[0] - a[1]; + a[0] += a[1]; + a[1] = xi; + LEAVE(30); + return; +} + + +/* -------- initializing routines -------- */ +static void +makewt ( const int nw, int* ip, float* w ) +{ + int j; + int nwh; + float x; + float y; + double delta; + + ENTER(31); + ip[0] = nw; + ip[1] = 1; + if ( nw > 2 ) { + nwh = nw >> 1; + delta = (M_PI/4) / nwh; + w[0] = 1.; + w[1] = 0.; + w[nwh] = COS (delta * nwh); + w[nwh + 1] = w[nwh]; + if ( nwh > 2 ) { + for ( j = 2; j < nwh; j += 2 ) { + x = COS (delta * j); + y = SIN (delta * j); + w[j] = x; + w[j + 1] = y; + w[nw - j] = y; + w[nw - j + 1] = x; + } + bitrv2 ( nw, ip + 2, w ); + } + } + LEAVE(31); + return; +} + + +static void +makect ( const int nc, int* ip, float* c ) +{ + int j; + int nch; + double delta; + + ENTER(32); + ip[1] = nc; + if ( nc > 1 ) { + nch = nc >> 1; + delta = (M_PI/4) / nch; + c[0] = COS (delta * nch); + c[nch] = 0.5f * c[0]; + for ( j = 1; j < nch; j++ ) { + c[j] = 0.5f * COS (delta * j); + c[nc - j] = 0.5f * SIN (delta * j); + } + } + LEAVE(32); + return; +} + + +/* -------- 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; + + ENTER(33); + 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; + } + } + } + LEAVE(33); + 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; + + ENTER(34); + 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 ); + } + LEAVE(34); + 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; + + ENTER(35); + 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 ); + LEAVE(35); + return; +} + +extern void Cdecl cftmdl_3DNow_1 ( const int n, const int l, float* a, float* w ); +extern void 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; + + ENTER(36); + 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); + } + LEAVE(36); + + ENTER(39); + 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 ); + } + LEAVE(39); + 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; + + ENTER(36); + cftmdl_3DNow_1 (n,l,a,w); + LEAVE(36); + + ENTER(39); + 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 ); + } + LEAVE(39); + 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; + + ENTER(37); + 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 ); + LEAVE(37); + return; +} + +/* end of fft4g.c */ Index: mppenc/branches/r2d/libmpcpsy/fft_routines.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/fft_routines.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/fft_routines.c (revision 59) @@ -0,0 +1,337 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include "mppenc.h" + +#define CX0 -1. +#define CX1 0.5 + +#define SX1 -1. +#define SX2 (2./9/ 1) +#define SX3 (2./9/ 4) +#define SX4 (2./9/ 10) +#define SX5 (2./9/ 20) +#define SX6 (2./9/ 35) +#define SX7 (2./9/ 56) +#define SX8 (2./9/ 84) +#define SX9 (2./9/120) +#define SX10 (2./9/165) + + +#ifdef EXTRA_DECONV +# define DECONV \ + { \ + tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \ + aix[ 0] -= CX0*tmp; \ + aix[ 2] -= CX1*tmp; \ + tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \ + aix[ 3] -= SX1*tmp; \ + aix[ 5] -= SX2*tmp; \ + aix[ 7] -= SX3*tmp; \ + aix[ 9] -= SX4*tmp; \ + aix[11] -= SX5*tmp; \ + } +#elif 0 +# define DECONV \ + { \ + float A[20]; \ + int i; \ + memcpy (A, aix, 20*sizeof(aix)); \ + tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \ + aix[ 0] -= CX0*tmp; \ + aix[ 2] -= CX1*tmp; \ + tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \ + aix[ 3] -= SX1*tmp; \ + aix[ 5] -= SX2*tmp; \ + aix[ 7] -= SX3*tmp; \ + aix[ 9] -= SX4*tmp; \ + aix[11] -= SX5*tmp; \ + for ( i=0; i<10; i++) \ + printf ("%u%9.0f%7.0f%9.0f%7.0f\n",i, A[i+i], A[i+i+1], aix[i+i], aix[i+i+1] ); \ + } +#else +# define DECONV +#endif + + +/* V A R I A B L E S */ +static int ip [4096]; // bitinverse for maximum 2048 FFT +static float w [4096]; // butterfly-coefficient for maximum 2048 FFT +static float a [4096]; // holds real input for FFT +static float Hann_256 [ 256]; +static float Hann_1024 [1024]; +static float Hann_1600 [1600]; + + +////////////////////////////// +// +// 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 ( float* window, unsigned int size, float alpha ) +{ + double sumvalue = 0.; + double scale; + int i; + + scale = 0.25 / sqrt (size); + for ( i = 0; i < (int)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 < (int)size/2; i++ ) + window [size-1-i] = window [i] = /*sqrt*/ ( window [i] / sumvalue ) * scale; +} + +static void +CosWindow ( float* window, unsigned int size ) +{ + double x; + double scale; + int i; + + scale = 0.25 / sqrt (size); + for ( i = 0; i < (int)size/2; i++ ) { + x = cos ( (i+0.5) * (M_PI / size) ); + window [size/2-1-i] = window [size/2+i] = scale * x * x; + } +} + +static void +Window ( float* window, unsigned int size, float alpha ) +{ + if ( alpha < 0. ) + CosWindow ( window, size ) ; + else + KBDWindow ( window, size, alpha ); +} + + +/* F U N C T I O N S */ +// generates FFT lookup-tables +void +Init_FFT ( void ) +{ + int n; + double x; + double scale; + + // normalized hann functions + Window ( Hann_256 , 256, KBD1 ); + Window ( Hann_1024, 1024, KBD2 ); + scale = 0.25 / sqrt (2048.); + for ( n = 0; n < 800; n++ ) + x = cos ((n+0.5) * (M_PI/1600)), Hann_1600 [799-n] = Hann_1600 [800+n] = (float)(x * x * scale); + + Generate_FFT_Tables ( 2048, ip, w ); +} + +// input : Signal *x +// output: energy spectrum *erg +void +PowSpec256 ( const float* x, float* erg ) +{ + const float* win = Hann_256; + float* aix = a; + int i; + + ENTER(40); + // windowing + i = 256; + while (i--) + *aix++ = *x++ * *win++; + + // perform FFT + rdft ( 256, a, ip, w ); + + // calculate power + aix = a; // reset pointer + i = 128; + while (i--) { + *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; + aix += 2; + } + LEAVE(40); +} + +// input : Signal *x +// output: energy spectrum *erg +void +PowSpec1024 ( const float* x, float* erg ) +{ + const float* win = Hann_1024; + float* aix = a; + int i; + + ENTER(41); + i = 1024; // windowing + while (i--) + *aix++ = *x++ * *win++; + +// for (i=0; i<1024; i++) +// a[i] = Hann_1024[i] * ((i==0 ? 0 : i-512) + 1000); + + rdft ( 1024, a, ip, w ); // perform FFT + + aix = a; // calculate power + i = 512; + + + DECONV; +// for (i = 0; i <= 512; i++ ) +// printf ("%3u %12.6f %12.6f\n", i, a[i+i], a[i+i+1]); +// exit(1); + while (i--) { + *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; + aix += 2; + } + LEAVE(41); +} + +// input : Signal *x +// output: energy spectrum *erg +void +PowSpec2048 ( const float* x, float* erg ) +{ + const float* win = Hann_1600; + float* aix = a; + int i; + + ENTER(42); + // windowing (only 1600 samples available -> centered in 2048!) + memset ( a , 0, 224*sizeof(*a) ); + aix = a + 224; + i = 1600; + while (i--) + *aix++ = *x++ * *win++; + memset ( a+1824, 0, 224*sizeof(*a) ); + + rdft ( 2048, a, ip, w ); // perform FFT + + aix = a; // calculate power + i = 1024; + while (i--) { + *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; + aix += 2; + } + LEAVE(42); +} + +#include "fastmath.h" + +// input : Signal *x +// output: energy spectrum *erg and phase spectrum *phs +void +PolarSpec1024 ( const float* x, float* erg, float* phs ) +{ + const float* win = Hann_1024; + float* aix = a; + int i; + + ENTER(43); + i = 1024; // windowing + while (i--) + *aix++ = *x++ * *win++; + + rdft ( 1024, a, ip, w ); // perform FFT + + // calculate power and phase + aix = a; // reset pointer + i = 512; + while (i--) { + *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; + *phs++ = ATAN2F (aix[1], aix[0]); + aix += 2; + } + LEAVE(43); +} + +// input : logarithmized energy spectrum *cep +// output: Cepstrum *cep (in-place) +void +Cepstrum2048 ( float* cep, const int MaxLine ) +{ + float* aix = cep; + float* bix = cep + 2048; + int i; + + ENTER(44); + // generate real, even spectrum (symmetric around 1024, cep[2048-i] = cep[i]) + for ( i = 0; i < 1024; i++ ) + *bix-- = *aix++; + + // perform IFFT + rdft ( 2048, cep, ip, w ); + + // only real part as outcome (all even indexes of cep[]) + aix = cep; + bix = cep; + i = MaxLine + 1; + while (i--) { + *aix = *bix * (float) (0.9888 / 2048.); +// *aix = *bix * 0.0004828125f; + aix ++; + bix += 2; + } + LEAVE(44); +} Index: mppenc/branches/r2d/libmpcpsy/psy.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/psy.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/psy.c (revision 59) @@ -0,0 +1,1309 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +/* + * Prediction + * Short-Block-detection with smooth inset + * revise CalcMSThreshold + * /dev/audio for Windows too + * revise PNS/IS + * CVS with smoother inset + * several files per call + * revise ANS with changing SCFs + + * No IS + * PNS estimation very rough, also IS should be used to reduce data rate in the side channel + * ANS problems at Frame boundaries when resolution changes + * ANS problems at Subframe boundaries when SCF changes + * CVS+ with smoother transition + +---------------------------------------- + +Optimize Tabelle[18] (use second table) +CVS+ + +- ANS is disregarded during the search for the best Res +- ANS messes up if res changes (each 36 samples) and/or SCF changes (each 12 samples) +- PNS not in difference signal + +- implement IS in decoder +- Experimental Quantizer with complete energy preservation + - 1D, calculated + - 2D, calculated + - 2D, manually modified, coeffs set to 1.f + + */ + +#include "mppenc.h" + +/* V A R I A B L E S */ +/* further switches for the psymodel */ +unsigned int CVD_used; // global flag for ClearVoiceDetection +float varLtq; // variable threshold in quiet +unsigned int tmpMask_used; // global flag for temporal masking +float ShortThr; // Factor to calculate the masking threshold with transients +float minSMR; // minimum SMR for all subbands + +float a [PART_LONG]; +float b [PART_LONG]; +float c [PART_LONG]; +float d [PART_LONG]; // Integrations for tmpMask +static float Xsave_L [3 * 512]; +static float Xsave_R [3 * 512]; // FFT-Amplitudes L/R +static float Ysave_L [3 * 512]; +static float Ysave_R [3 * 512]; // FFT-Phases L/R +float T_L [PART_LONG]; +float T_R [PART_LONG]; // time-constants for tmpMask +float pre_erg_L[2][PART_SHORT]; +float pre_erg_R[2][PART_SHORT]; // Preecho-control short +float PreThr_L [PART_LONG]; +float PreThr_R [PART_LONG]; // for Pre-Echo-control L/R +float tmp_Mask_L [PART_LONG]; +float tmp_Mask_R [PART_LONG]; // for Post-Masking L/R +int Vocal_L [MAX_CVD_LINE + 4]; +int Vocal_R [MAX_CVD_LINE + 4]; // FFT-Line belongs to harmonic? + +/* F U N C T I O N S */ +// Resets Arrays +void +Init_Psychoakustik ( void ) +{ + int i; + + ENTER(200); + // generate FFT lookup-tables with largest FFT-size of 1024 + Init_FFT (); + + // setting pre-echo variables to Ltq + for ( i = 0; i < PART_LONG; i++ ) { + pre_erg_L [0][i/3] = pre_erg_R [0][i/3] = + pre_erg_L [1][i/3] = pre_erg_R [1][i/3] = + tmp_Mask_L [i] = tmp_Mask_R [i] = + PreThr_L [i] = PreThr_R [i] = partLtq [i]; + } + + // initializing arrays with zero + memset ( Xsave_L, 0, sizeof Xsave_L ); + memset ( Xsave_R, 0, sizeof Xsave_R ); + memset ( Ysave_L, 0, sizeof Ysave_L ); + memset ( Ysave_R, 0, sizeof Ysave_R ); + memset ( a, 0, sizeof a ); + memset ( b, 0, sizeof b ); + memset ( c, 0, sizeof c ); + memset ( d, 0, sizeof d ); + memset ( T_L, 0, sizeof T_L ); + memset ( T_R, 0, sizeof T_R ); + memset ( Vocal_L, 0, sizeof Vocal_L ); + memset ( Vocal_R, 0, sizeof Vocal_R ); + + LEAVE(200); + return; +} + + +// VBRmode 1: Adjustment of all SMRs via a factor (offset of SMRoffset dB) +// VBRmode 2: SMRs have a minimum of minSMR dB +static void +RaiseSMR_Signal ( const int MaxBand, float* signal, float tmp ) +{ + int Band; + float z = 0.; + + for ( Band = MaxBand; Band >= 0; Band-- ) { + if ( z < signal [Band] ) z = signal [Band]; + if ( z > tmp ) z = tmp; + if ( signal [Band] < z ) signal [Band] = z; + } +} + + +void +RaiseSMR ( const int MaxBand, SMRTyp* smr ) +{ + float tmp = POW10 ( 0.1 * minSMR ); + + ENTER(201); + RaiseSMR_Signal ( MaxBand, smr->L, tmp ); + RaiseSMR_Signal ( MaxBand, smr->R, tmp ); + RaiseSMR_Signal ( MaxBand, smr->M, tmp ); + RaiseSMR_Signal ( MaxBand, smr->S, 0.5 * tmp ); + + LEAVE(201); + return; +} + +// input : *smr +// output: *smr, *ms, *x (only the entries for L/R contain relevant data) +// Check if either M/S- or L/R-coding has a lower perceptual entropy +// Choose the better mode, copy the appropriate data into the +// arrays that belong to L and R and set the ms-Flag accordingly. +void +MS_LR_Entscheidung ( const int MaxBand, unsigned char* ms, SMRTyp* smr, SubbandFloatTyp* x ) +{ + int Band; + int n; + float PE_MS; + float PE_LR; + float tmpM; + float tmpS; + float* l; + float* r; + + ENTER(202); + + for ( Band = 0; Band <= MaxBand; Band++ ) { // calculate perceptual entropy + PE_LR = PE_MS = 1.f; + if (smr->L[Band] > 1.) PE_LR *= smr->L[Band]; + if (smr->R[Band] > 1.) PE_LR *= smr->R[Band]; + if (smr->M[Band] > 1.) PE_MS *= smr->M[Band]; + if (smr->S[Band] > 1.) PE_MS *= smr->S[Band]; + + if ( PE_MS < PE_LR ) { + ms[Band] = 1; + + // calculate M/S-signal and copies it to L/R-array + l = x[Band].L; + r = x[Band].R; + for ( n = 0; n < 36; n++, l++, r++ ) { + tmpM = (*l + *r) * 0.5f; + tmpS = (*l - *r) * 0.5f; + *l = tmpM; + *r = tmpS; + } + + // copy M/S - SMR to L/R-fields + smr->L[Band] = smr->M[Band]; + smr->R[Band] = smr->S[Band]; + } + else { + ms[Band] = 0; + } + } + + LEAVE(202); + return; +} + +// input : FFT-spectrums *spec0 und *spec1 +// output: energy in the individual subbands *erg0 and *erg1 +// With Butfly[], you can calculate the results of aliasing during calculation +// of subband energy from the FFT-spectrums. +static void +SubbandEnergy ( const int MaxBand, + float* erg0, + float* erg1, + const float* spec0, + const float* spec1 ) +{ + int n; + int k; + int alias; + float tmp0; + float tmp1; + + ENTER(203); + + // Is this here correct for FFT-based data or is this calculation rule only for MDCTs??? + + for ( k = 0; k <= MaxBand; k++ ) { // subband index + tmp0 = tmp1 = 0.f; + for ( n = 0; n < 16; n++, spec0++, spec1++ ) { // spectral index + tmp0 += *spec0; + tmp1 += *spec1; + + // Consideration of Aliasing between the subbands + if ( n < +sizeof(Butfly)/sizeof(*Butfly) && k != 0 ) { + alias = -1 - (n<<1); + tmp0 += Butfly [n] * (spec0[alias] - *spec0); + tmp1 += Butfly [n] * (spec1[alias] - *spec1); + } + else if ( n > 15-sizeof(Butfly)/sizeof(*Butfly) && k != 31 ) { + alias = 31 - (n<<1); + tmp0 += Butfly [15-n] * (spec0[alias] - *spec0); + tmp1 += Butfly [15-n] * (spec1[alias] - *spec1); + } + } + *erg0++ = tmp0; + *erg1++ = tmp1; + } + + LEAVE(203); + return; +} + +// input : FFT-Spectrums *spec0 and *spec1 +// output: energy in the individual partitions *erg0 and *erg1 +static void +PartitionEnergy ( float* erg0, + float* erg1, + const float* spec0, + const float* spec1 ) +{ + unsigned int n; + unsigned int k; + float e0; + float e1; + + ENTER(204); + +#if 000000 + for ( n = 0; n < PART_LONG; n++ ) { + k = wh[n] - wl[n]; + e0 = *spec0++; + e1 = *spec1++; + while ( k-- ) { + e0 += *spec0++; + e1 += *spec1++; + } + *erg0++ = e0; + *erg1++ = e1; + } +#else + n = 0; + + for ( ; n < 23; n++ ) { // 11 or 23 + k = wh[n] - wl[n]; + e0 = *spec0++; + e1 = *spec1++; + while ( k-- ) { + e0 += *spec0++; + e1 += *spec1++; + } + *erg0++ = e0; + *erg1++ = e1; + } + + for ( ; n < 48; n++ ) { // 37 ... 46, 48, 57 + k = wh[n] - wl[n]; + e0 = sqrt (*spec0++); + e1 = sqrt (*spec1++); + while ( k-- ) { + e0 += sqrt (*spec0++); + e1 += sqrt (*spec1++); + } + *erg0++ = e0*e0 * iw[n]; + *erg1++ = e1*e1 * iw[n]; + } + + for ( ; n < PART_LONG; n++ ) { + k = wh[n] - wl[n]; + e0 = *spec0++; + e1 = *spec1++; + while ( k-- ) { + e0 += *spec0++; + e1 += *spec1++; + } + *erg0++ = e0; + *erg1++ = e1; + } + + +#endif + + LEAVE(204); + return; +} + + +// input : FFT-Spectrums *spec0, *spec1 and unpredictability *cw0 and *cw1 +// output: weighted energy in the individual partitions *erg0, *erg1 +static void +WeightedPartitionEnergy ( float* erg0, + float* erg1, + const float* spec0, + const float* spec1, + const float* cw0, + const float* cw1 ) +{ + unsigned int n; + unsigned int k; + float e0; + float e1; + + ENTER(205); + +#if 000000 + for ( n = 0; n < PART_LONG; n++ ) { + e0 = *spec0++ * *cw0++; + e1 = *spec1++ * *cw1++; + k = wh[n] - wl[n]; + while ( k-- ) { + e0 += *spec0++ * *cw0++; + e1 += *spec1++ * *cw1++; + } + *erg0++ = e0; + *erg1++ = e1; + } +#else + n = 0; + + for ( ; n < 23; n++ ) { + e0 = *spec0++ * *cw0++; + e1 = *spec1++ * *cw1++; + k = wh[n] - wl[n]; + while ( k-- ) { + e0 += *spec0++ * *cw0++; + e1 += *spec1++ * *cw1++; + } + *erg0++ = e0; + *erg1++ = e1; + } + + for ( ; n < 48; n++ ) { + e0 = sqrt (*spec0++ * *cw0++); + e1 = sqrt (*spec1++ * *cw1++); + k = wh[n] - wl[n]; + while ( k-- ) { + e0 += sqrt (*spec0++ * *cw0++); + e1 += sqrt (*spec1++ * *cw1++); + } + *erg0++ = e0*e0 * iw[n]; + *erg1++ = e1*e1 * iw[n]; + } + + for ( ; n < PART_LONG; n++ ) { + e0 = *spec0++ * *cw0++; + e1 = *spec1++ * *cw1++; + k = wh[n] - wl[n]; + while ( k-- ) { + e0 += *spec0++ * *cw0++; + e1 += *spec1++ * *cw1++; + } + *erg0++ = e0; + *erg1++ = e1; + } +#endif + + LEAVE(205); + return; +} + +// input : masking thresholds, first half of the arrays *shaped0 and *shaped1 +// output: masking thresholds, second half of the arrays *shaped0 and *shaped1 +// Considering the result of aliasing via InvButfly[] +// The input *thr0, *thr1 is gathered via address calculation from *shaped0, *shaped1 + +static void +AdaptThresholds ( const int MaxLine, float* shaped0, float* shaped1 ) +{ + int n; + int mod; + int alias; + float tmp; + const float* invb = InvButfly; + const float* thr0 = shaped0 - 512; + const float* thr1 = shaped1 - 512; + float tmp0; + float tmp1; + + ENTER(206); + + // should be able to optimize it with coasting. [ 9 ] + n * [ 7 + 7 + 2 ] + [ 7 ] + // Schleife Schl Schl Ausr Schleife + for ( n = 0; n < MaxLine; n++, thr0++, thr1++ ) { + mod = n & 15; // n%16 + tmp0 = *thr0; + tmp1 = *thr1; + + if ( mod < +sizeof(InvButfly)/sizeof(*InvButfly) && n > 12 ) { + alias = -1 - (mod<<1); + tmp = thr0[alias] * invb[mod]; + if ( tmp < tmp0 ) tmp0 = tmp; + tmp = thr1[alias] * invb[mod]; + if ( tmp < tmp1 ) tmp1 = tmp; + } + else if ( mod > 15-sizeof(InvButfly)/sizeof(*InvButfly) && n < 499 ) { + alias = 31 - (mod<<1); + tmp = thr0[alias] * invb[15-mod]; + if ( tmp < tmp0 ) tmp0 = tmp; + tmp = thr1[alias] * invb[15-mod]; + if ( tmp < tmp1 ) tmp1 = tmp; + } + *shaped0++ = tmp0; + *shaped1++ = tmp1; + } + + LEAVE(206); + return; +} + +#include "fastmath.h" + +// input : current spectrum in the form of power *spec and phase *phase, +// the last two earlier spectrums are at position +// 512 and 1024 of the corresponding Input-Arrays. +// Array *vocal, which can mark an FFT_Linie as harmonic +// output: current amplitude *amp and unpredictability *cw +static void +CalcUnpred ( const int MaxLine, + const float* spec, + const float* phase, + const int* vocal, + float* amp0, + float* phs0, + float* cw ) +{ + int n; + float amp; + float tmp; +#define amp1 ((amp0) + 512) // amp[ 512...1023] contains data of frame-1 +#define amp2 ((amp0) + 1024) // amp[1024...1535] contains data of frame-2 +#define phs1 ((phs0) + 512) // phs[ 512...1023] contains data of frame-1 +#define phs2 ((phs0) + 1024) // phs[1024...1535] contains data of frame-2 + + ENTER(207); + + for ( n = 0; n < MaxLine; n++ ) { + tmp = COSF ((phs0[n] = phase[n]) - 2*phs1[n] + phs2[n]); // copy phase to output-array, predict phase and calculate predictive error + amp0[n] = SQRTF (spec[n]); // calculate and set amplitude + amp = 2*amp1[n] - amp2[n]; // predict amplitude + + // calculate unpredictability + cw[n] = SQRTF (spec[n] + amp * (amp - 2*amp0[n] * tmp)) / (amp0[n] + FABS(amp)); + } + + // postprocessing of harmonic FFT-lines (*cw is set to CVD_UNPRED) + if ( CVD_used && vocal != NULL ) { + for ( n = 0; n < MAX_CVD_LINE; n++, cw++, vocal++ ) + if ( *vocal != 0 && *cw > CVD_UNPRED * 0.01 * *vocal ) + *cw = CVD_UNPRED * 0.01 * *vocal; + } + + LEAVE(207); + return; +} +#undef amp1 +#undef amp2 +#undef phs1 +#undef phs2 + + +// input : Energy *erg, calibrated energy *werg +// output: spread energy *res, spread weighted energy *wres +// SPRD describes the spreading function as calculated in psy_tab.c +static void +SpreadingSignal ( const float* erg, const float* werg, float* res, float* wres ) +{ + int n; + int k; + int start; + int stop; + const float* sprd; + float e; + float ew; + + ENTER(208); + + for (k=0; k +static float +AdaptLtq ( const float* erg0, const float* erg1 ) +{ + static float loud = 0.f; + float* weight = Loudness; + float sum = 0.f; + int n; + + // calculate loudness + for ( n = 0; n < PART_LONG; n++ ) + sum += (*erg0++ + *erg1++) * *weight++; + + // Utilization of the time constants (fast drop of Ltq T=5, slow rise of Ltq T=20) + //loud = (sum < loud) ? (4 * sum + loud)*0.2f : (19 * loud + sum)*0.05f; + loud = 0.98 * loud + 0.02 * (0.5 * sum); + + // calculate dynamic offset for threshold in quiet, 0...+20 dB, at 96 dB loudness, an offset of 20 dB is assumed + return 1.f + varLtq * loud * 5.023772e-08f; +} + +// input : simultaneous masking threshold *frqthr, +// previous masking threshold *tmpthr, +// Integrations *a (short-time) and *b (long-time) +// output: tracked Integrations *a and *b, time constant *tau +static void +CalcTemporalThreshold ( float* a, float* b, float* tau, float* frqthr, float* tmpthr ) +{ + int n; + float tmp; + + ENTER(220); + + for ( n = 0; n < PART_LONG; n++ ) { + // following calculations relative to threshold in quiet + frqthr[n] *= invLtq[n]; + tmpthr[n] *= invLtq[n]; + + // new post-masking 'tmp' via time constant tau, if old post-masking > Ltq (=1) + tmp = tmpthr[n] > 1.f ? POW ( tmpthr[n], tau[n] ) : 1.f; + + // calculate time constant for post-masking in next frame, + // if new time constant has to be calculated (new tmpMask < frqMask) + a[n] += 0.5f * (frqthr[n] - a[n]); // short time integrator + b[n] += 0.15f * (frqthr[n] - b[n]); // long time integrator + if (tmp < frqthr[n]) + tau[n] = a[n] <= b[n] ? 0.8f : 0.2f + b[n] / a[n] * 0.6f; + + // use post-masking of (Re-Normalization) + tmpthr[n] = maxf (frqthr[n], tmp) * partLtq[n]; + } + + LEAVE(220); + return; +} + +// input : L/R-Masking thresholds in Partitions *thrL, *thrR +// L/R-Subband energies *ergL, *ergR +// M/S-Subband energies *ergM, *ergS +// output: M/S-Masking thresholds in Partitions *thrM, *thrS +static void +CalcMSThreshold ( const float* const ergL, + const float* const ergR, + const float* const ergM, + const float* const ergS, + float* const thrL, + float* const thrR, + float* const thrM, + float* const thrS ) +{ + int n; + float norm; + float tmp; + + // All hardcoded numbers here should be pulled from somewhere, + // the "4.", the -2 dB, the 0.0625 and the 0.9375, as well as all bands where this is done + + for ( n = 0; n < PART_LONG; n++ ) { + // estimate M/S thresholds out of L/R thresholds and M/S and L/R energies + thrS[n] = thrM[n] = maxf (ergM[n], ergS[n]) / maxf (ergL[n], ergR[n]) * minf (thrL[n], thrR[n]); + + switch ( MS_Channelmode ) { // preserve 'near-mid' signal components + case 3: + if ( n > 0 ) { + double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n]; + double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n]; + if ( ratioMS < ratioLR ) { // MS + if ( ergM[n] > ergS[n] ) + thrS[n] = thrL[n] = thrR[n] = 1.e18f; + else + thrM[n] = thrL[n] = thrR[n] = 1.e18f; + } + else { // LR + if ( ergL[n] > ergR[n] ) + thrR[n] = thrM[n] = thrS[n] = 1.e18f; + else + thrL[n] = thrM[n] = thrS[n] = 1.e18f; + } + } + break; + case 4: + if ( n > 0 ) { + double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n]; + double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n]; + if ( ratioMS < ratioLR ) { // MS + if ( ergM[n] > ergS[n] ) + thrS[n] = 1.e18f; + else + thrM[n] = 1.e18f; + } + else { // LR + if ( ergL[n] > ergR[n] ) + thrR[n] = 1.e18f; + else + thrL[n] = 1.e18f; + } + } + break; + case 5: + thrS[n] *= 2.; // +3 dB + break; + case 6: + break; + default: + fprintf ( stderr, "Unknown stereo mode\n"); + case 10: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.70794578f * iw[n]; // -1.5 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = MS2SPAT1 * thrS[n] + (1.f-MS2SPAT1) * tmp; // raises masking threshold by up to 3 dB + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = MS2SPAT1 * thrM[n] + (1.f-MS2SPAT1) * tmp; + } + } + break; + case 11: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.63095734f * iw[n]; // -2.0 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = MS2SPAT2 * thrS[n] + (1.f-MS2SPAT2) * tmp; // raises masking threshold by up to 6 dB + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = MS2SPAT2 * thrM[n] + (1.f-MS2SPAT2) * tmp; + } + } + break; + case 12: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = MS2SPAT3 * thrS[n] + (1.f-MS2SPAT3) * tmp; // raises masking threshold by up to 9 dB + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = MS2SPAT3 * thrM[n] + (1.f-MS2SPAT3) * tmp; + } + } + break; + case 13: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = MS2SPAT4 * thrS[n] + (1.f-MS2SPAT4) * tmp; // raises masking threshold by up to 12 dB + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = MS2SPAT4 * thrM[n] + (1.f-MS2SPAT4) * tmp; + } + } + break; + case 15: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = tmp; // raises masking threshold by up to +oo dB an + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = tmp; + } + } + break; + case 22: + if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB + norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth + if ( ergM[n] > ergS[n] ) { + tmp = ergS[n] * norm; + if ( thrS[n] > tmp ) + thrS[n] = maxf (tmp, ergM[n]*iw[n]*0.025); // +/- 1.414° + } else if ( ergS[n] > ergM[n] ) { + tmp = ergM[n] * norm; + if ( thrM[n] > tmp ) + thrM[n] = maxf (tmp, ergS[n]*iw[n]*0.025); // +/- 1.414° + } + } + break; + } + } + + return; +} + +// input : Masking thresholds in Partitions *partThr0, *partThr1 +// level of threshold in quiet *ltq in FFT-resolution +// output: Masking thresholds in FFT-resolution *thr0, *thr1 +// inline, because it's called 4x +static void +ApplyLtq ( float* thr0, + float* thr1, + const float* partThr0, + const float* partThr1, + const float AdaptedLTQ, + int MSflag ) +{ + int n; + int k; + float ltq; + float tmp; + float ms = MSflag ? 0.125f * AdaptedLTQ : 0.25f * AdaptedLTQ ; + + for ( n = 0; n < PART_LONG; n++ ) { + for ( k = wl[n]; k <= wh[n]; k++, thr0++, thr1++ ) { // threshold in quiet (Partition) +#if 0 + ltq = AdaptedLTQ * fftLtq [k]; + *thr0 = maxf ( partThr0 [n], ltq ); + *thr1 = maxf ( partThr1 [n], ltq ); +#else + // Applies a much more gentle ATH rolloff + 6 dB more dynamic + ltq = sqrt (ms * fftLtq [k]); + tmp = sqrt (partThr0 [n]) + ltq; + *thr0 = tmp * tmp; + tmp = sqrt (partThr1 [n]) + ltq; + *thr1 = tmp * tmp; +#endif + } + } + return; +} + +// input : Subband energies *erg0, *erg1 +// Masking thresholds in FFT-resolution *thr0, *thr1 +// output: SMR per Subband *smr0, *smr1 +static void +CalculateSMR ( const int MaxBand, + const float* erg0, + const float* erg1, + const float* thr0, + const float* thr1, + float* smr0, + float* smr1 ) +{ + int n; + int k; + float tmp0; + float tmp1; + + // calculation of the masked thresholds in the subbands + for (n = 0; n <= MaxBand; n++ ) { + tmp0 = *thr0++; + tmp1 = *thr1++; + for (k=1; k<16; ++k, ++thr0, ++thr1) { + if (*thr0 < tmp0) tmp0 = *thr0; + if (*thr1 < tmp1) tmp1 = *thr1; + } + *smr0++ = 0.0625f * *erg0++ / tmp0; + *smr1++ = 0.0625f * *erg1++ / tmp1; + } + + return; +} + +// input : energy spectrums erg[4][128] (4 delayed FFTs) +// Energy of the last short block *preerg in short partitions +// PreechoFac declares allowed traved of the masking threshold +// output: masking threshold *thr in short partitions +// Energy of the last short block *preerg in short partitions +#if 0 +static void +CalcShortThreshold ( const float erg [] [128], + const float PreechoFac, + float* thr, + float preerg[2][PART_SHORT], + int* transient ) +{ + const int* lo = wl_short; // lower FFT-index + const int* hi = wh_short; // upper FFT-index + const float* iwidth = iw_short; // inverse partition-width + int k; + int n; + int m; + float tmp; + float enrg; + float th; + const float* ep; + + for ( k = 0; k < PART_SHORT; k++, lo++, hi++ ) { + transient[k] = 0; + th = 1.e20f; + for ( n = 0; n < 4; n++ ) { + ep = erg[n] + *lo; + m = *hi - *lo; + enrg = *ep++; + while (m--) + enrg += *ep++; + + // preecho prevention + tmp = enrg; + if (preerg[0][k] < enrg) + enrg = preerg[0][k]; + preerg[0][k] = tmp; + + // is signal transient? + if (tmp > TransDetect*enrg) transient[k] = 1; + + // assume short threshold = engr*PreechoFac + th = minf (th, enrg*PreechoFac); + } + thr[k] = th * *iwidth++; + } + + return; +} +#else +static void +CalcShortThreshold ( const float erg [4] [128], + const float ShortThr, + float* thr, + float old_erg [2][PART_SHORT], + int* transient ) +{ + const int* index_lo = wl_short; // lower FFT-index + const int* index_hi = wh_short; // upper FFT-index + const float* iwidth = iw_short; // inverse partition-width + int k; + int n; + int m; + float new_erg; + float th; + const float* ep; + + for ( k = 0; k < PART_SHORT; k++ ) { + transient [k] = 0; + th = old_erg [0][k]; + for ( n = 0; n < 4; n++ ) { + ep = erg[n] + index_lo [k]; + m = index_hi [k] - index_lo [k]; + + new_erg = *ep++; + while (m--) + new_erg += *ep++; // e = Short_Partition-energy in piece n + + if ( new_erg > old_erg [0][k] ) { // bigger than the old? + + if ( new_erg > old_erg [0][k] * TransDetect || + new_erg > old_erg [1][k] * TransDetect*2 ) // is signal transient? + transient [k] = 1; + } + else { + th = minf ( th, new_erg ); // assume short threshold = engr*PreechoFac + } + + old_erg [1][k] = old_erg [0][k]; + old_erg [0][k] = new_erg; // save the current one + } + thr [k] = th * ShortThr * *iwidth++; // pull out and multiply only when transient[k]=1 + } + + return; +} + +#endif + +// input : previous simultaneous masking threshold *preThr, +// current simultaneous masking threshold *simThr +// output: update of *preThr for next call, +// current masking threshold *partThr +static void +PreechoControl ( float* partThr0, + float* preThr0, + const float* simThr0, + float* partThr1, + float* preThr1, + const float* simThr1 ) +{ + int n; + + for ( n = 0; n < PART_LONG; n++ ) { + *partThr0++ = minf ( *simThr0, *preThr0 * PREFAC_LONG); + *partThr1++ = minf ( *simThr1, *preThr1 * PREFAC_LONG); + *preThr0++ = *simThr0++; + *preThr1++ = *simThr1++; + } + return; +} + + +void +TransientenCalc ( int* T, + const int* TL, + const int* TR ) +{ + int i; + int x1; + int x2; + + memset ( T, 0, 32*sizeof(*T) ); + + for ( i = 0; i < PART_SHORT; i++ ) + if ( TL[i] || TR[i] ) { + x1 = wl_short[i] >> 2; + x2 = wh_short[i] >> 2; + while ( x1 <= x2 ) + T [x1++] = 1; + } +} + + +// input : PCM-Data *data +// output: SMRs for the input data +SMRTyp +Psychoakustisches_Modell ( const int MaxBand, const PCMDataTyp* data, int* TransientL, int* TransientR ) +{ + float Xi_L[32], Xi_R[32]; // acoustic pressure per Subband L/R + float Xi_M[32], Xi_S[32]; // acoustic pressure per Subband M/S + float cw_L[512], cw_R[512]; // unpredictability (only L/R) + float erg0[512], erg1[512]; // holds energy spectrum of long FFT + float phs0[512], phs1[512]; // holds phase spectrum of long FFT + float Thr_L[2*512], Thr_R[2*512]; // masking thresholds L/R, second half for triangle swap + float Thr_M[2*512], Thr_S[2*512]; // masking thresholds M/S, second half for triangle swap + float F_256[4][128]; // holds energies of short FFTs (L/R only) + float Xerg[1024]; // holds energy spectrum of very long FFT + float Ls_L[PART_LONG], Ls_R[PART_LONG]; // acoustic pressure in Partition L/R + float Ls_M[PART_LONG], Ls_S[PART_LONG]; // acoustic pressure per each partition M/S + float PartThr_L[PART_LONG], PartThr_R[PART_LONG]; // masking thresholds L/R (Partition) + float PartThr_M[PART_LONG], PartThr_S[PART_LONG]; // masking thresholds M/S (Partition) + float sim_Mask_L[PART_LONG], sim_Mask_R[PART_LONG]; // simultaneous masking (only L/R) + float clow_L[PART_LONG], clow_R[PART_LONG]; // spread, weighted energy (only L/R) + float cLs_L[PART_LONG], cLs_R[PART_LONG]; // weighted partition energy (only L/R) + float shortThr_L[PART_SHORT],shortThr_R[PART_SHORT]; // threshold for short FFT (only L/R) + int n; + int MaxLine = (MaxBand+1)*16; // set FFT-resolution according to MaxBand + SMRTyp SMR0; + SMRTyp SMR1; // holds SMR's for first and second Analysis + int isvoc_L; + int isvoc_R; + float factorLTQ = 1.f; // Offset after variable LTQ + + ENTER(50); + // 'ClearVocalDetection'-Process + if ( CVD_used ) { + memset ( Vocal_L, 0, sizeof Vocal_L ); + memset ( Vocal_R, 0, sizeof Vocal_R ); + + // left channel + PowSpec2048 ( &data->L[0], Xerg ); + isvoc_L = CVD2048 ( Xerg, Vocal_L ); + // right channel + PowSpec2048 ( &data->R[0], Xerg ); + isvoc_R = CVD2048 ( Xerg, Vocal_R ); + } + + // calculation of the spectral energy via FFT + PolarSpec1024 ( &data->L[0], erg0, phs0 ); // left + PolarSpec1024 ( &data->R[0], erg1, phs1 ); // right + + // calculation of the acoustic pressures per each subband for L/R-signals + SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 ); + + // calculation of the acoustic pressures per each partition + PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 ); + + // calculate the predictability of the signal + // left + memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) ); + memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) ); + CalcUnpred ( MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L ); + // right + memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) ); + memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) ); + CalcUnpred ( MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R ); + + // calculation of the weighted acoustic pressures per each partition + WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R ); + + // Spreading Signal & weighted unpredictability-signal + // left + memset ( clow_L , 0, sizeof clow_L ); + memset ( sim_Mask_L, 0, sizeof sim_Mask_L ); + SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L ); + // right + memset ( clow_R , 0, sizeof clow_R ); + memset ( sim_Mask_R, 0, sizeof sim_Mask_R ); + SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R ); + + // Offset depending on tonality + ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R ); + + // handling of transient signals + // calculate four short FFTs (left) + PowSpec256 ( &data->L[ 0+SHORTFFT_OFFSET], F_256[0] ); + PowSpec256 ( &data->L[144+SHORTFFT_OFFSET], F_256[1] ); + PowSpec256 ( &data->L[288+SHORTFFT_OFFSET], F_256[2] ); + PowSpec256 ( &data->L[432+SHORTFFT_OFFSET], F_256[3] ); + // calculate short Threshold + CalcShortThreshold ( F_256, ShortThr, shortThr_L, pre_erg_L, TransientL ); + + // calculate four short FFTs (right) + PowSpec256 ( &data->R[ 0+SHORTFFT_OFFSET], F_256[0] ); + PowSpec256 ( &data->R[144+SHORTFFT_OFFSET], F_256[1] ); + PowSpec256 ( &data->R[288+SHORTFFT_OFFSET], F_256[2] ); + PowSpec256 ( &data->R[432+SHORTFFT_OFFSET], F_256[3] ); + // calculate short Threshold + CalcShortThreshold ( F_256, ShortThr, shortThr_R, pre_erg_R, TransientR ); + + // dynamic adjustment of the threshold in quiet to the loudness of the current sequence + if ( varLtq > 0. ) + factorLTQ = AdaptLtq ( Ls_L, Ls_R ); + + // utilization of the temporal post-masking + if ( tmpMask_used ) { + CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L ); + CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R ); + memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L ); + memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R ); + } + + // transient signal? + for ( n = 0; n < PART_SHORT; n++ ) { + if ( TransientL [n] ) { + sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] ); + sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] ); + sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] ); + } + if ( TransientR[n] ) { + sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] ); + sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] ); + sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] ); + } + } + + // Pre-Echo control + PreechoControl ( PartThr_L, PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R ); + + // utilization of the threshold in quiet + ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 ); + + // Consideration of aliasing between the subbands (noise is smeared) + // In: Thr[0..511], Out: Thr[512...1023] + AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 ); + memmove ( Thr_L, Thr_L+512, 512*sizeof(float) ); + memmove ( Thr_R, Thr_R+512, 512*sizeof(float) ); + + // calculation of the Signal-to-Mask-Ratio + CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR0.L, SMR0.R ); + + /***************************************************************************************/ + /***************************************************************************************/ + if ( MS_Channelmode > 0 ) { + // calculation of the spectral energy via FFT + PowSpec1024 ( &data->M[0], erg0 ); // mid + PowSpec1024 ( &data->S[0], erg1 ); // side + + // calculation of the acoustic pressures per each subband for M/S-signals + SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 ); + + // calculation of the acoustic pressures per each partition + PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 ); + + // calculate masking thresholds for M/S + CalcMSThreshold ( Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S ); + ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 ); + + // Consideration of aliasing between the subbands (noise is smeared) + // In: Thr[0..511], Out: Thr[512...1023] + AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 ); + memmove ( Thr_M, Thr_M+512, 512*sizeof(float) ); + memmove ( Thr_S, Thr_S+512, 512*sizeof(float) ); + + // calculation of the Signal-to-Mask-Ratio + CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR0.M, SMR0.S ); + } + + if ( NS_Order > 0 ) { // providing the Noise Shaping thresholds + memcpy ( ANSspec_L, Thr_L, sizeof ANSspec_L ); + memcpy ( ANSspec_R, Thr_R, sizeof ANSspec_R ); + memcpy ( ANSspec_M, Thr_M, sizeof ANSspec_M ); + memcpy ( ANSspec_S, Thr_S, sizeof ANSspec_S ); + } + /***************************************************************************************/ + /***************************************************************************************/ + + // + //-------- second model calculation via shifted FFT ------------------------ + // + // calculation of the spectral power via FFT + PolarSpec1024 ( &data->L[576], erg0, phs0 ); // left + PolarSpec1024 ( &data->R[576], erg1, phs1 ); // right + + // calculation of the acoustic pressures per each subband for L/R-signals + SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 ); + + // calculation of the acoustic pressures per each partition + PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 ); + + // calculate the predictability of the signal + // left + memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) ); + memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) ); + CalcUnpred ( MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L ); + // right + memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) ); + memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) ); + CalcUnpred ( MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R ); + + // calculation of the weighted acoustic pressure per each partition + WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R ); + + // Spreading Signal & weighted unpredictability-signal + // left + memset ( clow_L , 0, sizeof clow_L ); + memset ( sim_Mask_L, 0, sizeof sim_Mask_L ); + SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L ); + // right + memset ( clow_R , 0, sizeof clow_R ); + memset ( sim_Mask_R, 0, sizeof sim_Mask_R ); + SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R ); + + // Offset depending on tonality + ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R ); + + // Handling of transient signals + // calculate four short FFTs (left) + PowSpec256 ( &data->L[ 576+SHORTFFT_OFFSET], F_256[0] ); + PowSpec256 ( &data->L[ 720+SHORTFFT_OFFSET], F_256[1] ); + PowSpec256 ( &data->L[ 864+SHORTFFT_OFFSET], F_256[2] ); + PowSpec256 ( &data->L[1008+SHORTFFT_OFFSET], F_256[3] ); + // calculate short Threshold + CalcShortThreshold ( F_256, ShortThr, shortThr_L, pre_erg_L, TransientL ); + + // calculate four short FFTs (right) + PowSpec256 ( &data->R[ 576+SHORTFFT_OFFSET], F_256[0] ); + PowSpec256 ( &data->R[ 720+SHORTFFT_OFFSET], F_256[1] ); + PowSpec256 ( &data->R[ 864+SHORTFFT_OFFSET], F_256[2] ); + PowSpec256 ( &data->R[1008+SHORTFFT_OFFSET], F_256[3] ); + // calculate short Threshold + CalcShortThreshold ( F_256, ShortThr, shortThr_R, pre_erg_R, TransientR ); + + // dynamic adjustment of threshold in quiet to loudness of the current sequence + if ( varLtq > 0. ) + factorLTQ = AdaptLtq ( Ls_L, Ls_R ); + + // utilization of temporal post-masking + if (tmpMask_used) { + CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L ); + CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R ); + memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L ); + memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R ); + } + + // transient signal? + for ( n = 0; n < PART_SHORT; n++ ) { + if ( TransientL[n] ) { + sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] ); + sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] ); + sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] ); + } + if ( TransientR[n] ) { + sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] ); + sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] ); + sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] ); + } + } + + // Pre-Echo control + PreechoControl ( PartThr_L, PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R ); + + // utilization of threshold in quiet + ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 ); + + // Consideration of aliasing between the subbands (noise is smeared) + // In: Thr[0..511], Out: Thr[512...1023] + AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 ); + memmove ( Thr_L, Thr_L+512, 512*sizeof(float) ); + memmove ( Thr_R, Thr_R+512, 512*sizeof(float) ); + + // calculation of the Signal-to-Mask-Ratio + CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR1.L, SMR1.R ); + + /***************************************************************************************/ + /***************************************************************************************/ + if ( MS_Channelmode > 0 ) { + // calculation of the spectral energy via FFT + PowSpec1024 ( &data->M[576], erg0 ); // mid + PowSpec1024 ( &data->S[576], erg1 ); // side + + // calculation of the acoustic pressure per each subband for M/S-signals + SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 ); + + // calculation of the acoustic pressure per each partition + PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 ); + + // calculate masking thresholds for M/S + CalcMSThreshold ( Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S ); + ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 ); + + // Consideration of aliasing between the subbands (noise is smeared) + // In: Thr[0..511], Out: Thr[512...1023] + AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 ); + memmove ( Thr_M, Thr_M+512, 512*sizeof(float) ); + memmove ( Thr_S, Thr_S+512, 512*sizeof(float) ); + + // calculation of the Signal-to-Mask-Ratio + CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR1.M, SMR1.S ); + } + /***************************************************************************************/ + /***************************************************************************************/ + + if ( NS_Order > 0 ) { + for ( n = 0; n < MAX_ANS_LINES; n++ ) { // providing Noise Shaping thresholds + ANSspec_L [n] = minf ( ANSspec_L [n], Thr_L [n] ); + ANSspec_R [n] = minf ( ANSspec_R [n], Thr_R [n] ); + ANSspec_M [n] = minf ( ANSspec_M [n], Thr_M [n] ); + ANSspec_S [n] = minf ( ANSspec_S [n], Thr_S [n] ); + } + } + + for ( n = 0; n <= MaxBand; n++ ) { // choose 'worst case'-SMR from shifted analysis windows + SMR0.L[n] = maxf ( SMR0.L[n], SMR1.L[n] ); + SMR0.R[n] = maxf ( SMR0.R[n], SMR1.R[n] ); + SMR0.M[n] = maxf ( SMR0.M[n], SMR1.M[n] ); + SMR0.S[n] = maxf ( SMR0.S[n], SMR1.S[n] ); + } + + LEAVE(50); + return SMR0; +} Index: mppenc/branches/r2d/libmpcpsy/psy_tab.c =================================================================== --- mppenc/branches/r2d/libmpcpsy/psy_tab.c (revision 59) +++ mppenc/branches/r2d/libmpcpsy/psy_tab.c (revision 59) @@ -0,0 +1,462 @@ +/* + * Musepack audio compression + * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include "mppenc.h" + +// Antialiasing for calculation of the subband power +const float Butfly [7] = { 0.5f, 0.2776f, 0.1176f, 0.0361f, 0.0075f, 0.000948f, 0.0000598f }; + +// Antialiasing for calculation of the masking thresholds +const float InvButfly [7] = { 2.f, 3.6023f, 8.5034f, 27.701f, 133.33f, 1054.852f, 16722.408f }; + +// w_low for long 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 +const int wl [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 38, 41, 44, 47, 50, 54, 58, 62, 67, 72, 78, 84, 91, 98, 106, 115, 124, 134, 145, 157, 170, 184, 199, 216, 234, 254, 276, 301, 329, 360, 396, 437, 485 }; +const int wh [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 37, 40, 43, 46, 49, 53, 57, 61, 66, 71, 77, 83, 90, 97, 105, 114, 123, 133, 144, 156, 169, 183, 198, 215, 233, 253, 275, 300, 328, 359, 395, 436, 484, 511 }; +// Width: 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 4 5 5 6 6 7 7 8 9 9 10 11 12 13 14 15 17 18 20 22 25 28 31 36 41 48 27 + +// inverse partition-width for long +const float iw [PART_LONG] = { 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/4, 1.f/4, 1.f/4, 1.f/5, 1.f/5, 1.f/6, 1.f/6, 1.f/7, 1.f/7, 1.f/8, 1.f/9, 1.f/9, 1.f/10, 1.f/11, 1.f/12, 1.f/13, 1.f/14, 1.f/15, 1.f/17, 1.f/18, 1.f/20, 1.f/22, 1.f/25, 1.f/28, 1.f/31, 1.f/36, 1.f/41, 1.f/48, 1.f/27 }; + +// w_low for short 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +const int wl_short [PART_SHORT] = { 0, 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 23, 29, 36, 46, 59, 75, 99 }; +const int wh_short [PART_SHORT] = { 0, 1, 2, 3, 5, 6, 7, 9, 12, 14, 18, 23, 29, 36, 46, 58, 75, 99, 127 }; + +// inverse partition-width for short +const float iw_short [PART_SHORT] = { 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/4, 1.f/6, 1.f/7, 1.f/8, 1.f/11, 1.f/13, 1.f/17, 1.f/25, 1.f/29 }; + +/* +Nr. wl wh fl fh bl bh bm Nr. wl wh fl fh bl bh bm + 0: 0 0 0 0 0.000000 0.000000 0.000000 + 1: 1 1 43 43 0.425460 0.425460 0.425460 0: 0 0 0 0 0.000000 0.000000 0.000000 + 2: 2 2 86 86 0.850241 0.850241 0.850241 + + 3: 3 3 129 129 1.273448 1.273448 1.273448 + 4: 4 4 172 172 1.694205 1.694205 1.694205 1: 1 1 172 172 1.694205 1.694205 1.694205 + 5: 5 5 215 215 2.111672 2.111672 2.111672 + + 6: 6 6 258 258 2.525051 2.525051 2.525051 + 7: 7 7 301 301 2.933594 2.933594 2.933594 2: 2 2 345 345 3.336612 3.336612 3.336612 + 8: 8 8 345 345 3.336612 3.336612 3.336612 + + 9: 9 9 388 388 3.733479 3.733479 3.733479 +10: 10 10 431 431 4.123635 4.123635 4.123635 3: 3 3 517 517 4.881924 4.881924 4.881924 +11- 11 12 474 517 4.506591 4.881924 4.695234 + +12: 13 14 560 603 5.249283 5.608381 5.429880 +13: 15 16 646 689 5.958998 6.300971 6.131073 4: 4 5 689 861 6.300971 7.581073 6.958618 +14: 17 18 732 775 6.634195 6.958618 6.797509 + +15: 19 20 818 861 7.274232 7.581073 7.428745 +16: 21 22 904 947 7.879211 8.168753 8.025049 5: 5 6 861 1034 7.581073 8.722594 8.168753 +17: 23 24 991 1034 8.449828 8.722594 8.587239 + +18: 25 26 1077 1120 8.987223 9.243908 9.116546 +19: 27 28 1163 1206 9.492850 9.734263 9.614484 6: 6 7 1034 1206 8.722594 9.734263 9.243908 +20: 29 30 1249 1292 9.968365 10.195382 10.082745 + +21: 31 32 1335 1378 10.415539 10.629064 10.523116 +22: 33 34 1421 1464 10.836184 11.037125 10.937413 7: 8 9 1378 1550 10.629064 11.421352 11.037125 +23: 35 37 1507 1593 11.232108 11.605071 11.421352 + +24: 38 40 1637 1723 11.783474 12.125139 11.956764 +25: 41 43 1766 1852 12.288791 12.602659 12.447904 8: 10 12 1723 2067 12.125139 13.316883 12.753228 +26: 44 46 1895 1981 12.753228 13.042468 12.899777 + +27: 47 49 2024 2110 13.181453 13.448898 13.316883 +28: 50 53 2153 2283 13.577635 13.945465 13.764881 9: 12 14 2067 2412 13.316883 14.288198 13.825796 +29: 54 57 2326 2455 14.062349 14.397371 14.232693 + +30: 58 61 2498 2627 14.504172 14.811258 14.660130 +31: 62 66 2670 2842 14.909464 15.283564 15.100115 10: 15 18 2584 3101 14.711029 15.795819 15.283564 +32: 67 71 2885 3058 15.372757 15.714074 15.546390 + +33: 72 77 3101 3316 15.795819 16.185532 15.994471 +34: 78 83 3359 3575 16.259980 16.616871 16.441494 11: 18 23 3101 3962 15.795819 17.204658 16.547424 +35: 84 90 3618 3876 16.685418 17.079349 16.885941 + +36: 91 97 3919 4177 17.142352 17.506445 17.327264 +37: 98 105 4221 4522 17.564981 17.959646 17.765487 12: 23 29 3962 4996 17.204658 18.533945 17.904788 +38: 106 114 4565 4910 18.014031 18.433233 18.227034 + +39: 115 123 4953 5297 18.483782 18.874805 18.682185 +40: 124 133 5340 5728 18.922095 19.332992 19.130789 13: 29 36 4996 6202 18.533945 19.801451 19.198897 +41: 134 144 5771 6202 19.377073 19.801451 19.592946 + +42: 145 156 6245 6718 19.842285 20.272889 20.061808 +43: 157 169 6761 7278 20.310373 20.739167 20.529583 14: 36 46 6202 7924 19.801451 21.222342 20.565177 +44: 170 183 7321 7881 20.773175 21.191895 20.987911 + +45: 184 198 7924 8527 21.222342 21.623344 21.428652 +46: 199 215 8570 9259 21.650236 22.050787 21.857360 15: 46 58 7924 9991 21.222342 22.420001 21.882271 +47: 216 233 9302 10034 22.074042 22.440072 22.263795 + +48- 234 253 10078 10896 22.459969 22.807140 22.640652 +49: 254 275 10939 11843 22.823891 23.144847 22.991444 16: 59 75 10164 12920 22.499251 23.461146 23.044078 +50: 276 300 11886 12920 23.158772 23.461146 23.317264 + +51: 301 328 12963 14126 23.472530 23.748999 23.617861 +52: 329 359 14169 15461 23.758199 24.005540 23.888450 17: 75 99 12920 17054 23.461146 24.248491 23.920884 +53: 360 395 15504 17011 24.012922 24.242660 24.134368 + +54: 396 436 17054 18777 24.248491 24.454928 24.357873 +55: 437 484 18820 20844 24.459492 24.647977 24.559711 18: 99 127 17054 21878 24.248491 24.727775 24.524955 +56: 485 511 20887 22007 24.651498 24.737100 24.695685 +*/ + + +/* V A R I A B L E S */ +float MinVal [PART_LONG]; // contains minimum tonality soffsets +float Loudness [PART_LONG]; // weighting factors for loudness calculation +float SPRD [PART_LONG] [PART_LONG]; // tabulated spreading function +float O_MAX; +float O_MIN; +float FAC1; +float FAC2; // constants for offset calculation +float partLtq [PART_LONG]; // threshold in quiet (partitions) +float invLtq [PART_LONG]; // inverse threshold in quiet (partitions, long) +float fftLtq [512]; // threshold in quiet (FFT) +float Ltq_offset; // Offset for threshold in quiet +float Ltq_max; // maximum level for threshold in quiet +float TMN; +float NMT; +float TransDetect; +unsigned int EarModelFlag; +int MinValChoice; + + +/* + * Klemm 1994 and 1997. Experimental data. Sorry, data looks a little bit + * dodderly. Data below 30 Hz is extrapolated from other material, above 18 + * kHz the ATH is limited due to the original purpose (too much noise at + * ATH is not good even if it's theoretically inaudible). + */ + +static float +ATHformula_Frank ( float freq ) +{ + /* + * one value per 100 cent = 1 + * semitone = 1/4 + * third = 1/12 + * octave = 1/40 decade + * rest is linear interpolated, values are currently in millibel rel. 20 µPa + */ + static short tab [] = { + /* 10.0 */ 9669, 9669, 9626, 9512, + /* 12.6 */ 9353, 9113, 8882, 8676, + /* 15.8 */ 8469, 8243, 7997, 7748, + /* 20.0 */ 7492, 7239, 7000, 6762, + /* 25.1 */ 6529, 6302, 6084, 5900, + /* 31.6 */ 5717, 5534, 5351, 5167, + /* 39.8 */ 5004, 4812, 4638, 4466, + /* 50.1 */ 4310, 4173, 4050, 3922, + /* 63.1 */ 3723, 3577, 3451, 3281, + /* 79.4 */ 3132, 3036, 2902, 2760, + /* 100.0 */ 2658, 2591, 2441, 2301, + /* 125.9 */ 2212, 2125, 2018, 1900, + /* 158.5 */ 1770, 1682, 1594, 1512, + /* 199.5 */ 1430, 1341, 1260, 1198, + /* 251.2 */ 1136, 1057, 998, 943, + /* 316.2 */ 887, 846, 744, 712, + /* 398.1 */ 693, 668, 637, 606, + /* 501.2 */ 580, 555, 529, 502, + /* 631.0 */ 475, 448, 422, 398, + /* 794.3 */ 375, 351, 327, 322, + /* 1000.0 */ 312, 301, 291, 268, + /* 1258.9 */ 246, 215, 182, 146, + /* 1584.9 */ 107, 61, 13, -35, + /* 1995.3 */ -96, -156, -179, -235, + /* 2511.9 */ -295, -350, -401, -421, + /* 3162.3 */ -446, -499, -532, -535, + /* 3981.1 */ -513, -476, -431, -313, + /* 5011.9 */ -179, 8, 203, 403, + /* 6309.6 */ 580, 736, 881, 1022, + /* 7943.3 */ 1154, 1251, 1348, 1421, + /* 10000.0 */ 1479, 1399, 1285, 1193, + /* 12589.3 */ 1287, 1519, 1914, 2369, +#if 0 + /* 15848.9 */ 3352, 4865, 5942, 6177, + /* 19952.6 */ 6385, 6604, 6833, 7009, + /* 25118.9 */ 7066, 7127, 7191, 7260, +#else + /* 15848.9 */ 3352, 4352, 5352, 6352, + /* 19952.6 */ 7352, 8352, 9352, 9999, + /* 25118.9 */ 9999, 9999, 9999, 9999, +#endif + }; + double freq_log; + unsigned index; + + if ( freq < 10. ) freq = 10.; + if ( freq > 29853. ) freq = 29853.; + + freq_log = 40. * log10 (0.1 * freq); /* 4 steps per third, starting at 10 Hz */ + index = (unsigned) freq_log; + return 0.01 * (tab [index] * (1 + index - freq_log) + tab [index+1] * (freq_log - index)); +} + + +/* F U N C T I O N S */ +// calculation of the threshold in quiet in FFT-resolution +static void +Ruhehoerschwelle ( unsigned int EarModelFlag, + int Ltq_offset, + int Ltq_max ) +{ + int n; + int k; + float f; + float erg; + double tmp; + float absLtq [512]; + + for ( n = 0; n < 512; n++ ) { + f = (float) ( (n+1) * (float)(SampleFreq / 2000.) / 512 ); // Frequency in kHz + + switch ( EarModelFlag / 100 ) { + case 0: // ISO-threshold in quiet + tmp = 3.64*pow (f,-0.8) - 6.5*exp (-0.6*(f-3.3)*(f-3.3)) + 0.001*pow (f, 4.0); + break; + default: + case 1: // measured threshold in quiet (Nick Berglmeir, Andree Buschmann, Kopfhörer) + tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); + break; + case 2: // measured threshold in quiet (Filburt, Kopfhörer) + tmp = 9.00*pow (f,-0.5) - 15.0*exp (-0.1*(f-4.0)*(f-4.0)) + 0.0341796875*pow (f, 2.5) + 15.*exp (-(f-0.1)*(f-0.1)) - 18; + tmp = mind ( tmp, Ltq_max - 18 ); + break; + case 3: + tmp = ATHformula_Frank ( 1.e3 * f ); + break; + case 4: + tmp = ATHformula_Frank ( 1.e3 * f ); + if ( f > 4.8 ) { + tmp += 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); + tmp *= 0.5 ; + } + break; + case 5: + tmp = ATHformula_Frank ( 1.e3 * f ); + if ( f > 4.8 ) { + tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); + } + break; + } + + tmp -= f * f * (int)(EarModelFlag % 100 - 50) * 0.0015; // 00: +30 dB, 100: -30 dB @20 kHz + + tmp = mind ( tmp, Ltq_max ); // Limit ATH + tmp += Ltq_offset - 23; // Add chosen Offset + fftLtq[n] = absLtq[n] = POW10 ( 0.1 * tmp); // conversion into power + } + + // threshold in quiet in partitions (long) + for ( n = 0; n < PART_LONG; n++ ) { + erg = 1.e20f; + for ( k = wl[n]; k <= wh[n]; k++ ) + erg = minf (erg, absLtq[k]); + + partLtq[n] = erg; // threshold in quiet + invLtq [n] = 1.f / partLtq[n]; // Inverse + } +} + +#ifdef _WIN32 +static double +asinh ( double x ) +{ + return x >= 0 ? log (sqrt (x*x+1) + x) : -log (sqrt (x*x+1) - x); +} +#endif + + +static double +Freq2Bark ( double Hz ) // Klemm 2002 +{ + return 9.97074*asinh (1.1268e-3 * Hz) - 6.25817*asinh (0.197193e-3 * Hz) ; +} + +static double +Bark2Freq ( double Bark ) // Klemm 2002 +{ + return 956.86 * sinh (0.101561*Bark) + 11.7296 * sinh (0.304992*Bark) + 6.33622e-3*sinh (0.538621*Bark); +} + +static double +LongPart2Bark ( int Part ) +{ + return Freq2Bark ((wl [Part] + wh [Part]) * SampleFreq / 2048.); +} + +// calculating the table for loudness calculation based on absLtq = ank +static void +Loudness_Tabelle (void) +{ + int n; + float midfreq; + float tmp; + + // ca. dB(A) + for ( n = 0; n < PART_LONG; n++ ){ + midfreq = (wh[n] + wl[n] + 3) * (0.25 * SampleFreq / 512); // center frequency in kHz, why +3 ??? + tmp = LOG10 (midfreq) - 3.5f; // dB(A) + tmp = -10 * tmp * tmp + 3 - midfreq/3000; + Loudness [n] = POW10 ( 0.1 * tmp ); // conversion into power + } +} + + +static double +Bass ( float f, float TMN, float NMT, float bass ) +{ + static unsigned char lfe [11] = { 120, 100, 80, 60, 50, 40, 30, 20, 15, 10, 5 }; + int tmp = (int) ( 1024/44100. * f + 0.5 ); + + switch ( tmp ) { + case 0: + case 1: + case 2: + case 3: + case 4: + case 5: + case 6: + case 7: + case 8: + case 9: + case 10: + return TMN + bass * lfe [tmp]; + case 11: + case 12: + case 13: + case 14: + case 15: + case 16: + case 17: + case 18: + return TMN; + case 19: + case 20: + case 21: + case 22: + return TMN*0.75 + NMT*0.25; + case 23: + case 24: + return TMN*0.50 + NMT*0.50; + case 25: + case 26: + return TMN*0.25 + NMT*0.75; + default: + return NMT; + } +} + + +// calculating the coefficient for utilization of the tonality offset, depending on TMN und NMT +static void +Tonalitaetskoeffizienten ( void ) +{ + double tmp; + int n; + float bass; + + bass = 0.1/8 * NMT; + if ( MinValChoice <= 2 && bass > 0.1 ) + bass = 0.1f; + if ( MinValChoice <= 1 ) + bass = 0.0f; + + // alternative: calculation of the minval-values dependent on TMN and TMN + for ( n = 0; n < PART_LONG; n++ ) { + tmp = Bass ( (wl [n] + wh [n]) / 2048. * SampleFreq, TMN, NMT, bass ); + MinVal [n] = POW10 ( -0.1 * tmp ); // conversion into power + } + + // calculation of the constants for "tonality offset" + O_MAX = POW10 ( -0.1 * TMN ); + O_MIN = POW10 ( -0.1 * NMT ); + FAC1 = POW10 ( -0.1 * (NMT - (TMN - NMT) * 0.229) ) ; + FAC2 = (TMN - NMT) * (0.99011159 * 0.1); +} + + +// calculation of the spreading function +static void +Spread ( void ) +{ + int i; + int j; + float tmpx; + float tmpy; + float tmpz; + float x; + + // calculation of the spreading-function for all occuring values + for ( i = 0; i < PART_LONG; i++ ) { // i is masking Partition, Source + for ( j = 0; j < PART_LONG; j++ ) { // j is masking Partition, Target + tmpx = LongPart2Bark (j) - LongPart2Bark (i);// Difference of the partitions in Bark + tmpy = tmpz = 0.; // tmpz = 0: no dip + + if ( tmpx < 0 ) { // downwards (S1) + tmpy = -32.f * tmpx; // 32 dB per Bark, e33 (10) + } + else if ( tmpx > 0 ) { // upwards (S2) +#if 0 + x = (wl[i]+wh[i])/2 * (float)(SampleFreq / 2000)/512; // center frequency in kHz ??????? + if (i==0) x = 0.5f * (float)(SampleFreq / 2000)/512; // if first spectral line +#else + x = i ? wl[i]+wh[i] : 1; + x *= SampleFreq / 1000. / 2048; // center frequency in kHz +#endif + // dB/Bark + tmpy = (22.f + 0.23f / x) * tmpx; // e33 (10) + + // dip (up to 6 dB) + tmpz = 8 * minf ( (tmpx-0.5f) * (tmpx-0.5f) - 2 * (tmpx-0.5f), 0.f ); + } + + // calculate coefficient + SPRD[i][j] = POW10 ( -0.1 * (tmpy+tmpz) ); // [Source] [Target] + } + } + + // Normierung e33 (10) + for ( i = 0; i < PART_LONG; i++ ) { // i is masked Partition + float norm = 0.f; + for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition + norm += SPRD [j] [i]; + for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition + SPRD [j] [i] /= norm; + } +} + +// call all initialisation procedures +void +Init_Psychoakustiktabellen ( void ) +{ + Max_Band = (int) ( Bandwidth * 64. / SampleFreq ); + if ( Max_Band < 1 ) Max_Band = 1; + if ( Max_Band > 31 ) Max_Band = 31; + + Tonalitaetskoeffizienten (); + Ruhehoerschwelle ( EarModelFlag, Ltq_offset, Ltq_max ); + Loudness_Tabelle (); + Spread (); +} + +/* end of psy_tab.c */ Index: mppenc/branches/r2d/src/ans.c =================================================================== --- mppenc/branches/r2d/src/ans.c (revision 58) +++ (revision ) @@ -1,305 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -/* - * Depending on how transient it is, it can be further reduced (up to 0=No ANS). - * Estimate coefficient for feedback at Order=1 over Mask_fu - Mask_fo. - * 3 quantization routines: Order=0, Order=1, Order=2...6 - * Order doesn't specify the power of the noise shaping, but only the flexibility of the form. - * Don't reset utilization of the "remains" at the frame borders; - * "remains"-utilization as scalefactor-independent values, - * so that a utilization beyond Subframe/Frame Borders is even possible. - */ - -#include "mppenc.h" - - -static float InvFourier [MAX_NS_ORDER + 1] [16]; -static float Cos_Tab [16] [MAX_NS_ORDER + 1]; -static float Sin_Tab [16] [MAX_NS_ORDER + 1]; -unsigned int NS_Order; // Maximum order for ANS -unsigned int NS_Order_L [32]; -unsigned int NS_Order_R [32]; // frame-wise order of the Noiseshaping (0: off, 1...5: on) -float FIR_L [32] [MAX_NS_ORDER]; -float FIR_R [32] [MAX_NS_ORDER]; // contains FIR-Filter for NoiseShaping -float ANSspec_L [MAX_ANS_LINES]; -float ANSspec_R [MAX_ANS_LINES]; // L/R-masking thresholds for ANS -float ANSspec_M [MAX_ANS_LINES]; -float ANSspec_S [MAX_ANS_LINES]; // M/S-masking thresholds for ANS - - -void -Init_ANS ( void ) -{ - int n; - int k; - - // calculate Fourier tables - for ( k = 0; k <= MAX_NS_ORDER; k++ ) { - for ( n = 0; n < 16; n++ ) { - InvFourier [k] [n] = (float) cos ( +2*M_PI/64 * (2*n) * k ) / 16.; - Cos_Tab [n] [k] = (float) cos ( -2*M_PI/64 * (2*n+1) * (k+1) ); - Sin_Tab [n] [k] = (float) sin ( -2*M_PI/64 * (2*n+1) * (k+1) ); - } - } -} - - -// calculates optimal reflection coefficients and time response of a prediction filter in LPC analysis -static __inline void -durbin_akf_to_kh1( float* k, // out: reflection coefficients - float* h, // out: time response - const float* akf ) // in : autocorrelation function (0..1 used) -{ - h[0] = k[0] = akf [1] / akf [0]; -} - -static __inline void -durbin_akf_to_kh2( float* k, // out: reflection coefficients - float* h, // out: time response - const float* akf ) // in : autocorrelation function (0..2 used) -{ - float tk,e; - - tk = akf [1] / akf[0]; - e = akf[0] * (1. - tk*tk); - h[0] = k[0] = tk; - h[0] *= 1. - (h[1] = k[1] = tk = (akf[2] - h[0] * akf[1]) / e); -} - -static __inline void -durbin_akf_to_kh3( float* k, // out: reflection coefficients - float* h, // out: time response - const float* akf ) // in : autocorrelation function (0..3 used) -{ - float a,b,tk,e; - - tk = akf[1] / akf[0]; - e = akf[0] * (1. - tk*tk); - h[0] = k[0] = tk; - - tk = (akf[2] - h[0] * akf[1]) / e; - e *= 1. - tk*tk; - h[0] *= 1. - (h[1] = k[1] = tk); - h[2] = k[2] = tk = (akf[3] - h[0] * akf[2] - h[1] * akf[1]) / e; - - h[0] = (a=h[0]) - (b=h[1])*tk; - h[1] = b - a*tk; -} - - -static __inline void -durbin_akf_to_kh ( float* k, // out: reflection coefficients - float* h, // out: time response - float* akf, // in : autocorrelation function (0..n used) - const int n ) // in : number of parameters to calculate -{ - int i,j; - float s,a,b,tk,e; - float* p; - float* q; - - e = akf [0]; - for ( i = 0; i < n; i++ ) { - s = 0.f; - p = h; - q = akf+i; - j = i; - while ( j-- ) - s += *p++ * *q--; - - tk = (akf[i+1] - s) / e; - e *= 1. - tk*tk; - h[i] = k[i] = tk; - p = h; - q = h + i - 1; - - for ( ; p < q; p++, q-- ) { - a = *p; - b = *q; - *p = a - b*tk; - *q = b - a*tk; - } - if ( p == q ) - *p *= 1. - tk; - } -} - -static const unsigned char maxANSOrder [32] = { - 6, 5, 4, 3, 2, 2, 2, 2, - 2, 2, 2, 2, 1, 1, 1, 1, - 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, -}; - -static void -FindOptimalANS ( const int MaxBand, - const unsigned char* ms, - const float* spec0, - const float* spec1, - unsigned int* NS, - float* snr_comp, - float fir [] [MAX_NS_ORDER], - const float* smr0, - const float* smr1, - const int scf [] [3], - const int Transient [32] ) -{ - int Band; - int n; - int k; - int order; - float akf [MAX_NS_ORDER + 1]; - float h [MAX_NS_ORDER]; - float reflex [MAX_NS_ORDER]; - float spec [16]; - float invspec [16]; - float norm; - float ns_loss; - float min_spec; - float min_diff; - float re; - float im; - float ns_energy; - float gain; - float NS_Gain; - float actSMR; - int max; - const float* tmp; - - ENTER(235); - for ( Band = 0; Band <= MaxBand && maxANSOrder[Band]; Band++ ) { - - if ( scf[Band][0] != scf[Band][1] || scf[Band][1] != scf[Band][2] ) - continue; - - if ( Transient[Band] ) - continue; - - max = maxANSOrder [Band]; - - if ( ms[Band] ) { // setting pointer and SMR in relation to the M/S-flag - tmp = &spec1 [Band<<4]; // pointer to MS-data - actSMR = smr1 [Band]; // selecting SMR - } - else { - tmp = &spec0 [Band<<4]; // pointer to LR-data - actSMR = smr0 [Band]; // selecting SMR - } - - if ( actSMR >= 1. ) { - NS_Gain = 1.f; // reset gain - norm = 1.e-30f; - - // Selection of the masking threshold of the current subband, also considering frequency inversion in every 2nd subband - if ( Band & 1 ) - for ( n = 0, tmp += 15; n < 16; n++ ) - norm += spec[n] = *tmp--; - else - for ( n = 0; n < 16; n++ ) - norm += spec[n] = *tmp++; - - // Preprocessing: normalization of the the power of spec[] to 1, and search for minimum of masking threshold - norm = 16.f / norm; - min_spec = 1.e+12f; - for ( n = 0; n < 16; n++ ) { - invspec[n] = 1.f / (spec[n] *= norm); - if ( spec[n] < min_spec ) // normalize spec[] - min_spec = spec[n]; - } - - // Calculation of the auto-correlation function - tmp = InvFourier [0]; - for ( k = 0; k <= max; k++, tmp += 16 ) { - akf[k] = tmp[ 0]*invspec[ 0] + tmp[ 1]*invspec[ 1] + tmp[ 2]*invspec[ 2] + tmp[ 3]*invspec[ 3] + - tmp[ 4]*invspec[ 4] + tmp[ 5]*invspec[ 5] + tmp[ 6]*invspec[ 6] + tmp[ 7]*invspec[ 7] + - tmp[ 8]*invspec[ 8] + tmp[ 9]*invspec[ 9] + tmp[10]*invspec[10] + tmp[11]*invspec[11] + - tmp[12]*invspec[12] + tmp[13]*invspec[13] + tmp[14]*invspec[14] + tmp[15]*invspec[15]; - } - - // Searching for the noise-shaper with maximum gain - for ( order = 1; order <= max; order++ ) { - switch ( order ) { // calculating best FIR-Filter for the return - case 1: durbin_akf_to_kh1 (reflex, h, akf); break; - case 2: durbin_akf_to_kh2 (reflex, h, akf); break; - case 3: durbin_akf_to_kh3 (reflex, h, akf); break; - default: durbin_akf_to_kh (reflex, h, akf, order); break; - } - - ns_loss = 1.e-30f; // estimating the gain - min_diff = 1.e+12f; - for ( n = 0; n < 16; n++ ) { - re = 1.f; // calculating the obtained noise shaping - im = 0.f; - for ( k = 0; k < order; k++ ) { - re -= h[k] * Cos_Tab[n][k]; - im += h[k] * Sin_Tab[n][k]; - } - - ns_energy = re*re + im*im; // calculated spectral shaped noise - ns_loss += ns_energy; // noise energy increases with shaping - - if ( spec[n] < min_diff * ns_energy ) // Searching for minimum distance between the shaped noise and the masking threshold - min_diff = spec[n] / ns_energy; - } - - // Updating the Filter if new gain is bigger than old gain and if the extra noise power through shaping is smaller than the SMR of this band - gain = 16. * min_diff / (min_spec * ns_loss); - if ( gain > NS_Gain && ns_loss < actSMR ) { - NS [Band] = order; - NS_Gain = gain; - memcpy ( fir [Band], h, order * sizeof(*h) ); - } - } - - if ( NS_Gain > 1.f ) { // Activation of ANS if there is gain - snr_comp[Band] *= NS_Gain; - } - } - } - - LEAVE(235); - return; -} - - -// perform ANS-analysis (calculation of FIR-filter and gain) -void -NS_Analyse ( const int MaxBand, - const unsigned char* MSflag, - const SMRTyp smr, - const int* Transient ) -{ - ENTER(10); - - // for L or M, respectively - memset ( FIR_L, 0, sizeof FIR_L ); // reset FIR - memset ( NS_Order_L, 0, sizeof NS_Order_L ); // reset Flags - FindOptimalANS ( MaxBand, MSflag, ANSspec_L, ANSspec_M, NS_Order_L, SNR_comp_L, FIR_L, smr.L, smr.M, SCF_Index_L, Transient ); - - // for R or S, respectively - memset ( FIR_R, 0, sizeof FIR_R ); // reset FIR - memset ( NS_Order_R, 0, sizeof NS_Order_R ); // reset Flags - FindOptimalANS ( MaxBand, MSflag, ANSspec_R, ANSspec_S, NS_Order_R, SNR_comp_R, FIR_R, smr.R, smr.S, SCF_Index_R, Transient ); - - LEAVE(10); - return; -} - -/* end of ans.c */ Index: mppenc/branches/r2d/src/cvd.c =================================================================== --- mppenc/branches/r2d/src/cvd.c (revision 58) +++ (revision ) @@ -1,267 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#include "mppenc.h" - -/* C O N S T A N T S */ -// from MatLab-Simulation (Fourier-transforms of the Cos-Rolloff) -#if 0 -static const float Puls [11] = { - -0.02724753942504f, -0.10670808991329f, -0.06198987803623f, 0.18006206051664f, - 0.49549552704050f, 0.64201253447071f, 0.49549552704050f, 0.18006206051664f, - -0.06198987803623f, -0.10670808991329f, -0.02724753942504f -}; -#endif - -static const float Puls [ 9] = { - -0.10670808991329f, -0.06198987803623f, 0.18006206051664f, 0.49549552704050f, - 0.64201253447071f, 0.49549552704050f, 0.18006206051664f, -0.06198987803623f, - -0.10670808991329f -}; - -/* -// Generating the Cos-Rolloff of the Cepstral-analysis, Cos-Rolloff from 5512,5 Hz to 11025 Hz -// for ( k = 0; k <= 1024; k++ ) { -// if (k < 256) CosWin [k-256] = 1; -// else if (k < 512) CosWin [k-256] = 0.5 + 0.5*cos (M_PI*(k-256)/256); -// else CosWin [k-256] = 0; -// } -*/ -static const float CosWin [256] = { - 1.0000000000000000f, 0.9999623298645020f, 0.9998494386672974f, 0.9996612071990967f, 0.9993977546691895f, 0.9990590810775757f, 0.9986452460289002f, 0.9981563091278076f, 0.9975923895835877f, 0.9969534873962402f, 0.9962397813796997f, 0.9954513311386108f, 0.9945882558822632f, 0.9936507344245911f, 0.9926388263702393f, 0.9915527701377869f, 0.9903926253318787f, 0.9891586899757385f, 0.9878510832786560f, 0.9864699840545654f, 0.9850156307220459f, 0.9834882616996765f, 0.9818880558013916f, 0.9802152514457703f, 0.9784701466560364f, 0.9766530394554138f, 0.9747641086578369f, 0.9728036522865295f, 0.9707720279693604f, 0.9686695337295532f, 0.9664964079856873f, 0.9642530679702759f, 0.9619397521018982f, 0.9595569372177124f, 0.9571048617362976f, 0.9545840024948120f, 0.9519946575164795f, 0.9493372440338135f, 0.9466121792793274f, 0.9438198208808899f, 0.9409606456756592f, 0.9380350708961487f, 0.9350435137748718f, 0.9319864511489868f, 0.9288643002510071f, 0.9256775975227356f, 0.9224267601966858f, 0.9191123247146606f, 0.9157348275184631f, 0.9122946262359619f, 0.9087924361228943f, 0.9052286148071289f, 0.9016037583351135f, 0.8979184627532959f, 0.8941732048988342f, 0.8903686404228210f, 0.8865052461624146f, 0.8825836181640625f, 0.8786044120788574f, 0.8745682239532471f, 0.8704755902290344f, 0.8663271069526672f, 0.8621235489845276f, 0.8578653931617737f, - 0.8535534143447876f, 0.8491881489753723f, 0.8447702527046204f, 0.8403005003929138f, 0.8357794880867004f, 0.8312078714370728f, 0.8265864253044128f, 0.8219157457351685f, 0.8171966671943665f, 0.8124297261238098f, 0.8076158165931702f, 0.8027555346488953f, 0.7978496551513672f, 0.7928989529609680f, 0.7879040837287903f, 0.7828658819198608f, 0.7777851223945618f, 0.7726625204086304f, 0.7674987912178040f, 0.7622948288917542f, 0.7570513486862183f, 0.7517691850662231f, 0.7464491128921509f, 0.7410919070243835f, 0.7356983423233032f, 0.7302693724632263f, 0.7248056530952454f, 0.7193081378936768f, 0.7137775421142578f, 0.7082147598266602f, 0.7026206851005554f, 0.6969960331916809f, 0.6913416981697083f, 0.6856585741043091f, 0.6799474954605103f, 0.6742093563079834f, 0.6684449315071106f, 0.6626551747322083f, 0.6568408608436585f, 0.6510030031204224f, 0.6451423168182373f, 0.6392598152160645f, 0.6333563923835754f, 0.6274328231811523f, 0.6214900612831116f, 0.6155290603637695f, 0.6095505952835083f, 0.6035556793212891f, 0.5975451469421387f, 0.5915199518203735f, 0.5854809284210205f, 0.5794290900230408f, 0.5733652114868164f, 0.5672903656959534f, 0.5612053275108337f, 0.5551111102104187f, 0.5490085482597351f, 0.5428986549377441f, 0.5367822647094727f, 0.5306603908538818f, 0.5245338082313538f, 0.5184035897254944f, 0.5122706294059753f, 0.5061357617378235f, - 0.5000000000000000f, 0.4938642382621765f, 0.4877294003963471f, 0.4815963804721832f, 0.4754661619663239f, 0.4693396389484406f, 0.4632177054882050f, 0.4571013450622559f, 0.4509914219379425f, 0.4448888897895813f, 0.4387946724891663f, 0.4327096343040466f, 0.4266347587108612f, 0.4205709397792816f, 0.4145190417766571f, 0.4084800481796265f, 0.4024548530578613f, 0.3964443206787109f, 0.3904493749141693f, 0.3844709396362305f, 0.3785099089145660f, 0.3725671768188477f, 0.3666436076164246f, 0.3607401549816132f, 0.3548576533794403f, 0.3489970266819000f, 0.3431591391563416f, 0.3373448550701141f, 0.3315550684928894f, 0.3257906734943390f, 0.3200524747371674f, 0.3143413960933685f, 0.3086582720279694f, 0.3030039668083191f, 0.2973793447017670f, 0.2917852103710175f, 0.2862224578857422f, 0.2806918919086456f, 0.2751943469047546f, 0.2697306573390961f, 0.2643016278743744f, 0.2589081227779388f, 0.2535509169101715f, 0.2482308149337769f, 0.2429486215114594f, 0.2377051562070847f, 0.2325011938810349f, 0.2273375093936920f, 0.2222148776054382f, 0.2171340882778168f, 0.2120959013700485f, 0.2071010768413544f, 0.2021503448486328f, 0.1972444802522659f, 0.1923841983079910f, 0.1875702589750290f, 0.1828033626079559f, 0.1780842244625092f, 0.1734135746955872f, 0.1687921136617661f, 0.1642205268144608f, 0.1596994996070862f, 0.1552297323942184f, 0.1508118808269501f, - 0.1464466154575348f, 0.1421345919370651f, 0.1378764659166336f, 0.1336728632450104f, 0.1295244395732880f, 0.1254318058490753f, 0.1213955804705620f, 0.1174163669347763f, 0.1134947761893272f, 0.1096313893795013f, 0.1058267876505852f, 0.1020815446972847f, 0.0983962342143059f, 0.0947714000940323f, 0.0912075936794281f, 0.0877053514122963f, 0.0842651948332787f, 0.0808876454830170f, 0.0775732174515724f, 0.0743224024772644f, 0.0711356922984123f, 0.0680135712027550f, 0.0649565011262894f, 0.0619649514555931f, 0.0590393692255020f, 0.0561801902949810f, 0.0533878505229950f, 0.0506627671420574f, 0.0480053536593914f, 0.0454160086810589f, 0.0428951233625412f, 0.0404430739581585f, 0.0380602329969406f, 0.0357469581067562f, 0.0335035994648933f, 0.0313304923474789f, 0.0292279683053494f, 0.0271963365375996f, 0.0252359099686146f, 0.0233469791710377f, 0.0215298328548670f, 0.0197847411036491f, 0.0181119665503502f, 0.0165117643773556f, 0.0149843730032444f, 0.0135300243273377f, 0.0121489353477955f, 0.0108413146808743f, 0.0096073597669601f, 0.0084472559392452f, 0.0073611787520349f, 0.0063492907211185f, 0.0054117450490594f, 0.0045486823655665f, 0.0037602325901389f, 0.0030465149320662f, 0.0024076367262751f, 0.0018436938989908f, 0.0013547716662288f, 0.0009409435442649f, 0.0006022718735039f, 0.0003388077020645f, 0.0001505906548118f, 0.0000376490788767f, -}; - - -/* F U N C T I O N S */ -// sets all the harmonics -static void -SetVoiceLines ( int* VoiceLine, const float base, int val ) -{ - int n; - int max = (int) (MAX_CVD_LINE * base / 1024.f); // harmonics up to Index MAX_CVD_LINE (spectral lines outside of that don't make sense) - int line; - float frq = 1024.f / base; // frq = 1024./i is the Index of the basic harmonic - - // go through all harmonics - for ( n = 1; n <= max; n++ ) { - line = (int) (n * frq); - VoiceLine [line] = VoiceLine [line+1] = val; - } -} - - -// Analyze the Cepstrum, search for the basic harmonic -static void -CEP_Analyse2048 ( float* res1, - float* res2, - float* qual1, - float* qual2, - float* cep ) -{ - int n; - int line; - float cc [MAX_ANALYZED_IDX + 3]; // cross correlation - float ref; - float line_sum; - float sum; - float kkf; - float norm; - const float* x; - - // cross-correlation with pulse shape - // Calculate idx = MIN_ANALYZED_IDX-2 to MAX_ANALYZED_IDX+2, - // because they are read during search for maximum - // 50 -> 882 Hz, 700 -> 63 Hz base frequency - - *res1 = *res2 = 0. ; - memset ( cc, 0, sizeof cc ); - - for ( n = MIN_ANALYZED_IDX - 2; n <= MAX_ANALYZED_IDX + 2; n++ ) { - x = cep + n; - if ( x[0] > 0 ) { - norm = x[-4] * x[-4] + - x[-3] * x[-3] + - x[-2] * x[-2] + - x[-1] * x[-1] + - x[ 0] * x[ 0] + - x[ 1] * x[ 1] + - x[ 2] * x[ 2] + - x[ 3] * x[ 3] + - x[ 4] * x[ 4]; - kkf = x[-4] * Puls [0] + - x[-3] * Puls [1] + - x[-2] * Puls [2] + - x[-1] * Puls [3] + - x[ 0] * Puls [4] + - x[ 1] * Puls [5] + - x[ 2] * Puls [6] + - x[ 3] * Puls [7] + - x[ 4] * Puls [8]; - cc [n] = kkf * kkf / norm; // calculate the square of ncc to avoid sqrt() - } - } - - // search for the (relative) maximum - ref = 0.f; - line = MED_ANALYZED_IDX; - for ( n = MAX_ANALYZED_IDX; n >= MED_ANALYZED_IDX; n-- ) { - if ( - cc[n] * cep[n] * cep[n] > ref && - cc[n] > 0.40f && // e33 (02) 0.85 - cep[n] > 0.00f && // e33 (02) - cc[n ] >= cc[n+1] && - cc[n ] >= cc[n-1] && - cc[n+1] >= cc[n+2] && - cc[n-1] >= cc[n-2] - ) - { - ref = cc[n] * cep[n] * cep[n]; - line = n; - } - } - - // Calculating the center of the maximum (Interpolation) - x = cep + line; - sum = x[-3] + x[-2] + x[-1] + x[0] + x[1] + x[2] + x[3] + 1.e-30f; - line_sum = (x[1]-x[-1]) + 2 * (x[2]-x[-2]) + 3 * (x[3]-x[-3]) + sum * line + 1.e-30f; - - /* e33 (04) */ - ref = cc[line ] * cep[line ] * cep[line ] - + cc[line-1] * cep[line-1] * cep[line-1] - + cc[line+1] * cep[line+1] * cep[line+1]; - - //{ - // static unsigned int x = 0; - // - // printf ("%7.3f s ", (x/2)*1152./44100 ); - // x++; - //} - - //printf ("ref=%5.3f *res1=%7.3f f=%8.3f ", ref, line_sum / sum, 44100. / (line_sum / sum) ); - - *qual1 = ref; - if ( ref > 0.015f ) - *res1 = line_sum / sum; - - if ( CVD_used < 2 ) - return; - - // search for the (relative) maximum - ref = 0.f; - line = MIN_ANALYZED_IDX; - - for ( n = MED_ANALYZED_IDX + 1; n >= MIN_ANALYZED_IDX - 1; n-- ) { - cc [2*n ] += 0.5 * cc [n]; - cc [2*n+1] += 0.5 * (cc [n] + cc[n+1]); - cep [2*n ] += 0.5 * cep [n]; - cep [2*n+1] += 0.5 * (cep [n] + cep[n+1]); - } - - for ( n = 2*MED_ANALYZED_IDX; n >= 2*MIN_ANALYZED_IDX; n-- ) { - if ( - cc[n] * cep[n] * cep[n] > ref && - cc[n] > 0.85f && /* e33 (02) */ - cep[n] > 0.00f && /* e33 (02) */ - cc[n ] >= cc[n+1] && - cc[n ] >= cc[n-1] && - cc[n+1] >= cc[n+2] && - cc[n-1] >= cc[n-2] - ) - { - ref = cc[n] * cep[n] * cep[n]; - line = n; - } - } - - // Calculating the center of the maximum (Interpolation) - x = cep + line; - sum = x[-3] + x[-2] + x[-1] + x[0] + x[1] + x[2] + x[3] + 1.e-30f; - line_sum = (x[1]-x[-1]) + 2 * (x[2]-x[-2]) + 3 * (x[3]-x[-3]) + sum * line + 1.e-30f; - - /* e33 (04) */ - ref = cc[line ] * cep[line ] * cep[line ] - + cc[line-1] * cep[line-1] * cep[line-1] - + cc[line+1] * cep[line+1] * cep[line+1]; - - //printf ("ref=%5.3f *res2=%8.3f f=%8.3f\n", ref, 0.5 * line_sum / sum, 44100. / (0.5 * line_sum / sum) ); - - *qual2 = ref; - if ( ref >= 0.1f ) - *res2 = 0.5 * line_sum / sum; - - return; -} - -#ifndef CVD_FASTLOG -# define logfast(x) ((float) log (x)) -#else - -static __inline float /* This is a rough estimation with an accuracy of |x|<0.0037 */ -logfast ( float x ) -{ - double y = x * x; - y *= y; - y *= y; - return (((int*)(&y))[1] + (45127.5 - 1072693248.)) * ( M_LN2 / (1L<<23) ); -} - -#endif - -// ClearVoiceDetection for spectrum *spec -// input : Spectrum *spec -// output: Array *vocal contains information if the FFT-Line is a harmonic component -int -CVD2048 ( const float* spec, int* vocal ) -{ - static float cep [4096]; // cep[4096] -- array, which is also used for the 2048 FFT - const float* win = CosWin; // pointer to cos-roll-off - float res1; - float res2; - float qual1; - float qual2; - int n; - - ENTER(20); - // Calculating logarithmated, windowed spectrum cep[] - // cep[512...1024] = 0 -- cep[1025...2047] doesn't matter, because the first have to be filled by fft - for ( n = 0; n < 256; n++ ) - cep[n] = logfast (*spec++); - for ( n = 256; n < 512; n++ ) - cep[n] = logfast (*spec++) * *win++; - - memset ( cep+512, 0, 513*sizeof(*cep) ); - - // Calculating cepstrum of cep[] (the function Cepstrum() outputs the cepstrum in-place) - Cepstrum2048 ( cep, MAX_ANALYZED_IDX ); - - // search the harmonic - CEP_Analyse2048 ( &res1, &res2, &qual1, &qual2, cep ); -//#include "cvd.h" - if ( res1 > 0.f || res2 > 0.f ) { - if ( res1 > 0. ) SetVoiceLines ( vocal, res1, 100 ); - if ( res2 > 0. ) SetVoiceLines ( vocal, res2, 20 ); - LEAVE(20); - return 1; - } - LEAVE(20); - return 0; -} Index: mppenc/branches/r2d/src/cvd.h =================================================================== --- mppenc/branches/r2d/src/cvd.h (revision 58) +++ (revision ) @@ -1,9 +1,0 @@ -{ - static FILE* fp = NULL; - static int x = 0; - - if ( fp == NULL ) fp = fopen ( "cvd.txt", "a" ); - fprintf ( fp, "%7.3f %6.2f %7.2f\n", (x>>1)*1152./44100, res1, res2 ); - - x++; -} Index: mppenc/branches/r2d/src/fastmath.c =================================================================== --- mppenc/branches/r2d/src/fastmath.c (revision 58) +++ (revision ) @@ -1,85 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#include "mppenc.h" - -#ifdef FAST_MATH - -const float tabatan2 [ 2*TABSTEP+1] [2]; -const float tabcos [26*TABSTEP+1] [2]; -const float tabsqrt_ex [256]; -const float tabsqrt_m [ TABSTEP+1] [2]; - - -void Init_FastMath ( void ) -{ - int i; - float X; - float Y; - double xm; - double x0; - double xp; - double x; - double y; - float* p; - - p = (float*) tabatan2; - for ( i = -TABSTEP; i <= TABSTEP; i++ ) { - xm = atan ((i-0.5)/TABSTEP); - x0 = atan ((i+0.0)/TABSTEP); - xp = atan ((i+0.5)/TABSTEP); - x = x0/2 + (xm + xp)/4; - y = xp - xm; - *p++ = x; - *p++ = y; - } - - p = (float*) tabcos; - for ( i = -13*TABSTEP; i <= 13*TABSTEP; i++ ) { - xm = cos ((i-0.5)/TABSTEP); - x0 = cos ((i+0.0)/TABSTEP); - xp = cos ((i+0.5)/TABSTEP); - x = x0/2 + (xm + xp)/4; - y = xp - xm; - *p++ = x; - *p++ = y; - } - - p = (float*) tabsqrt_ex; - for ( i = 0; i < 255; i++ ) { - *(int*)&X = (i << 23); - *(int*)&Y = (i << 23) + (1<<23) - 1; - *p++ = sqrt(X); - } - *(int*)&X = (255 << 23) - 1; - *p++ = sqrt(X); - - p = (float*) tabsqrt_m; - for ( i = 1*TABSTEP; i <= 2*TABSTEP; i++ ) { - xm = sqrt ((i-0.5)/TABSTEP); - x0 = sqrt ((i+0.0)/TABSTEP); - xp = sqrt ((i+0.5)/TABSTEP); - x = x0/2 + (xm + xp)/4; - y = xp - xm; - *p++ = x; - *p++ = y; - } -} - -#endif Index: mppenc/branches/r2d/src/fastmath.h =================================================================== --- mppenc/branches/r2d/src/fastmath.h (revision 58) +++ (revision ) @@ -1,94 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#if 1 -# define ROUND32(x) ( floattmp = (x) + (int)0x00FD8000L, *(int*)(&floattmp) - (int)0x4B7D8000L ) -#else -# define ROUND32(x) ( (int) floor ((x) + 0.5) ) -#endif - -#ifdef FAST_MATH - -static __inline float -my_atan2 ( float x, float y ) -{ - float t; - int i; - float ret; - float floattmp; - - if ( (*(int*)&x & 0x7FFFFFFF) < (*(int*)&y & 0x7FFFFFFF) ) { - i = ROUND32 (t = TABSTEP * (x / y)); - ret = tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (t-i); - if ( *(int*)&y < 0 ) - ret = (float)(ret - M_PI); - } - else if ( *(int*)&x < 0) { - i = ROUND32 (t = TABSTEP * (y / x)); - 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)); - ret = + M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t); - } - else { - ret = 0.; - } - return ret; -} - - -static __inline float -my_cos ( float x ) -{ - float t; - int i; - float ret; - float floattmp; - - i = ROUND32 (t = TABSTEP * x); - ret = tabcos [13*TABSTEP+i][0] + tabcos [13*TABSTEP+i][1] * (t-i); - return ret; -} - - -static __inline int -my_ifloor ( float x ) -{ - x = x + (0x0C00000L + 0.500000001); - return *(int*)&x - 1262485505; -} - - -static __inline float -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 - ret *= tabsqrt_ex [ex]; - return ret; -} - -#endif Index: mppenc/branches/r2d/src/fft4g.c =================================================================== --- mppenc/branches/r2d/src/fft4g.c (revision 58) +++ (revision ) @@ -1,670 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#include "mppenc.h" - -/* F U N C T I O N S */ -static void makewt ( const int nw, int* ip, float* w ); -static void makect ( const int nc, int* ip, float* c ); -static __inline void bitrv2 ( const int n, int* ip, float* a ); // -static __inline void cftfsub ( const int n, float* a, float* w ); // -static __inline void rftfsub ( const int n, float* a, int nc, float* c ); // -static __inline void cft1st ( const int n, float* a, float* w ); // -static __inline void cftmdl_i386 ( const int n, const int l, float* a, float* w ); // 5648 -static __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 - -// generates lookup-tables -void -Generate_FFT_Tables ( const int n, int* ip, float* w ) -{ - int nw; - int nc; - - nw = n >> 2; - makewt ( nw, ip, w ); - - nc = n >> 2; - makect ( nc, ip, w + nw ); -} - - -// patched to only-forward -void -rdft ( const int n, float* a, int* ip, float* w ) -{ - float xi; - - ENTER(30); - if ( n > 4) { - bitrv2 ( n, ip + 2, a ); - cftfsub ( n, a, w ); - rftfsub ( n, a, ip[1], w + ip[0] ); - } - else if ( n == 4 ) { - cftfsub ( n, a, w ); - } - xi = a[0] - a[1]; - a[0] += a[1]; - a[1] = xi; - LEAVE(30); - return; -} - - -/* -------- initializing routines -------- */ -static void -makewt ( const int nw, int* ip, float* w ) -{ - int j; - int nwh; - float x; - float y; - double delta; - - ENTER(31); - ip[0] = nw; - ip[1] = 1; - if ( nw > 2 ) { - nwh = nw >> 1; - delta = (M_PI/4) / nwh; - w[0] = 1.; - w[1] = 0.; - w[nwh] = COS (delta * nwh); - w[nwh + 1] = w[nwh]; - if ( nwh > 2 ) { - for ( j = 2; j < nwh; j += 2 ) { - x = COS (delta * j); - y = SIN (delta * j); - w[j] = x; - w[j + 1] = y; - w[nw - j] = y; - w[nw - j + 1] = x; - } - bitrv2 ( nw, ip + 2, w ); - } - } - LEAVE(31); - return; -} - - -static void -makect ( const int nc, int* ip, float* c ) -{ - int j; - int nch; - double delta; - - ENTER(32); - ip[1] = nc; - if ( nc > 1 ) { - nch = nc >> 1; - delta = (M_PI/4) / nch; - c[0] = COS (delta * nch); - c[nch] = 0.5f * c[0]; - for ( j = 1; j < nch; j++ ) { - c[j] = 0.5f * COS (delta * j); - c[nc - j] = 0.5f * SIN (delta * j); - } - } - LEAVE(32); - return; -} - - -/* -------- 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; - - ENTER(33); - 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; - } - } - } - LEAVE(33); - 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; - - ENTER(34); - 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 ); - } - LEAVE(34); - 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; - - ENTER(35); - 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 ); - LEAVE(35); - return; -} - -extern void Cdecl cftmdl_3DNow_1 ( const int n, const int l, float* a, float* w ); -extern void 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; - - ENTER(36); - 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); - } - LEAVE(36); - - ENTER(39); - 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 ); - } - LEAVE(39); - 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; - - ENTER(36); - cftmdl_3DNow_1 (n,l,a,w); - LEAVE(36); - - ENTER(39); - 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 ); - } - LEAVE(39); - 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; - - ENTER(37); - 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 ); - LEAVE(37); - return; -} - -/* end of fft4g.c */ Index: mppenc/branches/r2d/src/fft_routines.c =================================================================== --- mppenc/branches/r2d/src/fft_routines.c (revision 58) +++ (revision ) @@ -1,337 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#include "mppenc.h" - -#define CX0 -1. -#define CX1 0.5 - -#define SX1 -1. -#define SX2 (2./9/ 1) -#define SX3 (2./9/ 4) -#define SX4 (2./9/ 10) -#define SX5 (2./9/ 20) -#define SX6 (2./9/ 35) -#define SX7 (2./9/ 56) -#define SX8 (2./9/ 84) -#define SX9 (2./9/120) -#define SX10 (2./9/165) - - -#ifdef EXTRA_DECONV -# define DECONV \ - { \ - tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \ - aix[ 0] -= CX0*tmp; \ - aix[ 2] -= CX1*tmp; \ - tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \ - aix[ 3] -= SX1*tmp; \ - aix[ 5] -= SX2*tmp; \ - aix[ 7] -= SX3*tmp; \ - aix[ 9] -= SX4*tmp; \ - aix[11] -= SX5*tmp; \ - } -#elif 0 -# define DECONV \ - { \ - float A[20]; \ - int i; \ - memcpy (A, aix, 20*sizeof(aix)); \ - tmp = (CX0*aix[0] + CX1*aix[2]) * (1./(CX0*CX0+CX1*CX1)); \ - aix[ 0] -= CX0*tmp; \ - aix[ 2] -= CX1*tmp; \ - tmp = (SX1*aix[3] + SX2*aix[5] + SX3*aix[7] + SX4*aix[9] + SX5*aix[11]) * (1./(SX1*SX1+SX2*SX2+SX3*SX3+SX4*SX4+SX5*SX5)); \ - aix[ 3] -= SX1*tmp; \ - aix[ 5] -= SX2*tmp; \ - aix[ 7] -= SX3*tmp; \ - aix[ 9] -= SX4*tmp; \ - aix[11] -= SX5*tmp; \ - for ( i=0; i<10; i++) \ - printf ("%u%9.0f%7.0f%9.0f%7.0f\n",i, A[i+i], A[i+i+1], aix[i+i], aix[i+i+1] ); \ - } -#else -# define DECONV -#endif - - -/* V A R I A B L E S */ -static int ip [4096]; // bitinverse for maximum 2048 FFT -static float w [4096]; // butterfly-coefficient for maximum 2048 FFT -static float a [4096]; // holds real input for FFT -static float Hann_256 [ 256]; -static float Hann_1024 [1024]; -static float Hann_1600 [1600]; - - -////////////////////////////// -// -// 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 ( float* window, unsigned int size, float alpha ) -{ - double sumvalue = 0.; - double scale; - int i; - - scale = 0.25 / sqrt (size); - for ( i = 0; i < (int)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 < (int)size/2; i++ ) - window [size-1-i] = window [i] = /*sqrt*/ ( window [i] / sumvalue ) * scale; -} - -static void -CosWindow ( float* window, unsigned int size ) -{ - double x; - double scale; - int i; - - scale = 0.25 / sqrt (size); - for ( i = 0; i < (int)size/2; i++ ) { - x = cos ( (i+0.5) * (M_PI / size) ); - window [size/2-1-i] = window [size/2+i] = scale * x * x; - } -} - -static void -Window ( float* window, unsigned int size, float alpha ) -{ - if ( alpha < 0. ) - CosWindow ( window, size ) ; - else - KBDWindow ( window, size, alpha ); -} - - -/* F U N C T I O N S */ -// generates FFT lookup-tables -void -Init_FFT ( void ) -{ - int n; - double x; - double scale; - - // normalized hann functions - Window ( Hann_256 , 256, KBD1 ); - Window ( Hann_1024, 1024, KBD2 ); - scale = 0.25 / sqrt (2048.); - for ( n = 0; n < 800; n++ ) - x = cos ((n+0.5) * (M_PI/1600)), Hann_1600 [799-n] = Hann_1600 [800+n] = (float)(x * x * scale); - - Generate_FFT_Tables ( 2048, ip, w ); -} - -// input : Signal *x -// output: energy spectrum *erg -void -PowSpec256 ( const float* x, float* erg ) -{ - const float* win = Hann_256; - float* aix = a; - int i; - - ENTER(40); - // windowing - i = 256; - while (i--) - *aix++ = *x++ * *win++; - - // perform FFT - rdft ( 256, a, ip, w ); - - // calculate power - aix = a; // reset pointer - i = 128; - while (i--) { - *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; - aix += 2; - } - LEAVE(40); -} - -// input : Signal *x -// output: energy spectrum *erg -void -PowSpec1024 ( const float* x, float* erg ) -{ - const float* win = Hann_1024; - float* aix = a; - int i; - - ENTER(41); - i = 1024; // windowing - while (i--) - *aix++ = *x++ * *win++; - -// for (i=0; i<1024; i++) -// a[i] = Hann_1024[i] * ((i==0 ? 0 : i-512) + 1000); - - rdft ( 1024, a, ip, w ); // perform FFT - - aix = a; // calculate power - i = 512; - - - DECONV; -// for (i = 0; i <= 512; i++ ) -// printf ("%3u %12.6f %12.6f\n", i, a[i+i], a[i+i+1]); -// exit(1); - while (i--) { - *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; - aix += 2; - } - LEAVE(41); -} - -// input : Signal *x -// output: energy spectrum *erg -void -PowSpec2048 ( const float* x, float* erg ) -{ - const float* win = Hann_1600; - float* aix = a; - int i; - - ENTER(42); - // windowing (only 1600 samples available -> centered in 2048!) - memset ( a , 0, 224*sizeof(*a) ); - aix = a + 224; - i = 1600; - while (i--) - *aix++ = *x++ * *win++; - memset ( a+1824, 0, 224*sizeof(*a) ); - - rdft ( 2048, a, ip, w ); // perform FFT - - aix = a; // calculate power - i = 1024; - while (i--) { - *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; - aix += 2; - } - LEAVE(42); -} - -#include "fastmath.h" - -// input : Signal *x -// output: energy spectrum *erg and phase spectrum *phs -void -PolarSpec1024 ( const float* x, float* erg, float* phs ) -{ - const float* win = Hann_1024; - float* aix = a; - int i; - - ENTER(43); - i = 1024; // windowing - while (i--) - *aix++ = *x++ * *win++; - - rdft ( 1024, a, ip, w ); // perform FFT - - // calculate power and phase - aix = a; // reset pointer - i = 512; - while (i--) { - *erg++ = aix[0]*aix[0] + aix[1]*aix[1]; - *phs++ = ATAN2F (aix[1], aix[0]); - aix += 2; - } - LEAVE(43); -} - -// input : logarithmized energy spectrum *cep -// output: Cepstrum *cep (in-place) -void -Cepstrum2048 ( float* cep, const int MaxLine ) -{ - float* aix = cep; - float* bix = cep + 2048; - int i; - - ENTER(44); - // generate real, even spectrum (symmetric around 1024, cep[2048-i] = cep[i]) - for ( i = 0; i < 1024; i++ ) - *bix-- = *aix++; - - // perform IFFT - rdft ( 2048, cep, ip, w ); - - // only real part as outcome (all even indexes of cep[]) - aix = cep; - bix = cep; - i = MaxLine + 1; - while (i--) { - *aix = *bix * (float) (0.9888 / 2048.); -// *aix = *bix * 0.0004828125f; - aix ++; - bix += 2; - } - LEAVE(44); -} Index: mppenc/branches/r2d/src/psy.c =================================================================== --- mppenc/branches/r2d/src/psy.c (revision 58) +++ (revision ) @@ -1,1309 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -/* - * Prediction - * Short-Block-detection with smooth inset - * revise CalcMSThreshold - * /dev/audio for Windows too - * revise PNS/IS - * CVS with smoother inset - * several files per call - * revise ANS with changing SCFs - - * No IS - * PNS estimation very rough, also IS should be used to reduce data rate in the side channel - * ANS problems at Frame boundaries when resolution changes - * ANS problems at Subframe boundaries when SCF changes - * CVS+ with smoother transition - ----------------------------------------- - -Optimize Tabelle[18] (use second table) -CVS+ - -- ANS is disregarded during the search for the best Res -- ANS messes up if res changes (each 36 samples) and/or SCF changes (each 12 samples) -- PNS not in difference signal - -- implement IS in decoder -- Experimental Quantizer with complete energy preservation - - 1D, calculated - - 2D, calculated - - 2D, manually modified, coeffs set to 1.f - - */ - -#include "mppenc.h" - -/* V A R I A B L E S */ -/* further switches for the psymodel */ -unsigned int CVD_used; // global flag for ClearVoiceDetection -float varLtq; // variable threshold in quiet -unsigned int tmpMask_used; // global flag for temporal masking -float ShortThr; // Factor to calculate the masking threshold with transients -float minSMR; // minimum SMR for all subbands - -float a [PART_LONG]; -float b [PART_LONG]; -float c [PART_LONG]; -float d [PART_LONG]; // Integrations for tmpMask -static float Xsave_L [3 * 512]; -static float Xsave_R [3 * 512]; // FFT-Amplitudes L/R -static float Ysave_L [3 * 512]; -static float Ysave_R [3 * 512]; // FFT-Phases L/R -float T_L [PART_LONG]; -float T_R [PART_LONG]; // time-constants for tmpMask -float pre_erg_L[2][PART_SHORT]; -float pre_erg_R[2][PART_SHORT]; // Preecho-control short -float PreThr_L [PART_LONG]; -float PreThr_R [PART_LONG]; // for Pre-Echo-control L/R -float tmp_Mask_L [PART_LONG]; -float tmp_Mask_R [PART_LONG]; // for Post-Masking L/R -int Vocal_L [MAX_CVD_LINE + 4]; -int Vocal_R [MAX_CVD_LINE + 4]; // FFT-Line belongs to harmonic? - -/* F U N C T I O N S */ -// Resets Arrays -void -Init_Psychoakustik ( void ) -{ - int i; - - ENTER(200); - // generate FFT lookup-tables with largest FFT-size of 1024 - Init_FFT (); - - // setting pre-echo variables to Ltq - for ( i = 0; i < PART_LONG; i++ ) { - pre_erg_L [0][i/3] = pre_erg_R [0][i/3] = - pre_erg_L [1][i/3] = pre_erg_R [1][i/3] = - tmp_Mask_L [i] = tmp_Mask_R [i] = - PreThr_L [i] = PreThr_R [i] = partLtq [i]; - } - - // initializing arrays with zero - memset ( Xsave_L, 0, sizeof Xsave_L ); - memset ( Xsave_R, 0, sizeof Xsave_R ); - memset ( Ysave_L, 0, sizeof Ysave_L ); - memset ( Ysave_R, 0, sizeof Ysave_R ); - memset ( a, 0, sizeof a ); - memset ( b, 0, sizeof b ); - memset ( c, 0, sizeof c ); - memset ( d, 0, sizeof d ); - memset ( T_L, 0, sizeof T_L ); - memset ( T_R, 0, sizeof T_R ); - memset ( Vocal_L, 0, sizeof Vocal_L ); - memset ( Vocal_R, 0, sizeof Vocal_R ); - - LEAVE(200); - return; -} - - -// VBRmode 1: Adjustment of all SMRs via a factor (offset of SMRoffset dB) -// VBRmode 2: SMRs have a minimum of minSMR dB -static void -RaiseSMR_Signal ( const int MaxBand, float* signal, float tmp ) -{ - int Band; - float z = 0.; - - for ( Band = MaxBand; Band >= 0; Band-- ) { - if ( z < signal [Band] ) z = signal [Band]; - if ( z > tmp ) z = tmp; - if ( signal [Band] < z ) signal [Band] = z; - } -} - - -void -RaiseSMR ( const int MaxBand, SMRTyp* smr ) -{ - float tmp = POW10 ( 0.1 * minSMR ); - - ENTER(201); - RaiseSMR_Signal ( MaxBand, smr->L, tmp ); - RaiseSMR_Signal ( MaxBand, smr->R, tmp ); - RaiseSMR_Signal ( MaxBand, smr->M, tmp ); - RaiseSMR_Signal ( MaxBand, smr->S, 0.5 * tmp ); - - LEAVE(201); - return; -} - -// input : *smr -// output: *smr, *ms, *x (only the entries for L/R contain relevant data) -// Check if either M/S- or L/R-coding has a lower perceptual entropy -// Choose the better mode, copy the appropriate data into the -// arrays that belong to L and R and set the ms-Flag accordingly. -void -MS_LR_Entscheidung ( const int MaxBand, unsigned char* ms, SMRTyp* smr, SubbandFloatTyp* x ) -{ - int Band; - int n; - float PE_MS; - float PE_LR; - float tmpM; - float tmpS; - float* l; - float* r; - - ENTER(202); - - for ( Band = 0; Band <= MaxBand; Band++ ) { // calculate perceptual entropy - PE_LR = PE_MS = 1.f; - if (smr->L[Band] > 1.) PE_LR *= smr->L[Band]; - if (smr->R[Band] > 1.) PE_LR *= smr->R[Band]; - if (smr->M[Band] > 1.) PE_MS *= smr->M[Band]; - if (smr->S[Band] > 1.) PE_MS *= smr->S[Band]; - - if ( PE_MS < PE_LR ) { - ms[Band] = 1; - - // calculate M/S-signal and copies it to L/R-array - l = x[Band].L; - r = x[Band].R; - for ( n = 0; n < 36; n++, l++, r++ ) { - tmpM = (*l + *r) * 0.5f; - tmpS = (*l - *r) * 0.5f; - *l = tmpM; - *r = tmpS; - } - - // copy M/S - SMR to L/R-fields - smr->L[Band] = smr->M[Band]; - smr->R[Band] = smr->S[Band]; - } - else { - ms[Band] = 0; - } - } - - LEAVE(202); - return; -} - -// input : FFT-spectrums *spec0 und *spec1 -// output: energy in the individual subbands *erg0 and *erg1 -// With Butfly[], you can calculate the results of aliasing during calculation -// of subband energy from the FFT-spectrums. -static void -SubbandEnergy ( const int MaxBand, - float* erg0, - float* erg1, - const float* spec0, - const float* spec1 ) -{ - int n; - int k; - int alias; - float tmp0; - float tmp1; - - ENTER(203); - - // Is this here correct for FFT-based data or is this calculation rule only for MDCTs??? - - for ( k = 0; k <= MaxBand; k++ ) { // subband index - tmp0 = tmp1 = 0.f; - for ( n = 0; n < 16; n++, spec0++, spec1++ ) { // spectral index - tmp0 += *spec0; - tmp1 += *spec1; - - // Consideration of Aliasing between the subbands - if ( n < +sizeof(Butfly)/sizeof(*Butfly) && k != 0 ) { - alias = -1 - (n<<1); - tmp0 += Butfly [n] * (spec0[alias] - *spec0); - tmp1 += Butfly [n] * (spec1[alias] - *spec1); - } - else if ( n > 15-sizeof(Butfly)/sizeof(*Butfly) && k != 31 ) { - alias = 31 - (n<<1); - tmp0 += Butfly [15-n] * (spec0[alias] - *spec0); - tmp1 += Butfly [15-n] * (spec1[alias] - *spec1); - } - } - *erg0++ = tmp0; - *erg1++ = tmp1; - } - - LEAVE(203); - return; -} - -// input : FFT-Spectrums *spec0 and *spec1 -// output: energy in the individual partitions *erg0 and *erg1 -static void -PartitionEnergy ( float* erg0, - float* erg1, - const float* spec0, - const float* spec1 ) -{ - unsigned int n; - unsigned int k; - float e0; - float e1; - - ENTER(204); - -#if 000000 - for ( n = 0; n < PART_LONG; n++ ) { - k = wh[n] - wl[n]; - e0 = *spec0++; - e1 = *spec1++; - while ( k-- ) { - e0 += *spec0++; - e1 += *spec1++; - } - *erg0++ = e0; - *erg1++ = e1; - } -#else - n = 0; - - for ( ; n < 23; n++ ) { // 11 or 23 - k = wh[n] - wl[n]; - e0 = *spec0++; - e1 = *spec1++; - while ( k-- ) { - e0 += *spec0++; - e1 += *spec1++; - } - *erg0++ = e0; - *erg1++ = e1; - } - - for ( ; n < 48; n++ ) { // 37 ... 46, 48, 57 - k = wh[n] - wl[n]; - e0 = sqrt (*spec0++); - e1 = sqrt (*spec1++); - while ( k-- ) { - e0 += sqrt (*spec0++); - e1 += sqrt (*spec1++); - } - *erg0++ = e0*e0 * iw[n]; - *erg1++ = e1*e1 * iw[n]; - } - - for ( ; n < PART_LONG; n++ ) { - k = wh[n] - wl[n]; - e0 = *spec0++; - e1 = *spec1++; - while ( k-- ) { - e0 += *spec0++; - e1 += *spec1++; - } - *erg0++ = e0; - *erg1++ = e1; - } - - -#endif - - LEAVE(204); - return; -} - - -// input : FFT-Spectrums *spec0, *spec1 and unpredictability *cw0 and *cw1 -// output: weighted energy in the individual partitions *erg0, *erg1 -static void -WeightedPartitionEnergy ( float* erg0, - float* erg1, - const float* spec0, - const float* spec1, - const float* cw0, - const float* cw1 ) -{ - unsigned int n; - unsigned int k; - float e0; - float e1; - - ENTER(205); - -#if 000000 - for ( n = 0; n < PART_LONG; n++ ) { - e0 = *spec0++ * *cw0++; - e1 = *spec1++ * *cw1++; - k = wh[n] - wl[n]; - while ( k-- ) { - e0 += *spec0++ * *cw0++; - e1 += *spec1++ * *cw1++; - } - *erg0++ = e0; - *erg1++ = e1; - } -#else - n = 0; - - for ( ; n < 23; n++ ) { - e0 = *spec0++ * *cw0++; - e1 = *spec1++ * *cw1++; - k = wh[n] - wl[n]; - while ( k-- ) { - e0 += *spec0++ * *cw0++; - e1 += *spec1++ * *cw1++; - } - *erg0++ = e0; - *erg1++ = e1; - } - - for ( ; n < 48; n++ ) { - e0 = sqrt (*spec0++ * *cw0++); - e1 = sqrt (*spec1++ * *cw1++); - k = wh[n] - wl[n]; - while ( k-- ) { - e0 += sqrt (*spec0++ * *cw0++); - e1 += sqrt (*spec1++ * *cw1++); - } - *erg0++ = e0*e0 * iw[n]; - *erg1++ = e1*e1 * iw[n]; - } - - for ( ; n < PART_LONG; n++ ) { - e0 = *spec0++ * *cw0++; - e1 = *spec1++ * *cw1++; - k = wh[n] - wl[n]; - while ( k-- ) { - e0 += *spec0++ * *cw0++; - e1 += *spec1++ * *cw1++; - } - *erg0++ = e0; - *erg1++ = e1; - } -#endif - - LEAVE(205); - return; -} - -// input : masking thresholds, first half of the arrays *shaped0 and *shaped1 -// output: masking thresholds, second half of the arrays *shaped0 and *shaped1 -// Considering the result of aliasing via InvButfly[] -// The input *thr0, *thr1 is gathered via address calculation from *shaped0, *shaped1 - -static void -AdaptThresholds ( const int MaxLine, float* shaped0, float* shaped1 ) -{ - int n; - int mod; - int alias; - float tmp; - const float* invb = InvButfly; - const float* thr0 = shaped0 - 512; - const float* thr1 = shaped1 - 512; - float tmp0; - float tmp1; - - ENTER(206); - - // should be able to optimize it with coasting. [ 9 ] + n * [ 7 + 7 + 2 ] + [ 7 ] - // Schleife Schl Schl Ausr Schleife - for ( n = 0; n < MaxLine; n++, thr0++, thr1++ ) { - mod = n & 15; // n%16 - tmp0 = *thr0; - tmp1 = *thr1; - - if ( mod < +sizeof(InvButfly)/sizeof(*InvButfly) && n > 12 ) { - alias = -1 - (mod<<1); - tmp = thr0[alias] * invb[mod]; - if ( tmp < tmp0 ) tmp0 = tmp; - tmp = thr1[alias] * invb[mod]; - if ( tmp < tmp1 ) tmp1 = tmp; - } - else if ( mod > 15-sizeof(InvButfly)/sizeof(*InvButfly) && n < 499 ) { - alias = 31 - (mod<<1); - tmp = thr0[alias] * invb[15-mod]; - if ( tmp < tmp0 ) tmp0 = tmp; - tmp = thr1[alias] * invb[15-mod]; - if ( tmp < tmp1 ) tmp1 = tmp; - } - *shaped0++ = tmp0; - *shaped1++ = tmp1; - } - - LEAVE(206); - return; -} - -#include "fastmath.h" - -// input : current spectrum in the form of power *spec and phase *phase, -// the last two earlier spectrums are at position -// 512 and 1024 of the corresponding Input-Arrays. -// Array *vocal, which can mark an FFT_Linie as harmonic -// output: current amplitude *amp and unpredictability *cw -static void -CalcUnpred ( const int MaxLine, - const float* spec, - const float* phase, - const int* vocal, - float* amp0, - float* phs0, - float* cw ) -{ - int n; - float amp; - float tmp; -#define amp1 ((amp0) + 512) // amp[ 512...1023] contains data of frame-1 -#define amp2 ((amp0) + 1024) // amp[1024...1535] contains data of frame-2 -#define phs1 ((phs0) + 512) // phs[ 512...1023] contains data of frame-1 -#define phs2 ((phs0) + 1024) // phs[1024...1535] contains data of frame-2 - - ENTER(207); - - for ( n = 0; n < MaxLine; n++ ) { - tmp = COSF ((phs0[n] = phase[n]) - 2*phs1[n] + phs2[n]); // copy phase to output-array, predict phase and calculate predictive error - amp0[n] = SQRTF (spec[n]); // calculate and set amplitude - amp = 2*amp1[n] - amp2[n]; // predict amplitude - - // calculate unpredictability - cw[n] = SQRTF (spec[n] + amp * (amp - 2*amp0[n] * tmp)) / (amp0[n] + FABS(amp)); - } - - // postprocessing of harmonic FFT-lines (*cw is set to CVD_UNPRED) - if ( CVD_used && vocal != NULL ) { - for ( n = 0; n < MAX_CVD_LINE; n++, cw++, vocal++ ) - if ( *vocal != 0 && *cw > CVD_UNPRED * 0.01 * *vocal ) - *cw = CVD_UNPRED * 0.01 * *vocal; - } - - LEAVE(207); - return; -} -#undef amp1 -#undef amp2 -#undef phs1 -#undef phs2 - - -// input : Energy *erg, calibrated energy *werg -// output: spread energy *res, spread weighted energy *wres -// SPRD describes the spreading function as calculated in psy_tab.c -static void -SpreadingSignal ( const float* erg, const float* werg, float* res, float* wres ) -{ - int n; - int k; - int start; - int stop; - const float* sprd; - float e; - float ew; - - ENTER(208); - - for (k=0; k -static float -AdaptLtq ( const float* erg0, const float* erg1 ) -{ - static float loud = 0.f; - float* weight = Loudness; - float sum = 0.f; - int n; - - // calculate loudness - for ( n = 0; n < PART_LONG; n++ ) - sum += (*erg0++ + *erg1++) * *weight++; - - // Utilization of the time constants (fast drop of Ltq T=5, slow rise of Ltq T=20) - //loud = (sum < loud) ? (4 * sum + loud)*0.2f : (19 * loud + sum)*0.05f; - loud = 0.98 * loud + 0.02 * (0.5 * sum); - - // calculate dynamic offset for threshold in quiet, 0...+20 dB, at 96 dB loudness, an offset of 20 dB is assumed - return 1.f + varLtq * loud * 5.023772e-08f; -} - -// input : simultaneous masking threshold *frqthr, -// previous masking threshold *tmpthr, -// Integrations *a (short-time) and *b (long-time) -// output: tracked Integrations *a and *b, time constant *tau -static void -CalcTemporalThreshold ( float* a, float* b, float* tau, float* frqthr, float* tmpthr ) -{ - int n; - float tmp; - - ENTER(220); - - for ( n = 0; n < PART_LONG; n++ ) { - // following calculations relative to threshold in quiet - frqthr[n] *= invLtq[n]; - tmpthr[n] *= invLtq[n]; - - // new post-masking 'tmp' via time constant tau, if old post-masking > Ltq (=1) - tmp = tmpthr[n] > 1.f ? POW ( tmpthr[n], tau[n] ) : 1.f; - - // calculate time constant for post-masking in next frame, - // if new time constant has to be calculated (new tmpMask < frqMask) - a[n] += 0.5f * (frqthr[n] - a[n]); // short time integrator - b[n] += 0.15f * (frqthr[n] - b[n]); // long time integrator - if (tmp < frqthr[n]) - tau[n] = a[n] <= b[n] ? 0.8f : 0.2f + b[n] / a[n] * 0.6f; - - // use post-masking of (Re-Normalization) - tmpthr[n] = maxf (frqthr[n], tmp) * partLtq[n]; - } - - LEAVE(220); - return; -} - -// input : L/R-Masking thresholds in Partitions *thrL, *thrR -// L/R-Subband energies *ergL, *ergR -// M/S-Subband energies *ergM, *ergS -// output: M/S-Masking thresholds in Partitions *thrM, *thrS -static void -CalcMSThreshold ( const float* const ergL, - const float* const ergR, - const float* const ergM, - const float* const ergS, - float* const thrL, - float* const thrR, - float* const thrM, - float* const thrS ) -{ - int n; - float norm; - float tmp; - - // All hardcoded numbers here should be pulled from somewhere, - // the "4.", the -2 dB, the 0.0625 and the 0.9375, as well as all bands where this is done - - for ( n = 0; n < PART_LONG; n++ ) { - // estimate M/S thresholds out of L/R thresholds and M/S and L/R energies - thrS[n] = thrM[n] = maxf (ergM[n], ergS[n]) / maxf (ergL[n], ergR[n]) * minf (thrL[n], thrR[n]); - - switch ( MS_Channelmode ) { // preserve 'near-mid' signal components - case 3: - if ( n > 0 ) { - double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n]; - double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n]; - if ( ratioMS < ratioLR ) { // MS - if ( ergM[n] > ergS[n] ) - thrS[n] = thrL[n] = thrR[n] = 1.e18f; - else - thrM[n] = thrL[n] = thrR[n] = 1.e18f; - } - else { // LR - if ( ergL[n] > ergR[n] ) - thrR[n] = thrM[n] = thrS[n] = 1.e18f; - else - thrL[n] = thrM[n] = thrS[n] = 1.e18f; - } - } - break; - case 4: - if ( n > 0 ) { - double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n]; - double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n]; - if ( ratioMS < ratioLR ) { // MS - if ( ergM[n] > ergS[n] ) - thrS[n] = 1.e18f; - else - thrM[n] = 1.e18f; - } - else { // LR - if ( ergL[n] > ergR[n] ) - thrR[n] = 1.e18f; - else - thrL[n] = 1.e18f; - } - } - break; - case 5: - thrS[n] *= 2.; // +3 dB - break; - case 6: - break; - default: - fprintf ( stderr, "Unknown stereo mode\n"); - case 10: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.70794578f * iw[n]; // -1.5 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = MS2SPAT1 * thrS[n] + (1.f-MS2SPAT1) * tmp; // raises masking threshold by up to 3 dB - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = MS2SPAT1 * thrM[n] + (1.f-MS2SPAT1) * tmp; - } - } - break; - case 11: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.63095734f * iw[n]; // -2.0 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = MS2SPAT2 * thrS[n] + (1.f-MS2SPAT2) * tmp; // raises masking threshold by up to 6 dB - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = MS2SPAT2 * thrM[n] + (1.f-MS2SPAT2) * tmp; - } - } - break; - case 12: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = MS2SPAT3 * thrS[n] + (1.f-MS2SPAT3) * tmp; // raises masking threshold by up to 9 dB - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = MS2SPAT3 * thrM[n] + (1.f-MS2SPAT3) * tmp; - } - } - break; - case 13: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = MS2SPAT4 * thrS[n] + (1.f-MS2SPAT4) * tmp; // raises masking threshold by up to 12 dB - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = MS2SPAT4 * thrM[n] + (1.f-MS2SPAT4) * tmp; - } - } - break; - case 15: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = tmp; // raises masking threshold by up to +oo dB an - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = tmp; - } - } - break; - case 22: - if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB - norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth - if ( ergM[n] > ergS[n] ) { - tmp = ergS[n] * norm; - if ( thrS[n] > tmp ) - thrS[n] = maxf (tmp, ergM[n]*iw[n]*0.025); // +/- 1.414° - } else if ( ergS[n] > ergM[n] ) { - tmp = ergM[n] * norm; - if ( thrM[n] > tmp ) - thrM[n] = maxf (tmp, ergS[n]*iw[n]*0.025); // +/- 1.414° - } - } - break; - } - } - - return; -} - -// input : Masking thresholds in Partitions *partThr0, *partThr1 -// level of threshold in quiet *ltq in FFT-resolution -// output: Masking thresholds in FFT-resolution *thr0, *thr1 -// inline, because it's called 4x -static void -ApplyLtq ( float* thr0, - float* thr1, - const float* partThr0, - const float* partThr1, - const float AdaptedLTQ, - int MSflag ) -{ - int n; - int k; - float ltq; - float tmp; - float ms = MSflag ? 0.125f * AdaptedLTQ : 0.25f * AdaptedLTQ ; - - for ( n = 0; n < PART_LONG; n++ ) { - for ( k = wl[n]; k <= wh[n]; k++, thr0++, thr1++ ) { // threshold in quiet (Partition) -#if 0 - ltq = AdaptedLTQ * fftLtq [k]; - *thr0 = maxf ( partThr0 [n], ltq ); - *thr1 = maxf ( partThr1 [n], ltq ); -#else - // Applies a much more gentle ATH rolloff + 6 dB more dynamic - ltq = sqrt (ms * fftLtq [k]); - tmp = sqrt (partThr0 [n]) + ltq; - *thr0 = tmp * tmp; - tmp = sqrt (partThr1 [n]) + ltq; - *thr1 = tmp * tmp; -#endif - } - } - return; -} - -// input : Subband energies *erg0, *erg1 -// Masking thresholds in FFT-resolution *thr0, *thr1 -// output: SMR per Subband *smr0, *smr1 -static void -CalculateSMR ( const int MaxBand, - const float* erg0, - const float* erg1, - const float* thr0, - const float* thr1, - float* smr0, - float* smr1 ) -{ - int n; - int k; - float tmp0; - float tmp1; - - // calculation of the masked thresholds in the subbands - for (n = 0; n <= MaxBand; n++ ) { - tmp0 = *thr0++; - tmp1 = *thr1++; - for (k=1; k<16; ++k, ++thr0, ++thr1) { - if (*thr0 < tmp0) tmp0 = *thr0; - if (*thr1 < tmp1) tmp1 = *thr1; - } - *smr0++ = 0.0625f * *erg0++ / tmp0; - *smr1++ = 0.0625f * *erg1++ / tmp1; - } - - return; -} - -// input : energy spectrums erg[4][128] (4 delayed FFTs) -// Energy of the last short block *preerg in short partitions -// PreechoFac declares allowed traved of the masking threshold -// output: masking threshold *thr in short partitions -// Energy of the last short block *preerg in short partitions -#if 0 -static void -CalcShortThreshold ( const float erg [] [128], - const float PreechoFac, - float* thr, - float preerg[2][PART_SHORT], - int* transient ) -{ - const int* lo = wl_short; // lower FFT-index - const int* hi = wh_short; // upper FFT-index - const float* iwidth = iw_short; // inverse partition-width - int k; - int n; - int m; - float tmp; - float enrg; - float th; - const float* ep; - - for ( k = 0; k < PART_SHORT; k++, lo++, hi++ ) { - transient[k] = 0; - th = 1.e20f; - for ( n = 0; n < 4; n++ ) { - ep = erg[n] + *lo; - m = *hi - *lo; - enrg = *ep++; - while (m--) - enrg += *ep++; - - // preecho prevention - tmp = enrg; - if (preerg[0][k] < enrg) - enrg = preerg[0][k]; - preerg[0][k] = tmp; - - // is signal transient? - if (tmp > TransDetect*enrg) transient[k] = 1; - - // assume short threshold = engr*PreechoFac - th = minf (th, enrg*PreechoFac); - } - thr[k] = th * *iwidth++; - } - - return; -} -#else -static void -CalcShortThreshold ( const float erg [4] [128], - const float ShortThr, - float* thr, - float old_erg [2][PART_SHORT], - int* transient ) -{ - const int* index_lo = wl_short; // lower FFT-index - const int* index_hi = wh_short; // upper FFT-index - const float* iwidth = iw_short; // inverse partition-width - int k; - int n; - int m; - float new_erg; - float th; - const float* ep; - - for ( k = 0; k < PART_SHORT; k++ ) { - transient [k] = 0; - th = old_erg [0][k]; - for ( n = 0; n < 4; n++ ) { - ep = erg[n] + index_lo [k]; - m = index_hi [k] - index_lo [k]; - - new_erg = *ep++; - while (m--) - new_erg += *ep++; // e = Short_Partition-energy in piece n - - if ( new_erg > old_erg [0][k] ) { // bigger than the old? - - if ( new_erg > old_erg [0][k] * TransDetect || - new_erg > old_erg [1][k] * TransDetect*2 ) // is signal transient? - transient [k] = 1; - } - else { - th = minf ( th, new_erg ); // assume short threshold = engr*PreechoFac - } - - old_erg [1][k] = old_erg [0][k]; - old_erg [0][k] = new_erg; // save the current one - } - thr [k] = th * ShortThr * *iwidth++; // pull out and multiply only when transient[k]=1 - } - - return; -} - -#endif - -// input : previous simultaneous masking threshold *preThr, -// current simultaneous masking threshold *simThr -// output: update of *preThr for next call, -// current masking threshold *partThr -static void -PreechoControl ( float* partThr0, - float* preThr0, - const float* simThr0, - float* partThr1, - float* preThr1, - const float* simThr1 ) -{ - int n; - - for ( n = 0; n < PART_LONG; n++ ) { - *partThr0++ = minf ( *simThr0, *preThr0 * PREFAC_LONG); - *partThr1++ = minf ( *simThr1, *preThr1 * PREFAC_LONG); - *preThr0++ = *simThr0++; - *preThr1++ = *simThr1++; - } - return; -} - - -void -TransientenCalc ( int* T, - const int* TL, - const int* TR ) -{ - int i; - int x1; - int x2; - - memset ( T, 0, 32*sizeof(*T) ); - - for ( i = 0; i < PART_SHORT; i++ ) - if ( TL[i] || TR[i] ) { - x1 = wl_short[i] >> 2; - x2 = wh_short[i] >> 2; - while ( x1 <= x2 ) - T [x1++] = 1; - } -} - - -// input : PCM-Data *data -// output: SMRs for the input data -SMRTyp -Psychoakustisches_Modell ( const int MaxBand, const PCMDataTyp* data, int* TransientL, int* TransientR ) -{ - float Xi_L[32], Xi_R[32]; // acoustic pressure per Subband L/R - float Xi_M[32], Xi_S[32]; // acoustic pressure per Subband M/S - float cw_L[512], cw_R[512]; // unpredictability (only L/R) - float erg0[512], erg1[512]; // holds energy spectrum of long FFT - float phs0[512], phs1[512]; // holds phase spectrum of long FFT - float Thr_L[2*512], Thr_R[2*512]; // masking thresholds L/R, second half for triangle swap - float Thr_M[2*512], Thr_S[2*512]; // masking thresholds M/S, second half for triangle swap - float F_256[4][128]; // holds energies of short FFTs (L/R only) - float Xerg[1024]; // holds energy spectrum of very long FFT - float Ls_L[PART_LONG], Ls_R[PART_LONG]; // acoustic pressure in Partition L/R - float Ls_M[PART_LONG], Ls_S[PART_LONG]; // acoustic pressure per each partition M/S - float PartThr_L[PART_LONG], PartThr_R[PART_LONG]; // masking thresholds L/R (Partition) - float PartThr_M[PART_LONG], PartThr_S[PART_LONG]; // masking thresholds M/S (Partition) - float sim_Mask_L[PART_LONG], sim_Mask_R[PART_LONG]; // simultaneous masking (only L/R) - float clow_L[PART_LONG], clow_R[PART_LONG]; // spread, weighted energy (only L/R) - float cLs_L[PART_LONG], cLs_R[PART_LONG]; // weighted partition energy (only L/R) - float shortThr_L[PART_SHORT],shortThr_R[PART_SHORT]; // threshold for short FFT (only L/R) - int n; - int MaxLine = (MaxBand+1)*16; // set FFT-resolution according to MaxBand - SMRTyp SMR0; - SMRTyp SMR1; // holds SMR's for first and second Analysis - int isvoc_L; - int isvoc_R; - float factorLTQ = 1.f; // Offset after variable LTQ - - ENTER(50); - // 'ClearVocalDetection'-Process - if ( CVD_used ) { - memset ( Vocal_L, 0, sizeof Vocal_L ); - memset ( Vocal_R, 0, sizeof Vocal_R ); - - // left channel - PowSpec2048 ( &data->L[0], Xerg ); - isvoc_L = CVD2048 ( Xerg, Vocal_L ); - // right channel - PowSpec2048 ( &data->R[0], Xerg ); - isvoc_R = CVD2048 ( Xerg, Vocal_R ); - } - - // calculation of the spectral energy via FFT - PolarSpec1024 ( &data->L[0], erg0, phs0 ); // left - PolarSpec1024 ( &data->R[0], erg1, phs1 ); // right - - // calculation of the acoustic pressures per each subband for L/R-signals - SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 ); - - // calculation of the acoustic pressures per each partition - PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 ); - - // calculate the predictability of the signal - // left - memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) ); - memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) ); - CalcUnpred ( MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L ); - // right - memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) ); - memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) ); - CalcUnpred ( MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R ); - - // calculation of the weighted acoustic pressures per each partition - WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R ); - - // Spreading Signal & weighted unpredictability-signal - // left - memset ( clow_L , 0, sizeof clow_L ); - memset ( sim_Mask_L, 0, sizeof sim_Mask_L ); - SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L ); - // right - memset ( clow_R , 0, sizeof clow_R ); - memset ( sim_Mask_R, 0, sizeof sim_Mask_R ); - SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R ); - - // Offset depending on tonality - ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R ); - - // handling of transient signals - // calculate four short FFTs (left) - PowSpec256 ( &data->L[ 0+SHORTFFT_OFFSET], F_256[0] ); - PowSpec256 ( &data->L[144+SHORTFFT_OFFSET], F_256[1] ); - PowSpec256 ( &data->L[288+SHORTFFT_OFFSET], F_256[2] ); - PowSpec256 ( &data->L[432+SHORTFFT_OFFSET], F_256[3] ); - // calculate short Threshold - CalcShortThreshold ( F_256, ShortThr, shortThr_L, pre_erg_L, TransientL ); - - // calculate four short FFTs (right) - PowSpec256 ( &data->R[ 0+SHORTFFT_OFFSET], F_256[0] ); - PowSpec256 ( &data->R[144+SHORTFFT_OFFSET], F_256[1] ); - PowSpec256 ( &data->R[288+SHORTFFT_OFFSET], F_256[2] ); - PowSpec256 ( &data->R[432+SHORTFFT_OFFSET], F_256[3] ); - // calculate short Threshold - CalcShortThreshold ( F_256, ShortThr, shortThr_R, pre_erg_R, TransientR ); - - // dynamic adjustment of the threshold in quiet to the loudness of the current sequence - if ( varLtq > 0. ) - factorLTQ = AdaptLtq ( Ls_L, Ls_R ); - - // utilization of the temporal post-masking - if ( tmpMask_used ) { - CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L ); - CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R ); - memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L ); - memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R ); - } - - // transient signal? - for ( n = 0; n < PART_SHORT; n++ ) { - if ( TransientL [n] ) { - sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] ); - sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] ); - sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] ); - } - if ( TransientR[n] ) { - sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] ); - sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] ); - sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] ); - } - } - - // Pre-Echo control - PreechoControl ( PartThr_L, PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R ); - - // utilization of the threshold in quiet - ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 ); - - // Consideration of aliasing between the subbands (noise is smeared) - // In: Thr[0..511], Out: Thr[512...1023] - AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 ); - memmove ( Thr_L, Thr_L+512, 512*sizeof(float) ); - memmove ( Thr_R, Thr_R+512, 512*sizeof(float) ); - - // calculation of the Signal-to-Mask-Ratio - CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR0.L, SMR0.R ); - - /***************************************************************************************/ - /***************************************************************************************/ - if ( MS_Channelmode > 0 ) { - // calculation of the spectral energy via FFT - PowSpec1024 ( &data->M[0], erg0 ); // mid - PowSpec1024 ( &data->S[0], erg1 ); // side - - // calculation of the acoustic pressures per each subband for M/S-signals - SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 ); - - // calculation of the acoustic pressures per each partition - PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 ); - - // calculate masking thresholds for M/S - CalcMSThreshold ( Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S ); - ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 ); - - // Consideration of aliasing between the subbands (noise is smeared) - // In: Thr[0..511], Out: Thr[512...1023] - AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 ); - memmove ( Thr_M, Thr_M+512, 512*sizeof(float) ); - memmove ( Thr_S, Thr_S+512, 512*sizeof(float) ); - - // calculation of the Signal-to-Mask-Ratio - CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR0.M, SMR0.S ); - } - - if ( NS_Order > 0 ) { // providing the Noise Shaping thresholds - memcpy ( ANSspec_L, Thr_L, sizeof ANSspec_L ); - memcpy ( ANSspec_R, Thr_R, sizeof ANSspec_R ); - memcpy ( ANSspec_M, Thr_M, sizeof ANSspec_M ); - memcpy ( ANSspec_S, Thr_S, sizeof ANSspec_S ); - } - /***************************************************************************************/ - /***************************************************************************************/ - - // - //-------- second model calculation via shifted FFT ------------------------ - // - // calculation of the spectral power via FFT - PolarSpec1024 ( &data->L[576], erg0, phs0 ); // left - PolarSpec1024 ( &data->R[576], erg1, phs1 ); // right - - // calculation of the acoustic pressures per each subband for L/R-signals - SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 ); - - // calculation of the acoustic pressures per each partition - PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 ); - - // calculate the predictability of the signal - // left - memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) ); - memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) ); - CalcUnpred ( MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L ); - // right - memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) ); - memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) ); - CalcUnpred ( MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R ); - - // calculation of the weighted acoustic pressure per each partition - WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R ); - - // Spreading Signal & weighted unpredictability-signal - // left - memset ( clow_L , 0, sizeof clow_L ); - memset ( sim_Mask_L, 0, sizeof sim_Mask_L ); - SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L ); - // right - memset ( clow_R , 0, sizeof clow_R ); - memset ( sim_Mask_R, 0, sizeof sim_Mask_R ); - SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R ); - - // Offset depending on tonality - ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R ); - - // Handling of transient signals - // calculate four short FFTs (left) - PowSpec256 ( &data->L[ 576+SHORTFFT_OFFSET], F_256[0] ); - PowSpec256 ( &data->L[ 720+SHORTFFT_OFFSET], F_256[1] ); - PowSpec256 ( &data->L[ 864+SHORTFFT_OFFSET], F_256[2] ); - PowSpec256 ( &data->L[1008+SHORTFFT_OFFSET], F_256[3] ); - // calculate short Threshold - CalcShortThreshold ( F_256, ShortThr, shortThr_L, pre_erg_L, TransientL ); - - // calculate four short FFTs (right) - PowSpec256 ( &data->R[ 576+SHORTFFT_OFFSET], F_256[0] ); - PowSpec256 ( &data->R[ 720+SHORTFFT_OFFSET], F_256[1] ); - PowSpec256 ( &data->R[ 864+SHORTFFT_OFFSET], F_256[2] ); - PowSpec256 ( &data->R[1008+SHORTFFT_OFFSET], F_256[3] ); - // calculate short Threshold - CalcShortThreshold ( F_256, ShortThr, shortThr_R, pre_erg_R, TransientR ); - - // dynamic adjustment of threshold in quiet to loudness of the current sequence - if ( varLtq > 0. ) - factorLTQ = AdaptLtq ( Ls_L, Ls_R ); - - // utilization of temporal post-masking - if (tmpMask_used) { - CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L ); - CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R ); - memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L ); - memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R ); - } - - // transient signal? - for ( n = 0; n < PART_SHORT; n++ ) { - if ( TransientL[n] ) { - sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] ); - sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] ); - sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] ); - } - if ( TransientR[n] ) { - sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] ); - sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] ); - sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] ); - } - } - - // Pre-Echo control - PreechoControl ( PartThr_L, PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R ); - - // utilization of threshold in quiet - ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 ); - - // Consideration of aliasing between the subbands (noise is smeared) - // In: Thr[0..511], Out: Thr[512...1023] - AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 ); - memmove ( Thr_L, Thr_L+512, 512*sizeof(float) ); - memmove ( Thr_R, Thr_R+512, 512*sizeof(float) ); - - // calculation of the Signal-to-Mask-Ratio - CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR1.L, SMR1.R ); - - /***************************************************************************************/ - /***************************************************************************************/ - if ( MS_Channelmode > 0 ) { - // calculation of the spectral energy via FFT - PowSpec1024 ( &data->M[576], erg0 ); // mid - PowSpec1024 ( &data->S[576], erg1 ); // side - - // calculation of the acoustic pressure per each subband for M/S-signals - SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 ); - - // calculation of the acoustic pressure per each partition - PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 ); - - // calculate masking thresholds for M/S - CalcMSThreshold ( Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S ); - ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 ); - - // Consideration of aliasing between the subbands (noise is smeared) - // In: Thr[0..511], Out: Thr[512...1023] - AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 ); - memmove ( Thr_M, Thr_M+512, 512*sizeof(float) ); - memmove ( Thr_S, Thr_S+512, 512*sizeof(float) ); - - // calculation of the Signal-to-Mask-Ratio - CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR1.M, SMR1.S ); - } - /***************************************************************************************/ - /***************************************************************************************/ - - if ( NS_Order > 0 ) { - for ( n = 0; n < MAX_ANS_LINES; n++ ) { // providing Noise Shaping thresholds - ANSspec_L [n] = minf ( ANSspec_L [n], Thr_L [n] ); - ANSspec_R [n] = minf ( ANSspec_R [n], Thr_R [n] ); - ANSspec_M [n] = minf ( ANSspec_M [n], Thr_M [n] ); - ANSspec_S [n] = minf ( ANSspec_S [n], Thr_S [n] ); - } - } - - for ( n = 0; n <= MaxBand; n++ ) { // choose 'worst case'-SMR from shifted analysis windows - SMR0.L[n] = maxf ( SMR0.L[n], SMR1.L[n] ); - SMR0.R[n] = maxf ( SMR0.R[n], SMR1.R[n] ); - SMR0.M[n] = maxf ( SMR0.M[n], SMR1.M[n] ); - SMR0.S[n] = maxf ( SMR0.S[n], SMR1.S[n] ); - } - - LEAVE(50); - return SMR0; -} Index: mppenc/branches/r2d/src/psy_tab.c =================================================================== --- mppenc/branches/r2d/src/psy_tab.c (revision 58) +++ (revision ) @@ -1,462 +1,0 @@ -/* - * Musepack audio compression - * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - */ - -#include "mppenc.h" - -// Antialiasing for calculation of the subband power -const float Butfly [7] = { 0.5f, 0.2776f, 0.1176f, 0.0361f, 0.0075f, 0.000948f, 0.0000598f }; - -// Antialiasing for calculation of the masking thresholds -const float InvButfly [7] = { 2.f, 3.6023f, 8.5034f, 27.701f, 133.33f, 1054.852f, 16722.408f }; - -// w_low for long 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 -const int wl [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 38, 41, 44, 47, 50, 54, 58, 62, 67, 72, 78, 84, 91, 98, 106, 115, 124, 134, 145, 157, 170, 184, 199, 216, 234, 254, 276, 301, 329, 360, 396, 437, 485 }; -const int wh [PART_LONG] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 37, 40, 43, 46, 49, 53, 57, 61, 66, 71, 77, 83, 90, 97, 105, 114, 123, 133, 144, 156, 169, 183, 198, 215, 233, 253, 275, 300, 328, 359, 395, 436, 484, 511 }; -// Width: 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 4 4 4 5 5 6 6 7 7 8 9 9 10 11 12 13 14 15 17 18 20 22 25 28 31 36 41 48 27 - -// inverse partition-width for long -const float iw [PART_LONG] = { 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/3, 1.f/4, 1.f/4, 1.f/4, 1.f/5, 1.f/5, 1.f/6, 1.f/6, 1.f/7, 1.f/7, 1.f/8, 1.f/9, 1.f/9, 1.f/10, 1.f/11, 1.f/12, 1.f/13, 1.f/14, 1.f/15, 1.f/17, 1.f/18, 1.f/20, 1.f/22, 1.f/25, 1.f/28, 1.f/31, 1.f/36, 1.f/41, 1.f/48, 1.f/27 }; - -// w_low for short 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 -const int wl_short [PART_SHORT] = { 0, 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 23, 29, 36, 46, 59, 75, 99 }; -const int wh_short [PART_SHORT] = { 0, 1, 2, 3, 5, 6, 7, 9, 12, 14, 18, 23, 29, 36, 46, 58, 75, 99, 127 }; - -// inverse partition-width for short -const float iw_short [PART_SHORT] = { 1.f, 1.f, 1.f, 1.f, 1.f/2, 1.f/2, 1.f/2, 1.f/2, 1.f/3, 1.f/3, 1.f/4, 1.f/6, 1.f/7, 1.f/8, 1.f/11, 1.f/13, 1.f/17, 1.f/25, 1.f/29 }; - -/* -Nr. wl wh fl fh bl bh bm Nr. wl wh fl fh bl bh bm - 0: 0 0 0 0 0.000000 0.000000 0.000000 - 1: 1 1 43 43 0.425460 0.425460 0.425460 0: 0 0 0 0 0.000000 0.000000 0.000000 - 2: 2 2 86 86 0.850241 0.850241 0.850241 - - 3: 3 3 129 129 1.273448 1.273448 1.273448 - 4: 4 4 172 172 1.694205 1.694205 1.694205 1: 1 1 172 172 1.694205 1.694205 1.694205 - 5: 5 5 215 215 2.111672 2.111672 2.111672 - - 6: 6 6 258 258 2.525051 2.525051 2.525051 - 7: 7 7 301 301 2.933594 2.933594 2.933594 2: 2 2 345 345 3.336612 3.336612 3.336612 - 8: 8 8 345 345 3.336612 3.336612 3.336612 - - 9: 9 9 388 388 3.733479 3.733479 3.733479 -10: 10 10 431 431 4.123635 4.123635 4.123635 3: 3 3 517 517 4.881924 4.881924 4.881924 -11- 11 12 474 517 4.506591 4.881924 4.695234 - -12: 13 14 560 603 5.249283 5.608381 5.429880 -13: 15 16 646 689 5.958998 6.300971 6.131073 4: 4 5 689 861 6.300971 7.581073 6.958618 -14: 17 18 732 775 6.634195 6.958618 6.797509 - -15: 19 20 818 861 7.274232 7.581073 7.428745 -16: 21 22 904 947 7.879211 8.168753 8.025049 5: 5 6 861 1034 7.581073 8.722594 8.168753 -17: 23 24 991 1034 8.449828 8.722594 8.587239 - -18: 25 26 1077 1120 8.987223 9.243908 9.116546 -19: 27 28 1163 1206 9.492850 9.734263 9.614484 6: 6 7 1034 1206 8.722594 9.734263 9.243908 -20: 29 30 1249 1292 9.968365 10.195382 10.082745 - -21: 31 32 1335 1378 10.415539 10.629064 10.523116 -22: 33 34 1421 1464 10.836184 11.037125 10.937413 7: 8 9 1378 1550 10.629064 11.421352 11.037125 -23: 35 37 1507 1593 11.232108 11.605071 11.421352 - -24: 38 40 1637 1723 11.783474 12.125139 11.956764 -25: 41 43 1766 1852 12.288791 12.602659 12.447904 8: 10 12 1723 2067 12.125139 13.316883 12.753228 -26: 44 46 1895 1981 12.753228 13.042468 12.899777 - -27: 47 49 2024 2110 13.181453 13.448898 13.316883 -28: 50 53 2153 2283 13.577635 13.945465 13.764881 9: 12 14 2067 2412 13.316883 14.288198 13.825796 -29: 54 57 2326 2455 14.062349 14.397371 14.232693 - -30: 58 61 2498 2627 14.504172 14.811258 14.660130 -31: 62 66 2670 2842 14.909464 15.283564 15.100115 10: 15 18 2584 3101 14.711029 15.795819 15.283564 -32: 67 71 2885 3058 15.372757 15.714074 15.546390 - -33: 72 77 3101 3316 15.795819 16.185532 15.994471 -34: 78 83 3359 3575 16.259980 16.616871 16.441494 11: 18 23 3101 3962 15.795819 17.204658 16.547424 -35: 84 90 3618 3876 16.685418 17.079349 16.885941 - -36: 91 97 3919 4177 17.142352 17.506445 17.327264 -37: 98 105 4221 4522 17.564981 17.959646 17.765487 12: 23 29 3962 4996 17.204658 18.533945 17.904788 -38: 106 114 4565 4910 18.014031 18.433233 18.227034 - -39: 115 123 4953 5297 18.483782 18.874805 18.682185 -40: 124 133 5340 5728 18.922095 19.332992 19.130789 13: 29 36 4996 6202 18.533945 19.801451 19.198897 -41: 134 144 5771 6202 19.377073 19.801451 19.592946 - -42: 145 156 6245 6718 19.842285 20.272889 20.061808 -43: 157 169 6761 7278 20.310373 20.739167 20.529583 14: 36 46 6202 7924 19.801451 21.222342 20.565177 -44: 170 183 7321 7881 20.773175 21.191895 20.987911 - -45: 184 198 7924 8527 21.222342 21.623344 21.428652 -46: 199 215 8570 9259 21.650236 22.050787 21.857360 15: 46 58 7924 9991 21.222342 22.420001 21.882271 -47: 216 233 9302 10034 22.074042 22.440072 22.263795 - -48- 234 253 10078 10896 22.459969 22.807140 22.640652 -49: 254 275 10939 11843 22.823891 23.144847 22.991444 16: 59 75 10164 12920 22.499251 23.461146 23.044078 -50: 276 300 11886 12920 23.158772 23.461146 23.317264 - -51: 301 328 12963 14126 23.472530 23.748999 23.617861 -52: 329 359 14169 15461 23.758199 24.005540 23.888450 17: 75 99 12920 17054 23.461146 24.248491 23.920884 -53: 360 395 15504 17011 24.012922 24.242660 24.134368 - -54: 396 436 17054 18777 24.248491 24.454928 24.357873 -55: 437 484 18820 20844 24.459492 24.647977 24.559711 18: 99 127 17054 21878 24.248491 24.727775 24.524955 -56: 485 511 20887 22007 24.651498 24.737100 24.695685 -*/ - - -/* V A R I A B L E S */ -float MinVal [PART_LONG]; // contains minimum tonality soffsets -float Loudness [PART_LONG]; // weighting factors for loudness calculation -float SPRD [PART_LONG] [PART_LONG]; // tabulated spreading function -float O_MAX; -float O_MIN; -float FAC1; -float FAC2; // constants for offset calculation -float partLtq [PART_LONG]; // threshold in quiet (partitions) -float invLtq [PART_LONG]; // inverse threshold in quiet (partitions, long) -float fftLtq [512]; // threshold in quiet (FFT) -float Ltq_offset; // Offset for threshold in quiet -float Ltq_max; // maximum level for threshold in quiet -float TMN; -float NMT; -float TransDetect; -unsigned int EarModelFlag; -int MinValChoice; - - -/* - * Klemm 1994 and 1997. Experimental data. Sorry, data looks a little bit - * dodderly. Data below 30 Hz is extrapolated from other material, above 18 - * kHz the ATH is limited due to the original purpose (too much noise at - * ATH is not good even if it's theoretically inaudible). - */ - -static float -ATHformula_Frank ( float freq ) -{ - /* - * one value per 100 cent = 1 - * semitone = 1/4 - * third = 1/12 - * octave = 1/40 decade - * rest is linear interpolated, values are currently in millibel rel. 20 µPa - */ - static short tab [] = { - /* 10.0 */ 9669, 9669, 9626, 9512, - /* 12.6 */ 9353, 9113, 8882, 8676, - /* 15.8 */ 8469, 8243, 7997, 7748, - /* 20.0 */ 7492, 7239, 7000, 6762, - /* 25.1 */ 6529, 6302, 6084, 5900, - /* 31.6 */ 5717, 5534, 5351, 5167, - /* 39.8 */ 5004, 4812, 4638, 4466, - /* 50.1 */ 4310, 4173, 4050, 3922, - /* 63.1 */ 3723, 3577, 3451, 3281, - /* 79.4 */ 3132, 3036, 2902, 2760, - /* 100.0 */ 2658, 2591, 2441, 2301, - /* 125.9 */ 2212, 2125, 2018, 1900, - /* 158.5 */ 1770, 1682, 1594, 1512, - /* 199.5 */ 1430, 1341, 1260, 1198, - /* 251.2 */ 1136, 1057, 998, 943, - /* 316.2 */ 887, 846, 744, 712, - /* 398.1 */ 693, 668, 637, 606, - /* 501.2 */ 580, 555, 529, 502, - /* 631.0 */ 475, 448, 422, 398, - /* 794.3 */ 375, 351, 327, 322, - /* 1000.0 */ 312, 301, 291, 268, - /* 1258.9 */ 246, 215, 182, 146, - /* 1584.9 */ 107, 61, 13, -35, - /* 1995.3 */ -96, -156, -179, -235, - /* 2511.9 */ -295, -350, -401, -421, - /* 3162.3 */ -446, -499, -532, -535, - /* 3981.1 */ -513, -476, -431, -313, - /* 5011.9 */ -179, 8, 203, 403, - /* 6309.6 */ 580, 736, 881, 1022, - /* 7943.3 */ 1154, 1251, 1348, 1421, - /* 10000.0 */ 1479, 1399, 1285, 1193, - /* 12589.3 */ 1287, 1519, 1914, 2369, -#if 0 - /* 15848.9 */ 3352, 4865, 5942, 6177, - /* 19952.6 */ 6385, 6604, 6833, 7009, - /* 25118.9 */ 7066, 7127, 7191, 7260, -#else - /* 15848.9 */ 3352, 4352, 5352, 6352, - /* 19952.6 */ 7352, 8352, 9352, 9999, - /* 25118.9 */ 9999, 9999, 9999, 9999, -#endif - }; - double freq_log; - unsigned index; - - if ( freq < 10. ) freq = 10.; - if ( freq > 29853. ) freq = 29853.; - - freq_log = 40. * log10 (0.1 * freq); /* 4 steps per third, starting at 10 Hz */ - index = (unsigned) freq_log; - return 0.01 * (tab [index] * (1 + index - freq_log) + tab [index+1] * (freq_log - index)); -} - - -/* F U N C T I O N S */ -// calculation of the threshold in quiet in FFT-resolution -static void -Ruhehoerschwelle ( unsigned int EarModelFlag, - int Ltq_offset, - int Ltq_max ) -{ - int n; - int k; - float f; - float erg; - double tmp; - float absLtq [512]; - - for ( n = 0; n < 512; n++ ) { - f = (float) ( (n+1) * (float)(SampleFreq / 2000.) / 512 ); // Frequency in kHz - - switch ( EarModelFlag / 100 ) { - case 0: // ISO-threshold in quiet - tmp = 3.64*pow (f,-0.8) - 6.5*exp (-0.6*(f-3.3)*(f-3.3)) + 0.001*pow (f, 4.0); - break; - default: - case 1: // measured threshold in quiet (Nick Berglmeir, Andree Buschmann, Kopfhörer) - tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); - break; - case 2: // measured threshold in quiet (Filburt, Kopfhörer) - tmp = 9.00*pow (f,-0.5) - 15.0*exp (-0.1*(f-4.0)*(f-4.0)) + 0.0341796875*pow (f, 2.5) + 15.*exp (-(f-0.1)*(f-0.1)) - 18; - tmp = mind ( tmp, Ltq_max - 18 ); - break; - case 3: - tmp = ATHformula_Frank ( 1.e3 * f ); - break; - case 4: - tmp = ATHformula_Frank ( 1.e3 * f ); - if ( f > 4.8 ) { - tmp += 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); - tmp *= 0.5 ; - } - break; - case 5: - tmp = ATHformula_Frank ( 1.e3 * f ); - if ( f > 4.8 ) { - tmp = 3.00*pow (f,-0.8) - 5.0*exp (-0.1*(f-3.0)*(f-3.0)) + 0.0000015022693846297*pow (f, 6.0) + 10.*exp (-(f-0.1)*(f-0.1)); - } - break; - } - - tmp -= f * f * (int)(EarModelFlag % 100 - 50) * 0.0015; // 00: +30 dB, 100: -30 dB @20 kHz - - tmp = mind ( tmp, Ltq_max ); // Limit ATH - tmp += Ltq_offset - 23; // Add chosen Offset - fftLtq[n] = absLtq[n] = POW10 ( 0.1 * tmp); // conversion into power - } - - // threshold in quiet in partitions (long) - for ( n = 0; n < PART_LONG; n++ ) { - erg = 1.e20f; - for ( k = wl[n]; k <= wh[n]; k++ ) - erg = minf (erg, absLtq[k]); - - partLtq[n] = erg; // threshold in quiet - invLtq [n] = 1.f / partLtq[n]; // Inverse - } -} - -#ifdef _WIN32 -static double -asinh ( double x ) -{ - return x >= 0 ? log (sqrt (x*x+1) + x) : -log (sqrt (x*x+1) - x); -} -#endif - - -static double -Freq2Bark ( double Hz ) // Klemm 2002 -{ - return 9.97074*asinh (1.1268e-3 * Hz) - 6.25817*asinh (0.197193e-3 * Hz) ; -} - -static double -Bark2Freq ( double Bark ) // Klemm 2002 -{ - return 956.86 * sinh (0.101561*Bark) + 11.7296 * sinh (0.304992*Bark) + 6.33622e-3*sinh (0.538621*Bark); -} - -static double -LongPart2Bark ( int Part ) -{ - return Freq2Bark ((wl [Part] + wh [Part]) * SampleFreq / 2048.); -} - -// calculating the table for loudness calculation based on absLtq = ank -static void -Loudness_Tabelle (void) -{ - int n; - float midfreq; - float tmp; - - // ca. dB(A) - for ( n = 0; n < PART_LONG; n++ ){ - midfreq = (wh[n] + wl[n] + 3) * (0.25 * SampleFreq / 512); // center frequency in kHz, why +3 ??? - tmp = LOG10 (midfreq) - 3.5f; // dB(A) - tmp = -10 * tmp * tmp + 3 - midfreq/3000; - Loudness [n] = POW10 ( 0.1 * tmp ); // conversion into power - } -} - - -static double -Bass ( float f, float TMN, float NMT, float bass ) -{ - static unsigned char lfe [11] = { 120, 100, 80, 60, 50, 40, 30, 20, 15, 10, 5 }; - int tmp = (int) ( 1024/44100. * f + 0.5 ); - - switch ( tmp ) { - case 0: - case 1: - case 2: - case 3: - case 4: - case 5: - case 6: - case 7: - case 8: - case 9: - case 10: - return TMN + bass * lfe [tmp]; - case 11: - case 12: - case 13: - case 14: - case 15: - case 16: - case 17: - case 18: - return TMN; - case 19: - case 20: - case 21: - case 22: - return TMN*0.75 + NMT*0.25; - case 23: - case 24: - return TMN*0.50 + NMT*0.50; - case 25: - case 26: - return TMN*0.25 + NMT*0.75; - default: - return NMT; - } -} - - -// calculating the coefficient for utilization of the tonality offset, depending on TMN und NMT -static void -Tonalitaetskoeffizienten ( void ) -{ - double tmp; - int n; - float bass; - - bass = 0.1/8 * NMT; - if ( MinValChoice <= 2 && bass > 0.1 ) - bass = 0.1f; - if ( MinValChoice <= 1 ) - bass = 0.0f; - - // alternative: calculation of the minval-values dependent on TMN and TMN - for ( n = 0; n < PART_LONG; n++ ) { - tmp = Bass ( (wl [n] + wh [n]) / 2048. * SampleFreq, TMN, NMT, bass ); - MinVal [n] = POW10 ( -0.1 * tmp ); // conversion into power - } - - // calculation of the constants for "tonality offset" - O_MAX = POW10 ( -0.1 * TMN ); - O_MIN = POW10 ( -0.1 * NMT ); - FAC1 = POW10 ( -0.1 * (NMT - (TMN - NMT) * 0.229) ) ; - FAC2 = (TMN - NMT) * (0.99011159 * 0.1); -} - - -// calculation of the spreading function -static void -Spread ( void ) -{ - int i; - int j; - float tmpx; - float tmpy; - float tmpz; - float x; - - // calculation of the spreading-function for all occuring values - for ( i = 0; i < PART_LONG; i++ ) { // i is masking Partition, Source - for ( j = 0; j < PART_LONG; j++ ) { // j is masking Partition, Target - tmpx = LongPart2Bark (j) - LongPart2Bark (i);// Difference of the partitions in Bark - tmpy = tmpz = 0.; // tmpz = 0: no dip - - if ( tmpx < 0 ) { // downwards (S1) - tmpy = -32.f * tmpx; // 32 dB per Bark, e33 (10) - } - else if ( tmpx > 0 ) { // upwards (S2) -#if 0 - x = (wl[i]+wh[i])/2 * (float)(SampleFreq / 2000)/512; // center frequency in kHz ??????? - if (i==0) x = 0.5f * (float)(SampleFreq / 2000)/512; // if first spectral line -#else - x = i ? wl[i]+wh[i] : 1; - x *= SampleFreq / 1000. / 2048; // center frequency in kHz -#endif - // dB/Bark - tmpy = (22.f + 0.23f / x) * tmpx; // e33 (10) - - // dip (up to 6 dB) - tmpz = 8 * minf ( (tmpx-0.5f) * (tmpx-0.5f) - 2 * (tmpx-0.5f), 0.f ); - } - - // calculate coefficient - SPRD[i][j] = POW10 ( -0.1 * (tmpy+tmpz) ); // [Source] [Target] - } - } - - // Normierung e33 (10) - for ( i = 0; i < PART_LONG; i++ ) { // i is masked Partition - float norm = 0.f; - for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition - norm += SPRD [j] [i]; - for ( j = 0; j < PART_LONG; j++ ) // j is masking Partition - SPRD [j] [i] /= norm; - } -} - -// call all initialisation procedures -void -Init_Psychoakustiktabellen ( void ) -{ - Max_Band = (int) ( Bandwidth * 64. / SampleFreq ); - if ( Max_Band < 1 ) Max_Band = 1; - if ( Max_Band > 31 ) Max_Band = 31; - - Tonalitaetskoeffizienten (); - Ruhehoerschwelle ( EarModelFlag, Ltq_offset, Ltq_max ); - Loudness_Tabelle (); - Spread (); -} - -/* end of psy_tab.c */