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