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