psy.c 35 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207
  1. /********************************************************************
  2. * *
  3. * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
  4. * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
  5. * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
  6. * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
  7. * *
  8. * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2010 *
  9. * by the Xiph.Org Foundation http://www.xiph.org/ *
  10. * *
  11. ********************************************************************
  12. function: psychoacoustics not including preecho
  13. last mod: $Id: psy.c 18077 2011-09-02 02:49:00Z giles $
  14. ********************************************************************/
  15. #include <stdlib.h>
  16. #include <math.h>
  17. #include <string.h>
  18. #include "vorbis/codec.h"
  19. #include "codec_internal.h"
  20. #include "masking.h"
  21. #include "psy.h"
  22. #include "os.h"
  23. #include "lpc.h"
  24. #include "smallft.h"
  25. #include "scales.h"
  26. #include "misc.h"
  27. #define NEGINF -9999.f
  28. static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
  29. static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10};
  30. vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
  31. codec_setup_info *ci=vi->codec_setup;
  32. vorbis_info_psy_global *gi=&ci->psy_g_param;
  33. vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
  34. look->channels=vi->channels;
  35. look->ampmax=-9999.;
  36. look->gi=gi;
  37. return(look);
  38. }
  39. void _vp_global_free(vorbis_look_psy_global *look){
  40. if(look){
  41. memset(look,0,sizeof(*look));
  42. _ogg_free(look);
  43. }
  44. }
  45. void _vi_gpsy_free(vorbis_info_psy_global *i){
  46. if(i){
  47. memset(i,0,sizeof(*i));
  48. _ogg_free(i);
  49. }
  50. }
  51. void _vi_psy_free(vorbis_info_psy *i){
  52. if(i){
  53. memset(i,0,sizeof(*i));
  54. _ogg_free(i);
  55. }
  56. }
  57. static void min_curve(float *c,
  58. float *c2){
  59. int i;
  60. for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
  61. }
  62. static void max_curve(float *c,
  63. float *c2){
  64. int i;
  65. for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
  66. }
  67. static void attenuate_curve(float *c,float att){
  68. int i;
  69. for(i=0;i<EHMER_MAX;i++)
  70. c[i]+=att;
  71. }
  72. static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
  73. float center_boost, float center_decay_rate){
  74. int i,j,k,m;
  75. float ath[EHMER_MAX];
  76. float workc[P_BANDS][P_LEVELS][EHMER_MAX];
  77. float athc[P_LEVELS][EHMER_MAX];
  78. float *brute_buffer=alloca(n*sizeof(*brute_buffer));
  79. float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);
  80. memset(workc,0,sizeof(workc));
  81. for(i=0;i<P_BANDS;i++){
  82. /* we add back in the ATH to avoid low level curves falling off to
  83. -infinity and unnecessarily cutting off high level curves in the
  84. curve limiting (last step). */
  85. /* A half-band's settings must be valid over the whole band, and
  86. it's better to mask too little than too much */
  87. int ath_offset=i*4;
  88. for(j=0;j<EHMER_MAX;j++){
  89. float min=999.;
  90. for(k=0;k<4;k++)
  91. if(j+k+ath_offset<MAX_ATH){
  92. if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
  93. }else{
  94. if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
  95. }
  96. ath[j]=min;
  97. }
  98. /* copy curves into working space, replicate the 50dB curve to 30
  99. and 40, replicate the 100dB curve to 110 */
  100. for(j=0;j<6;j++)
  101. memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
  102. memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
  103. memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
  104. /* apply centered curve boost/decay */
  105. for(j=0;j<P_LEVELS;j++){
  106. for(k=0;k<EHMER_MAX;k++){
  107. float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
  108. if(adj<0. && center_boost>0)adj=0.;
  109. if(adj>0. && center_boost<0)adj=0.;
  110. workc[i][j][k]+=adj;
  111. }
  112. }
  113. /* normalize curves so the driving amplitude is 0dB */
  114. /* make temp curves with the ATH overlayed */
  115. for(j=0;j<P_LEVELS;j++){
  116. attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0);
  117. memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
  118. attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0);
  119. max_curve(athc[j],workc[i][j]);
  120. }
  121. /* Now limit the louder curves.
  122. the idea is this: We don't know what the playback attenuation
  123. will be; 0dB SL moves every time the user twiddles the volume
  124. knob. So that means we have to use a single 'most pessimal' curve
  125. for all masking amplitudes, right? Wrong. The *loudest* sound
  126. can be in (we assume) a range of ...+100dB] SL. However, sounds
  127. 20dB down will be in a range ...+80], 40dB down is from ...+60],
  128. etc... */
  129. for(j=1;j<P_LEVELS;j++){
  130. min_curve(athc[j],athc[j-1]);
  131. min_curve(workc[i][j],athc[j]);
  132. }
  133. }
  134. for(i=0;i<P_BANDS;i++){
  135. int hi_curve,lo_curve,bin;
  136. ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);
  137. /* low frequency curves are measured with greater resolution than
  138. the MDCT/FFT will actually give us; we want the curve applied
  139. to the tone data to be pessimistic and thus apply the minimum
  140. masking possible for a given bin. That means that a single bin
  141. could span more than one octave and that the curve will be a
  142. composite of multiple octaves. It also may mean that a single
  143. bin may span > an eighth of an octave and that the eighth
  144. octave values may also be composited. */
  145. /* which octave curves will we be compositing? */
  146. bin=floor(fromOC(i*.5)/binHz);
  147. lo_curve= ceil(toOC(bin*binHz+1)*2);
  148. hi_curve= floor(toOC((bin+1)*binHz)*2);
  149. if(lo_curve>i)lo_curve=i;
  150. if(lo_curve<0)lo_curve=0;
  151. if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;
  152. for(m=0;m<P_LEVELS;m++){
  153. ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));
  154. for(j=0;j<n;j++)brute_buffer[j]=999.;
  155. /* render the curve into bins, then pull values back into curve.
  156. The point is that any inherent subsampling aliasing results in
  157. a safe minimum */
  158. for(k=lo_curve;k<=hi_curve;k++){
  159. int l=0;
  160. for(j=0;j<EHMER_MAX;j++){
  161. int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz;
  162. int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1;
  163. if(lo_bin<0)lo_bin=0;
  164. if(lo_bin>n)lo_bin=n;
  165. if(lo_bin<l)l=lo_bin;
  166. if(hi_bin<0)hi_bin=0;
  167. if(hi_bin>n)hi_bin=n;
  168. for(;l<hi_bin && l<n;l++)
  169. if(brute_buffer[l]>workc[k][m][j])
  170. brute_buffer[l]=workc[k][m][j];
  171. }
  172. for(;l<n;l++)
  173. if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
  174. brute_buffer[l]=workc[k][m][EHMER_MAX-1];
  175. }
  176. /* be equally paranoid about being valid up to next half ocatve */
  177. if(i+1<P_BANDS){
  178. int l=0;
  179. k=i+1;
  180. for(j=0;j<EHMER_MAX;j++){
  181. int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz;
  182. int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1;
  183. if(lo_bin<0)lo_bin=0;
  184. if(lo_bin>n)lo_bin=n;
  185. if(lo_bin<l)l=lo_bin;
  186. if(hi_bin<0)hi_bin=0;
  187. if(hi_bin>n)hi_bin=n;
  188. for(;l<hi_bin && l<n;l++)
  189. if(brute_buffer[l]>workc[k][m][j])
  190. brute_buffer[l]=workc[k][m][j];
  191. }
  192. for(;l<n;l++)
  193. if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
  194. brute_buffer[l]=workc[k][m][EHMER_MAX-1];
  195. }
  196. for(j=0;j<EHMER_MAX;j++){
  197. int bin=fromOC(j*.125+i*.5-2.)/binHz;
  198. if(bin<0){
  199. ret[i][m][j+2]=-999.;
  200. }else{
  201. if(bin>=n){
  202. ret[i][m][j+2]=-999.;
  203. }else{
  204. ret[i][m][j+2]=brute_buffer[bin];
  205. }
  206. }
  207. }
  208. /* add fenceposts */
  209. for(j=0;j<EHMER_OFFSET;j++)
  210. if(ret[i][m][j+2]>-200.f)break;
  211. ret[i][m][0]=j;
  212. for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
  213. if(ret[i][m][j+2]>-200.f)
  214. break;
  215. ret[i][m][1]=j;
  216. }
  217. }
  218. return(ret);
  219. }
  220. void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
  221. vorbis_info_psy_global *gi,int n,long rate){
  222. long i,j,lo=-99,hi=1;
  223. long maxoc;
  224. memset(p,0,sizeof(*p));
  225. p->eighth_octave_lines=gi->eighth_octave_lines;
  226. p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
  227. p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
  228. maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f;
  229. p->total_octave_lines=maxoc-p->firstoc+1;
  230. p->ath=_ogg_malloc(n*sizeof(*p->ath));
  231. p->octave=_ogg_malloc(n*sizeof(*p->octave));
  232. p->bark=_ogg_malloc(n*sizeof(*p->bark));
  233. p->vi=vi;
  234. p->n=n;
  235. p->rate=rate;
  236. /* AoTuV HF weighting */
  237. p->m_val = 1.;
  238. if(rate < 26000) p->m_val = 0;
  239. else if(rate < 38000) p->m_val = .94; /* 32kHz */
  240. else if(rate > 46000) p->m_val = 1.275; /* 48kHz */
  241. /* set up the lookups for a given blocksize and sample rate */
  242. for(i=0,j=0;i<MAX_ATH-1;i++){
  243. int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
  244. float base=ATH[i];
  245. if(j<endpos){
  246. float delta=(ATH[i+1]-base)/(endpos-j);
  247. for(;j<endpos && j<n;j++){
  248. p->ath[j]=base+100.;
  249. base+=delta;
  250. }
  251. }
  252. }
  253. for(;j<n;j++){
  254. p->ath[j]=p->ath[j-1];
  255. }
  256. for(i=0;i<n;i++){
  257. float bark=toBARK(rate/(2*n)*i);
  258. for(;lo+vi->noisewindowlomin<i &&
  259. toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
  260. for(;hi<=n && (hi<i+vi->noisewindowhimin ||
  261. toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
  262. p->bark[i]=((lo-1)<<16)+(hi-1);
  263. }
  264. for(i=0;i<n;i++)
  265. p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
  266. p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
  267. vi->tone_centerboost,vi->tone_decay);
  268. /* set up rolling noise median */
  269. p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
  270. for(i=0;i<P_NOISECURVES;i++)
  271. p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
  272. for(i=0;i<n;i++){
  273. float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
  274. int inthalfoc;
  275. float del;
  276. if(halfoc<0)halfoc=0;
  277. if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
  278. inthalfoc=(int)halfoc;
  279. del=halfoc-inthalfoc;
  280. for(j=0;j<P_NOISECURVES;j++)
  281. p->noiseoffset[j][i]=
  282. p->vi->noiseoff[j][inthalfoc]*(1.-del) +
  283. p->vi->noiseoff[j][inthalfoc+1]*del;
  284. }
  285. #if 0
  286. {
  287. static int ls=0;
  288. _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
  289. _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
  290. _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
  291. }
  292. #endif
  293. }
  294. void _vp_psy_clear(vorbis_look_psy *p){
  295. int i,j;
  296. if(p){
  297. if(p->ath)_ogg_free(p->ath);
  298. if(p->octave)_ogg_free(p->octave);
  299. if(p->bark)_ogg_free(p->bark);
  300. if(p->tonecurves){
  301. for(i=0;i<P_BANDS;i++){
  302. for(j=0;j<P_LEVELS;j++){
  303. _ogg_free(p->tonecurves[i][j]);
  304. }
  305. _ogg_free(p->tonecurves[i]);
  306. }
  307. _ogg_free(p->tonecurves);
  308. }
  309. if(p->noiseoffset){
  310. for(i=0;i<P_NOISECURVES;i++){
  311. _ogg_free(p->noiseoffset[i]);
  312. }
  313. _ogg_free(p->noiseoffset);
  314. }
  315. memset(p,0,sizeof(*p));
  316. }
  317. }
  318. /* octave/(8*eighth_octave_lines) x scale and dB y scale */
  319. static void seed_curve(float *seed,
  320. const float **curves,
  321. float amp,
  322. int oc, int n,
  323. int linesper,float dBoffset){
  324. int i,post1;
  325. int seedptr;
  326. const float *posts,*curve;
  327. int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
  328. choice=max(choice,0);
  329. choice=min(choice,P_LEVELS-1);
  330. posts=curves[choice];
  331. curve=posts+2;
  332. post1=(int)posts[1];
  333. seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
  334. for(i=posts[0];i<post1;i++){
  335. if(seedptr>0){
  336. float lin=amp+curve[i];
  337. if(seed[seedptr]<lin)seed[seedptr]=lin;
  338. }
  339. seedptr+=linesper;
  340. if(seedptr>=n)break;
  341. }
  342. }
  343. static void seed_loop(vorbis_look_psy *p,
  344. const float ***curves,
  345. const float *f,
  346. const float *flr,
  347. float *seed,
  348. float specmax){
  349. vorbis_info_psy *vi=p->vi;
  350. long n=p->n,i;
  351. float dBoffset=vi->max_curve_dB-specmax;
  352. /* prime the working vector with peak values */
  353. for(i=0;i<n;i++){
  354. float max=f[i];
  355. long oc=p->octave[i];
  356. while(i+1<n && p->octave[i+1]==oc){
  357. i++;
  358. if(f[i]>max)max=f[i];
  359. }
  360. if(max+6.f>flr[i]){
  361. oc=oc>>p->shiftoc;
  362. if(oc>=P_BANDS)oc=P_BANDS-1;
  363. if(oc<0)oc=0;
  364. seed_curve(seed,
  365. curves[oc],
  366. max,
  367. p->octave[i]-p->firstoc,
  368. p->total_octave_lines,
  369. p->eighth_octave_lines,
  370. dBoffset);
  371. }
  372. }
  373. }
  374. static void seed_chase(float *seeds, int linesper, long n){
  375. long *posstack=alloca(n*sizeof(*posstack));
  376. float *ampstack=alloca(n*sizeof(*ampstack));
  377. long stack=0;
  378. long pos=0;
  379. long i;
  380. for(i=0;i<n;i++){
  381. if(stack<2){
  382. posstack[stack]=i;
  383. ampstack[stack++]=seeds[i];
  384. }else{
  385. while(1){
  386. if(seeds[i]<ampstack[stack-1]){
  387. posstack[stack]=i;
  388. ampstack[stack++]=seeds[i];
  389. break;
  390. }else{
  391. if(i<posstack[stack-1]+linesper){
  392. if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
  393. i<posstack[stack-2]+linesper){
  394. /* we completely overlap, making stack-1 irrelevant. pop it */
  395. stack--;
  396. continue;
  397. }
  398. }
  399. posstack[stack]=i;
  400. ampstack[stack++]=seeds[i];
  401. break;
  402. }
  403. }
  404. }
  405. }
  406. /* the stack now contains only the positions that are relevant. Scan
  407. 'em straight through */
  408. for(i=0;i<stack;i++){
  409. long endpos;
  410. if(i<stack-1 && ampstack[i+1]>ampstack[i]){
  411. endpos=posstack[i+1];
  412. }else{
  413. endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
  414. discarded in short frames */
  415. }
  416. if(endpos>n)endpos=n;
  417. for(;pos<endpos;pos++)
  418. seeds[pos]=ampstack[i];
  419. }
  420. /* there. Linear time. I now remember this was on a problem set I
  421. had in Grad Skool... I didn't solve it at the time ;-) */
  422. }
  423. /* bleaugh, this is more complicated than it needs to be */
  424. #include<stdio.h>
  425. static void max_seeds(vorbis_look_psy *p,
  426. float *seed,
  427. float *flr){
  428. long n=p->total_octave_lines;
  429. int linesper=p->eighth_octave_lines;
  430. long linpos=0;
  431. long pos;
  432. seed_chase(seed,linesper,n); /* for masking */
  433. pos=p->octave[0]-p->firstoc-(linesper>>1);
  434. while(linpos+1<p->n){
  435. float minV=seed[pos];
  436. long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
  437. if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
  438. while(pos+1<=end){
  439. pos++;
  440. if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
  441. minV=seed[pos];
  442. }
  443. end=pos+p->firstoc;
  444. for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
  445. if(flr[linpos]<minV)flr[linpos]=minV;
  446. }
  447. {
  448. float minV=seed[p->total_octave_lines-1];
  449. for(;linpos<p->n;linpos++)
  450. if(flr[linpos]<minV)flr[linpos]=minV;
  451. }
  452. }
  453. static void bark_noise_hybridmp(int n,const long *b,
  454. const float *f,
  455. float *noise,
  456. const float offset,
  457. const int fixed){
  458. float *N=alloca(n*sizeof(*N));
  459. float *X=alloca(n*sizeof(*N));
  460. float *XX=alloca(n*sizeof(*N));
  461. float *Y=alloca(n*sizeof(*N));
  462. float *XY=alloca(n*sizeof(*N));
  463. float tN, tX, tXX, tY, tXY;
  464. int i;
  465. int lo, hi;
  466. float R=0.f;
  467. float A=0.f;
  468. float B=0.f;
  469. float D=1.f;
  470. float w, x, y;
  471. tN = tX = tXX = tY = tXY = 0.f;
  472. y = f[0] + offset;
  473. if (y < 1.f) y = 1.f;
  474. w = y * y * .5;
  475. tN += w;
  476. tX += w;
  477. tY += w * y;
  478. N[0] = tN;
  479. X[0] = tX;
  480. XX[0] = tXX;
  481. Y[0] = tY;
  482. XY[0] = tXY;
  483. for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
  484. y = f[i] + offset;
  485. if (y < 1.f) y = 1.f;
  486. w = y * y;
  487. tN += w;
  488. tX += w * x;
  489. tXX += w * x * x;
  490. tY += w * y;
  491. tXY += w * x * y;
  492. N[i] = tN;
  493. X[i] = tX;
  494. XX[i] = tXX;
  495. Y[i] = tY;
  496. XY[i] = tXY;
  497. }
  498. for (i = 0, x = 0.f;; i++, x += 1.f) {
  499. lo = b[i] >> 16;
  500. if( lo>=0 ) break;
  501. hi = b[i] & 0xffff;
  502. tN = N[hi] + N[-lo];
  503. tX = X[hi] - X[-lo];
  504. tXX = XX[hi] + XX[-lo];
  505. tY = Y[hi] + Y[-lo];
  506. tXY = XY[hi] - XY[-lo];
  507. A = tY * tXX - tX * tXY;
  508. B = tN * tXY - tX * tY;
  509. D = tN * tXX - tX * tX;
  510. R = (A + x * B) / D;
  511. if (R < 0.f)
  512. R = 0.f;
  513. noise[i] = R - offset;
  514. }
  515. for ( ;; i++, x += 1.f) {
  516. lo = b[i] >> 16;
  517. hi = b[i] & 0xffff;
  518. if(hi>=n)break;
  519. tN = N[hi] - N[lo];
  520. tX = X[hi] - X[lo];
  521. tXX = XX[hi] - XX[lo];
  522. tY = Y[hi] - Y[lo];
  523. tXY = XY[hi] - XY[lo];
  524. A = tY * tXX - tX * tXY;
  525. B = tN * tXY - tX * tY;
  526. D = tN * tXX - tX * tX;
  527. R = (A + x * B) / D;
  528. if (R < 0.f) R = 0.f;
  529. noise[i] = R - offset;
  530. }
  531. for ( ; i < n; i++, x += 1.f) {
  532. R = (A + x * B) / D;
  533. if (R < 0.f) R = 0.f;
  534. noise[i] = R - offset;
  535. }
  536. if (fixed <= 0) return;
  537. for (i = 0, x = 0.f;; i++, x += 1.f) {
  538. hi = i + fixed / 2;
  539. lo = hi - fixed;
  540. if(lo>=0)break;
  541. tN = N[hi] + N[-lo];
  542. tX = X[hi] - X[-lo];
  543. tXX = XX[hi] + XX[-lo];
  544. tY = Y[hi] + Y[-lo];
  545. tXY = XY[hi] - XY[-lo];
  546. A = tY * tXX - tX * tXY;
  547. B = tN * tXY - tX * tY;
  548. D = tN * tXX - tX * tX;
  549. R = (A + x * B) / D;
  550. if (R - offset < noise[i]) noise[i] = R - offset;
  551. }
  552. for ( ;; i++, x += 1.f) {
  553. hi = i + fixed / 2;
  554. lo = hi - fixed;
  555. if(hi>=n)break;
  556. tN = N[hi] - N[lo];
  557. tX = X[hi] - X[lo];
  558. tXX = XX[hi] - XX[lo];
  559. tY = Y[hi] - Y[lo];
  560. tXY = XY[hi] - XY[lo];
  561. A = tY * tXX - tX * tXY;
  562. B = tN * tXY - tX * tY;
  563. D = tN * tXX - tX * tX;
  564. R = (A + x * B) / D;
  565. if (R - offset < noise[i]) noise[i] = R - offset;
  566. }
  567. for ( ; i < n; i++, x += 1.f) {
  568. R = (A + x * B) / D;
  569. if (R - offset < noise[i]) noise[i] = R - offset;
  570. }
  571. }
  572. void _vp_noisemask(vorbis_look_psy *p,
  573. float *logmdct,
  574. float *logmask){
  575. int i,n=p->n;
  576. float *work=alloca(n*sizeof(*work));
  577. bark_noise_hybridmp(n,p->bark,logmdct,logmask,
  578. 140.,-1);
  579. for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
  580. bark_noise_hybridmp(n,p->bark,work,logmask,0.,
  581. p->vi->noisewindowfixed);
  582. for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
  583. #if 0
  584. {
  585. static int seq=0;
  586. float work2[n];
  587. for(i=0;i<n;i++){
  588. work2[i]=logmask[i]+work[i];
  589. }
  590. if(seq&1)
  591. _analysis_output("median2R",seq/2,work,n,1,0,0);
  592. else
  593. _analysis_output("median2L",seq/2,work,n,1,0,0);
  594. if(seq&1)
  595. _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
  596. else
  597. _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
  598. seq++;
  599. }
  600. #endif
  601. for(i=0;i<n;i++){
  602. int dB=logmask[i]+.5;
  603. if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
  604. if(dB<0)dB=0;
  605. logmask[i]= work[i]+p->vi->noisecompand[dB];
  606. }
  607. }
  608. void _vp_tonemask(vorbis_look_psy *p,
  609. float *logfft,
  610. float *logmask,
  611. float global_specmax,
  612. float local_specmax){
  613. int i,n=p->n;
  614. float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
  615. float att=local_specmax+p->vi->ath_adjatt;
  616. for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
  617. /* set the ATH (floating below localmax, not global max by a
  618. specified att) */
  619. if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
  620. for(i=0;i<n;i++)
  621. logmask[i]=p->ath[i]+att;
  622. /* tone masking */
  623. seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
  624. max_seeds(p,seed,logmask);
  625. }
  626. void _vp_offset_and_mix(vorbis_look_psy *p,
  627. float *noise,
  628. float *tone,
  629. int offset_select,
  630. float *logmask,
  631. float *mdct,
  632. float *logmdct){
  633. int i,n=p->n;
  634. float de, coeffi, cx;/* AoTuV */
  635. float toneatt=p->vi->tone_masteratt[offset_select];
  636. cx = p->m_val;
  637. for(i=0;i<n;i++){
  638. float val= noise[i]+p->noiseoffset[offset_select][i];
  639. if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
  640. logmask[i]=max(val,tone[i]+toneatt);
  641. /* AoTuV */
  642. /** @ M1 **
  643. The following codes improve a noise problem.
  644. A fundamental idea uses the value of masking and carries out
  645. the relative compensation of the MDCT.
  646. However, this code is not perfect and all noise problems cannot be solved.
  647. by Aoyumi @ 2004/04/18
  648. */
  649. if(offset_select == 1) {
  650. coeffi = -17.2; /* coeffi is a -17.2dB threshold */
  651. val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
  652. if(val > coeffi){
  653. /* mdct value is > -17.2 dB below floor */
  654. de = 1.0-((val-coeffi)*0.005*cx);
  655. /* pro-rated attenuation:
  656. -0.00 dB boost if mdct value is -17.2dB (relative to floor)
  657. -0.77 dB boost if mdct value is 0dB (relative to floor)
  658. -1.64 dB boost if mdct value is +17.2dB (relative to floor)
  659. etc... */
  660. if(de < 0) de = 0.0001;
  661. }else
  662. /* mdct value is <= -17.2 dB below floor */
  663. de = 1.0-((val-coeffi)*0.0003*cx);
  664. /* pro-rated attenuation:
  665. +0.00 dB atten if mdct value is -17.2dB (relative to floor)
  666. +0.45 dB atten if mdct value is -34.4dB (relative to floor)
  667. etc... */
  668. mdct[i] *= de;
  669. }
  670. }
  671. }
  672. float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
  673. vorbis_info *vi=vd->vi;
  674. codec_setup_info *ci=vi->codec_setup;
  675. vorbis_info_psy_global *gi=&ci->psy_g_param;
  676. int n=ci->blocksizes[vd->W]/2;
  677. float secs=(float)n/vi->rate;
  678. amp+=secs*gi->ampmax_att_per_sec;
  679. if(amp<-9999)amp=-9999;
  680. return(amp);
  681. }
  682. static float FLOOR1_fromdB_LOOKUP[256]={
  683. 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
  684. 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
  685. 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
  686. 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
  687. 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
  688. 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
  689. 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
  690. 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
  691. 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
  692. 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
  693. 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
  694. 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
  695. 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
  696. 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
  697. 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
  698. 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
  699. 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
  700. 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
  701. 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
  702. 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
  703. 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
  704. 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
  705. 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
  706. 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
  707. 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
  708. 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
  709. 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
  710. 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
  711. 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
  712. 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
  713. 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
  714. 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
  715. 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
  716. 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
  717. 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
  718. 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
  719. 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
  720. 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
  721. 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
  722. 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
  723. 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
  724. 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
  725. 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
  726. 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
  727. 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
  728. 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
  729. 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
  730. 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
  731. 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
  732. 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
  733. 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
  734. 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
  735. 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
  736. 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
  737. 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
  738. 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
  739. 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
  740. 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
  741. 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
  742. 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
  743. 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
  744. 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
  745. 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
  746. 0.82788260F, 0.88168307F, 0.9389798F, 1.F,
  747. };
  748. /* this is for per-channel noise normalization */
  749. static int apsort(const void *a, const void *b){
  750. float f1=**(float**)a;
  751. float f2=**(float**)b;
  752. return (f1<f2)-(f1>f2);
  753. }
  754. static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
  755. float *floor, int *flag, int i, int jn){
  756. int j;
  757. for(j=0;j<jn;j++){
  758. float point = j>=limit-i ? postpoint : prepoint;
  759. float r = fabs(mdct[j])/floor[j];
  760. if(r<point)
  761. flag[j]=0;
  762. else
  763. flag[j]=1;
  764. }
  765. }
  766. /* Overload/Side effect: On input, the *q vector holds either the
  767. quantized energy (for elements with the flag set) or the absolute
  768. values of the *r vector (for elements with flag unset). On output,
  769. *q holds the quantized energy for all elements */
  770. static float noise_normalize(vorbis_look_psy *p, int limit, float *r, float *q, float *f, int *flags, float acc, int i, int n, int *out){
  771. vorbis_info_psy *vi=p->vi;
  772. float **sort = alloca(n*sizeof(*sort));
  773. int j,count=0;
  774. int start = (vi->normal_p ? vi->normal_start-i : n);
  775. if(start>n)start=n;
  776. /* force classic behavior where only energy in the current band is considered */
  777. acc=0.f;
  778. /* still responsible for populating *out where noise norm not in
  779. effect. There's no need to [re]populate *q in these areas */
  780. for(j=0;j<start;j++){
  781. if(!flags || !flags[j]){ /* lossless coupling already quantized.
  782. Don't touch; requantizing based on
  783. energy would be incorrect. */
  784. float ve = q[j]/f[j];
  785. if(r[j]<0)
  786. out[j] = -rint(sqrt(ve));
  787. else
  788. out[j] = rint(sqrt(ve));
  789. }
  790. }
  791. /* sort magnitudes for noise norm portion of partition */
  792. for(;j<n;j++){
  793. if(!flags || !flags[j]){ /* can't noise norm elements that have
  794. already been loslessly coupled; we can
  795. only account for their energy error */
  796. float ve = q[j]/f[j];
  797. /* Despite all the new, more capable coupling code, for now we
  798. implement noise norm as it has been up to this point. Only
  799. consider promotions to unit magnitude from 0. In addition
  800. the only energy error counted is quantizations to zero. */
  801. /* also-- the original point code only applied noise norm at > pointlimit */
  802. if(ve<.25f && (!flags || j>=limit-i)){
  803. acc += ve;
  804. sort[count++]=q+j; /* q is fabs(r) for unflagged element */
  805. }else{
  806. /* For now: no acc adjustment for nonzero quantization. populate *out and q as this value is final. */
  807. if(r[j]<0)
  808. out[j] = -rint(sqrt(ve));
  809. else
  810. out[j] = rint(sqrt(ve));
  811. q[j] = out[j]*out[j]*f[j];
  812. }
  813. }/* else{
  814. again, no energy adjustment for error in nonzero quant-- for now
  815. }*/
  816. }
  817. if(count){
  818. /* noise norm to do */
  819. qsort(sort,count,sizeof(*sort),apsort);
  820. for(j=0;j<count;j++){
  821. int k=sort[j]-q;
  822. if(acc>=vi->normal_thresh){
  823. out[k]=unitnorm(r[k]);
  824. acc-=1.f;
  825. q[k]=f[k];
  826. }else{
  827. out[k]=0;
  828. q[k]=0.f;
  829. }
  830. }
  831. }
  832. return acc;
  833. }
  834. /* Noise normalization, quantization and coupling are not wholly
  835. seperable processes in depth>1 coupling. */
  836. void _vp_couple_quantize_normalize(int blobno,
  837. vorbis_info_psy_global *g,
  838. vorbis_look_psy *p,
  839. vorbis_info_mapping0 *vi,
  840. float **mdct,
  841. int **iwork,
  842. int *nonzero,
  843. int sliding_lowpass,
  844. int ch){
  845. int i;
  846. int n = p->n;
  847. int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
  848. int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
  849. float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
  850. float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
  851. #if 0
  852. float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
  853. #endif
  854. /* mdct is our raw mdct output, floor not removed. */
  855. /* inout passes in the ifloor, passes back quantized result */
  856. /* unquantized energy (negative indicates amplitude has negative sign) */
  857. float **raw = alloca(ch*sizeof(*raw));
  858. /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
  859. float **quant = alloca(ch*sizeof(*quant));
  860. /* floor energy */
  861. float **floor = alloca(ch*sizeof(*floor));
  862. /* flags indicating raw/quantized status of elements in raw vector */
  863. int **flag = alloca(ch*sizeof(*flag));
  864. /* non-zero flag working vector */
  865. int *nz = alloca(ch*sizeof(*nz));
  866. /* energy surplus/defecit tracking */
  867. float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc));
  868. /* The threshold of a stereo is changed with the size of n */
  869. if(n > 1000)
  870. postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
  871. raw[0] = alloca(ch*partition*sizeof(**raw));
  872. quant[0] = alloca(ch*partition*sizeof(**quant));
  873. floor[0] = alloca(ch*partition*sizeof(**floor));
  874. flag[0] = alloca(ch*partition*sizeof(**flag));
  875. for(i=1;i<ch;i++){
  876. raw[i] = &raw[0][partition*i];
  877. quant[i] = &quant[0][partition*i];
  878. floor[i] = &floor[0][partition*i];
  879. flag[i] = &flag[0][partition*i];
  880. }
  881. for(i=0;i<ch+vi->coupling_steps;i++)
  882. acc[i]=0.f;
  883. for(i=0;i<n;i+=partition){
  884. int k,j,jn = partition > n-i ? n-i : partition;
  885. int step,track = 0;
  886. memcpy(nz,nonzero,sizeof(*nz)*ch);
  887. /* prefill */
  888. memset(flag[0],0,ch*partition*sizeof(**flag));
  889. for(k=0;k<ch;k++){
  890. int *iout = &iwork[k][i];
  891. if(nz[k]){
  892. for(j=0;j<jn;j++)
  893. floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
  894. flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
  895. for(j=0;j<jn;j++){
  896. quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
  897. if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
  898. floor[k][j]*=floor[k][j];
  899. }
  900. acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
  901. }else{
  902. for(j=0;j<jn;j++){
  903. floor[k][j] = 1e-10f;
  904. raw[k][j] = 0.f;
  905. quant[k][j] = 0.f;
  906. flag[k][j] = 0;
  907. iout[j]=0;
  908. }
  909. acc[track]=0.f;
  910. }
  911. track++;
  912. }
  913. /* coupling */
  914. for(step=0;step<vi->coupling_steps;step++){
  915. int Mi = vi->coupling_mag[step];
  916. int Ai = vi->coupling_ang[step];
  917. int *iM = &iwork[Mi][i];
  918. int *iA = &iwork[Ai][i];
  919. float *reM = raw[Mi];
  920. float *reA = raw[Ai];
  921. float *qeM = quant[Mi];
  922. float *qeA = quant[Ai];
  923. float *floorM = floor[Mi];
  924. float *floorA = floor[Ai];
  925. int *fM = flag[Mi];
  926. int *fA = flag[Ai];
  927. if(nz[Mi] || nz[Ai]){
  928. nz[Mi] = nz[Ai] = 1;
  929. for(j=0;j<jn;j++){
  930. if(j<sliding_lowpass-i){
  931. if(fM[j] || fA[j]){
  932. /* lossless coupling */
  933. reM[j] = fabs(reM[j])+fabs(reA[j]);
  934. qeM[j] = qeM[j]+qeA[j];
  935. fM[j]=fA[j]=1;
  936. /* couple iM/iA */
  937. {
  938. int A = iM[j];
  939. int B = iA[j];
  940. if(abs(A)>abs(B)){
  941. iA[j]=(A>0?A-B:B-A);
  942. }else{
  943. iA[j]=(B>0?A-B:B-A);
  944. iM[j]=B;
  945. }
  946. /* collapse two equivalent tuples to one */
  947. if(iA[j]>=abs(iM[j])*2){
  948. iA[j]= -iA[j];
  949. iM[j]= -iM[j];
  950. }
  951. }
  952. }else{
  953. /* lossy (point) coupling */
  954. if(j<limit-i){
  955. /* dipole */
  956. reM[j] += reA[j];
  957. qeM[j] = fabs(reM[j]);
  958. }else{
  959. #if 0
  960. /* AoTuV */
  961. /** @ M2 **
  962. The boost problem by the combination of noise normalization and point stereo is eased.
  963. However, this is a temporary patch.
  964. by Aoyumi @ 2004/04/18
  965. */
  966. float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
  967. /* elliptical */
  968. if(reM[j]+reA[j]<0){
  969. reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
  970. }else{
  971. reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
  972. }
  973. #else
  974. /* elliptical */
  975. if(reM[j]+reA[j]<0){
  976. reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
  977. }else{
  978. reM[j] = (qeM[j] = fabs(reM[j])+fabs(reA[j]));
  979. }
  980. #endif
  981. }
  982. reA[j]=qeA[j]=0.f;
  983. fA[j]=1;
  984. iA[j]=0;
  985. }
  986. }
  987. floorM[j]=floorA[j]=floorM[j]+floorA[j];
  988. }
  989. /* normalize the resulting mag vector */
  990. acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
  991. track++;
  992. }
  993. }
  994. }
  995. for(i=0;i<vi->coupling_steps;i++){
  996. /* make sure coupling a zero and a nonzero channel results in two
  997. nonzero channels. */
  998. if(nonzero[vi->coupling_mag[i]] ||
  999. nonzero[vi->coupling_ang[i]]){
  1000. nonzero[vi->coupling_mag[i]]=1;
  1001. nonzero[vi->coupling_ang[i]]=1;
  1002. }
  1003. }
  1004. }