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 "vorbis_codec.h"
18 #include "vorbis_codec_internal.h"
19
20 #include "vorbis_masking.h"
21 #include "vorbis_psy.h"
22 #include "vorbis_os.h"
23 #include "vorbis_lpc.h"
24 #include "vorbis_smallft.h"
25 #include "vorbis_scales.h"
26 #include "vorbis_misc.h"
27
28 #define NEGINF -9999.f
29 static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
30 static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10};
31
_vp_global_look(vorbis_info * vi)32 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
33 codec_setup_info *ci = (codec_setup_info *) vi->codec_setup;
34 vorbis_info_psy_global *gi=&ci->psy_g_param;
35 vorbis_look_psy_global *look = (vorbis_look_psy_global *)_Melder_calloc(1,sizeof(*look));
36
37 look->channels=vi->channels;
38
39 look->ampmax=-9999.;
40 look->gi=gi;
41 return(look);
42 }
43
_vp_global_free(vorbis_look_psy_global * look)44 void _vp_global_free(vorbis_look_psy_global *look){
45 if(look){
46 memset(look,0,sizeof(*look));
47 Melder_free(look);
48 }
49 }
50
_vi_gpsy_free(vorbis_info_psy_global * i)51 void _vi_gpsy_free(vorbis_info_psy_global *i){
52 if(i){
53 memset(i,0,sizeof(*i));
54 Melder_free(i);
55 }
56 }
57
_vi_psy_free(vorbis_info_psy * i)58 void _vi_psy_free(vorbis_info_psy *i){
59 if(i){
60 memset(i,0,sizeof(*i));
61 Melder_free(i);
62 }
63 }
64
min_curve(float * c,float * c2)65 static void min_curve(float *c,
66 float *c2){
67 int i;
68 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
69 }
max_curve(float * c,float * c2)70 static void max_curve(float *c,
71 float *c2){
72 int i;
73 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
74 }
75
attenuate_curve(float * c,float att)76 static void attenuate_curve(float *c,float att){
77 int i;
78 for(i=0;i<EHMER_MAX;i++)
79 c[i]+=att;
80 }
81
setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,float center_boost,float center_decay_rate)82 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
83 float center_boost, float center_decay_rate){
84 int i,j,k,m;
85 float ath[EHMER_MAX];
86 float workc[P_BANDS][P_LEVELS][EHMER_MAX];
87 float athc[P_LEVELS][EHMER_MAX];
88 float *brute_buffer = (float *) alloca(n*sizeof(*brute_buffer));
89
90 float ***ret = (float ***) _Melder_malloc(sizeof(*ret)*P_BANDS);
91
92 memset(workc,0,sizeof(workc));
93
94 for(i=0;i<P_BANDS;i++){
95 /* we add back in the ATH to avoid low level curves falling off to
96 -infinity and unnecessarily cutting off high level curves in the
97 curve limiting (last step). */
98
99 /* A half-band's settings must be valid over the whole band, and
100 it's better to mask too little than too much */
101 int ath_offset=i*4;
102 for(j=0;j<EHMER_MAX;j++){
103 float min=999.;
104 for(k=0;k<4;k++)
105 if(j+k+ath_offset<MAX_ATH){
106 if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
107 }else{
108 if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
109 }
110 ath[j]=min;
111 }
112
113 /* copy curves into working space, replicate the 50dB curve to 30
114 and 40, replicate the 100dB curve to 110 */
115 for(j=0;j<6;j++)
116 memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
117 memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
118 memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
119
120 /* apply centered curve boost/decay */
121 for(j=0;j<P_LEVELS;j++){
122 for(k=0;k<EHMER_MAX;k++){
123 float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
124 if(adj<0. && center_boost>0)adj=0.;
125 if(adj>0. && center_boost<0)adj=0.;
126 workc[i][j][k]+=adj;
127 }
128 }
129
130 /* normalize curves so the driving amplitude is 0dB */
131 /* make temp curves with the ATH overlayed */
132 for(j=0;j<P_LEVELS;j++){
133 attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0);
134 memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
135 attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0);
136 max_curve(athc[j],workc[i][j]);
137 }
138
139 /* Now limit the louder curves.
140
141 the idea is this: We don't know what the playback attenuation
142 will be; 0dB SL moves every time the user twiddles the volume
143 knob. So that means we have to use a single 'most pessimal' curve
144 for all masking amplitudes, right? Wrong. The *loudest* sound
145 can be in (we assume) a range of ...+100dB] SL. However, sounds
146 20dB down will be in a range ...+80], 40dB down is from ...+60],
147 etc... */
148
149 for(j=1;j<P_LEVELS;j++){
150 min_curve(athc[j],athc[j-1]);
151 min_curve(workc[i][j],athc[j]);
152 }
153 }
154
155 for(i=0;i<P_BANDS;i++){
156 int hi_curve,lo_curve,bin;
157 ret[i] = (float **) _Melder_malloc(sizeof(**ret)*P_LEVELS);
158
159 /* low frequency curves are measured with greater resolution than
160 the MDCT/FFT will actually give us; we want the curve applied
161 to the tone data to be pessimistic and thus apply the minimum
162 masking possible for a given bin. That means that a single bin
163 could span more than one octave and that the curve will be a
164 composite of multiple octaves. It also may mean that a single
165 bin may span > an eighth of an octave and that the eighth
166 octave values may also be composited. */
167
168 /* which octave curves will we be compositing? */
169 bin=floor(fromOC(i*.5)/binHz);
170 lo_curve= ceil(toOC(bin*binHz+1)*2);
171 hi_curve= floor(toOC((bin+1)*binHz)*2);
172 if(lo_curve>i)lo_curve=i;
173 if(lo_curve<0)lo_curve=0;
174 if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;
175
176 for(m=0;m<P_LEVELS;m++){
177 ret[i][m] = (float *) _Melder_malloc(sizeof(***ret)*(EHMER_MAX+2));
178
179 for(j=0;j<n;j++)brute_buffer[j]=999.;
180
181 /* render the curve into bins, then pull values back into curve.
182 The point is that any inherent subsampling aliasing results in
183 a safe minimum */
184 for(k=lo_curve;k<=hi_curve;k++){
185 int l=0;
186
187 for(j=0;j<EHMER_MAX;j++){
188 int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz;
189 int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1;
190
191 if(lo_bin<0)lo_bin=0;
192 if(lo_bin>n)lo_bin=n;
193 if(lo_bin<l)l=lo_bin;
194 if(hi_bin<0)hi_bin=0;
195 if(hi_bin>n)hi_bin=n;
196
197 for(;l<hi_bin && l<n;l++)
198 if(brute_buffer[l]>workc[k][m][j])
199 brute_buffer[l]=workc[k][m][j];
200 }
201
202 for(;l<n;l++)
203 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
204 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
205
206 }
207
208 /* be equally paranoid about being valid up to next half ocatve */
209 if(i+1<P_BANDS){
210 int l=0;
211 k=i+1;
212 for(j=0;j<EHMER_MAX;j++){
213 int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz;
214 int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1;
215
216 if(lo_bin<0)lo_bin=0;
217 if(lo_bin>n)lo_bin=n;
218 if(lo_bin<l)l=lo_bin;
219 if(hi_bin<0)hi_bin=0;
220 if(hi_bin>n)hi_bin=n;
221
222 for(;l<hi_bin && l<n;l++)
223 if(brute_buffer[l]>workc[k][m][j])
224 brute_buffer[l]=workc[k][m][j];
225 }
226
227 for(;l<n;l++)
228 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
229 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
230
231 }
232
233
234 for(j=0;j<EHMER_MAX;j++){
235 int binn=fromOC(j*.125+i*.5-2.)/binHz;
236 if(binn<0){
237 ret[i][m][j+2]=-999.;
238 }else{
239 if(binn>=n){
240 ret[i][m][j+2]=-999.;
241 }else{
242 ret[i][m][j+2]=brute_buffer[binn];
243 }
244 }
245 }
246
247 /* add fenceposts */
248 for(j=0;j<EHMER_OFFSET;j++)
249 if(ret[i][m][j+2]>-200.f)break;
250 ret[i][m][0]=j;
251
252 for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
253 if(ret[i][m][j+2]>-200.f)
254 break;
255 ret[i][m][1]=j;
256
257 }
258 }
259
260 return(ret);
261 }
262
_vp_psy_init(vorbis_look_psy * p,vorbis_info_psy * vi,vorbis_info_psy_global * gi,int n,long rate)263 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
264 vorbis_info_psy_global *gi,int n,long rate){
265 long i,j,lo=-99,hi=1;
266 long maxoc;
267 memset(p,0,sizeof(*p));
268
269 p->eighth_octave_lines=gi->eighth_octave_lines;
270 p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
271
272 p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
273 maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f;
274 p->total_octave_lines=maxoc-p->firstoc+1;
275 p->ath= (float *) _Melder_malloc(n*sizeof(*p->ath));
276
277 p->octave = (long *)_Melder_malloc(n*sizeof(*p->octave));
278 p->bark = (long *) _Melder_malloc(n*sizeof(*p->bark));
279 p->vi=vi;
280 p->n=n;
281 p->rate=rate;
282
283 /* AoTuV HF weighting */
284 p->m_val = 1.;
285 if(rate < 26000) p->m_val = 0;
286 else if(rate < 38000) p->m_val = .94; /* 32kHz */
287 else if(rate > 46000) p->m_val = 1.275; /* 48kHz */
288
289 /* set up the lookups for a given blocksize and sample rate */
290
291 for(i=0,j=0;i<MAX_ATH-1;i++){
292 int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
293 float base=ATH[i];
294 if(j<endpos){
295 float delta=(ATH[i+1]-base)/(endpos-j);
296 for(;j<endpos && j<n;j++){
297 p->ath[j]=base+100.;
298 base+=delta;
299 }
300 }
301 }
302
303 for(;j<n;j++){
304 p->ath[j]=p->ath[j-1];
305 }
306
307 for(i=0;i<n;i++){
308 float bark=toBARK(rate/(2*n)*i);
309
310 for(;lo+vi->noisewindowlomin<i &&
311 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
312
313 for(;hi<=n && (hi<i+vi->noisewindowhimin ||
314 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
315
316 p->bark[i]=((lo-1)<<16)+(hi-1);
317
318 }
319
320 for(i=0;i<n;i++)
321 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
322
323 p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
324 vi->tone_centerboost,vi->tone_decay);
325
326 /* set up rolling noise median */
327 p->noiseoffset = (float **) _Melder_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
328 for(i=0;i<P_NOISECURVES;i++)
329 p->noiseoffset[i]= (float *) _Melder_malloc(n*sizeof(**p->noiseoffset));
330
331 for(i=0;i<n;i++){
332 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
333 int inthalfoc;
334 float del;
335
336 if(halfoc<0)halfoc=0;
337 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
338 inthalfoc=(int)halfoc;
339 del=halfoc-inthalfoc;
340
341 for(j=0;j<P_NOISECURVES;j++)
342 p->noiseoffset[j][i]=
343 p->vi->noiseoff[j][inthalfoc]*(1.-del) +
344 p->vi->noiseoff[j][inthalfoc+1]*del;
345
346 }
347 #if 0
348 {
349 static int ls=0;
350 _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
351 _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
352 _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
353 }
354 #endif
355 }
356
_vp_psy_clear(vorbis_look_psy * p)357 void _vp_psy_clear(vorbis_look_psy *p){
358 int i,j;
359 if(p){
360 if(p->ath)Melder_free(p->ath);
361 if(p->octave)Melder_free(p->octave);
362 if(p->bark)Melder_free(p->bark);
363 if(p->tonecurves){
364 for(i=0;i<P_BANDS;i++){
365 for(j=0;j<P_LEVELS;j++){
366 Melder_free(p->tonecurves[i][j]);
367 }
368 Melder_free(p->tonecurves[i]);
369 }
370 Melder_free(p->tonecurves);
371 }
372 if(p->noiseoffset){
373 for(i=0;i<P_NOISECURVES;i++){
374 Melder_free(p->noiseoffset[i]);
375 }
376 Melder_free(p->noiseoffset);
377 }
378 memset(p,0,sizeof(*p));
379 }
380 }
381
382 /* 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)383 static void seed_curve(float *seed,
384 const float **curves,
385 float amp,
386 int oc, int n,
387 int linesper,float dBoffset){
388 int i,post1;
389 int seedptr;
390 const float *posts,*curve;
391
392 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
393 choice=std::max(choice,0);
394 choice=std::min(choice,P_LEVELS-1);
395 posts=curves[choice];
396 curve=posts+2;
397 post1=(int)posts[1];
398 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
399
400 for(i=posts[0];i<post1;i++){
401 if(seedptr>0){
402 float lin=amp+curve[i];
403 if(seed[seedptr]<lin)seed[seedptr]=lin;
404 }
405 seedptr+=linesper;
406 if(seedptr>=n)break;
407 }
408 }
409
seed_loop(vorbis_look_psy * p,const float *** curves,const float * f,const float * flr,float * seed,float specmax)410 static void seed_loop(vorbis_look_psy *p,
411 const float ***curves,
412 const float *f,
413 const float *flr,
414 float *seed,
415 float specmax){
416 vorbis_info_psy *vi=p->vi;
417 long n=p->n,i;
418 float dBoffset=vi->max_curve_dB-specmax;
419
420 /* prime the working vector with peak values */
421
422 for(i=0;i<n;i++){
423 float max=f[i];
424 long oc=p->octave[i];
425 while(i+1<n && p->octave[i+1]==oc){
426 i++;
427 if(f[i]>max)max=f[i];
428 }
429
430 if(max+6.f>flr[i]){
431 oc=oc>>p->shiftoc;
432
433 if(oc>=P_BANDS)oc=P_BANDS-1;
434 if(oc<0)oc=0;
435
436 seed_curve(seed,
437 curves[oc],
438 max,
439 p->octave[i]-p->firstoc,
440 p->total_octave_lines,
441 p->eighth_octave_lines,
442 dBoffset);
443 }
444 }
445 }
446
seed_chase(float * seeds,int linesper,long n)447 static void seed_chase(float *seeds, int linesper, long n){
448 long *posstack = (long *) alloca(n*sizeof(*posstack));
449 float *ampstack = (float *) alloca(n*sizeof(*ampstack));
450 long stack=0;
451 long pos=0;
452 long i;
453
454 for(i=0;i<n;i++){
455 if(stack<2){
456 posstack[stack]=i;
457 ampstack[stack++]=seeds[i];
458 }else{
459 while(1){
460 if(seeds[i]<ampstack[stack-1]){
461 posstack[stack]=i;
462 ampstack[stack++]=seeds[i];
463 break;
464 }else{
465 if(i<posstack[stack-1]+linesper){
466 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
467 i<posstack[stack-2]+linesper){
468 /* we completely overlap, making stack-1 irrelevant. pop it */
469 stack--;
470 continue;
471 }
472 }
473 posstack[stack]=i;
474 ampstack[stack++]=seeds[i];
475 break;
476
477 }
478 }
479 }
480 }
481
482 /* the stack now contains only the positions that are relevant. Scan
483 'em straight through */
484
485 for(i=0;i<stack;i++){
486 long endpos;
487 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
488 endpos=posstack[i+1];
489 }else{
490 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
491 discarded in short frames */
492 }
493 if(endpos>n)endpos=n;
494 for(;pos<endpos;pos++)
495 seeds[pos]=ampstack[i];
496 }
497
498 /* there. Linear time. I now remember this was on a problem set I
499 had in Grad Skool... I didn't solve it at the time ;-) */
500
501 }
502
503 /* bleaugh, this is more complicated than it needs to be */
504 #include<stdio.h>
max_seeds(vorbis_look_psy * p,float * seed,float * flr)505 static void max_seeds(vorbis_look_psy *p,
506 float *seed,
507 float *flr){
508 long n=p->total_octave_lines;
509 int linesper=p->eighth_octave_lines;
510 long linpos=0;
511 long pos;
512
513 seed_chase(seed,linesper,n); /* for masking */
514
515 pos=p->octave[0]-p->firstoc-(linesper>>1);
516
517 while(linpos+1<p->n){
518 float minV=seed[pos];
519 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
520 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
521 while(pos+1<=end){
522 pos++;
523 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
524 minV=seed[pos];
525 }
526
527 end=pos+p->firstoc;
528 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
529 if(flr[linpos]<minV)flr[linpos]=minV;
530 }
531
532 {
533 float minV=seed[p->total_octave_lines-1];
534 for(;linpos<p->n;linpos++)
535 if(flr[linpos]<minV)flr[linpos]=minV;
536 }
537
538 }
539
bark_noise_hybridmp(int n,const long * b,const float * f,float * noise,const float offset,const int fixed)540 static void bark_noise_hybridmp(int n,const long *b,
541 const float *f,
542 float *noise,
543 const float offset,
544 const int fixed){
545
546 float *N = (float *) alloca(n*sizeof(*N));
547 float *X = (float *) alloca(n*sizeof(*N));
548 float *XX = (float *) alloca(n*sizeof(*N));
549 float *Y = (float *) alloca(n*sizeof(*N));
550 float *XY = (float *) alloca(n*sizeof(*N));
551
552 float tN, tX, tXX, tY, tXY;
553 int i;
554
555 int lo, hi;
556 float R=0.f;
557 float A=0.f;
558 float B=0.f;
559 float D=1.f;
560 float w, x, y;
561
562 tN = tX = tXX = tY = tXY = 0.f;
563
564 y = f[0] + offset;
565 if (y < 1.f) y = 1.f;
566
567 w = y * y * .5;
568
569 tN += w;
570 tX += w;
571 tY += w * y;
572
573 N[0] = tN;
574 X[0] = tX;
575 XX[0] = tXX;
576 Y[0] = tY;
577 XY[0] = tXY;
578
579 for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
580
581 y = f[i] + offset;
582 if (y < 1.f) y = 1.f;
583
584 w = y * y;
585
586 tN += w;
587 tX += w * x;
588 tXX += w * x * x;
589 tY += w * y;
590 tXY += w * x * y;
591
592 N[i] = tN;
593 X[i] = tX;
594 XX[i] = tXX;
595 Y[i] = tY;
596 XY[i] = tXY;
597 }
598
599 for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
600
601 lo = b[i] >> 16;
602 hi = b[i] & 0xffff;
603 if( lo>=0 || -lo>=n ) break;
604 if( hi>=n ) break;
605
606 tN = N[hi] + N[-lo];
607 tX = X[hi] - X[-lo];
608 tXX = XX[hi] + XX[-lo];
609 tY = Y[hi] + Y[-lo];
610 tXY = XY[hi] - XY[-lo];
611
612 A = tY * tXX - tX * tXY;
613 B = tN * tXY - tX * tY;
614 D = tN * tXX - tX * tX;
615 R = (A + x * B) / D;
616 if (R < 0.f) R = 0.f;
617
618 noise[i] = R - offset;
619 }
620
621 for ( ; i < n; i++, x += 1.f) {
622
623 lo = b[i] >> 16;
624 hi = b[i] & 0xffff;
625 if( lo<0 || lo>=n ) break;
626 if( hi>=n ) break;
627
628 tN = N[hi] - N[lo];
629 tX = X[hi] - X[lo];
630 tXX = XX[hi] - XX[lo];
631 tY = Y[hi] - Y[lo];
632 tXY = XY[hi] - XY[lo];
633
634 A = tY * tXX - tX * tXY;
635 B = tN * tXY - tX * tY;
636 D = tN * tXX - tX * tX;
637 R = (A + x * B) / D;
638 if (R < 0.f) R = 0.f;
639
640 noise[i] = R - offset;
641 }
642
643 for ( ; i < n; i++, x += 1.f) {
644
645 R = (A + x * B) / D;
646 if (R < 0.f) R = 0.f;
647
648 noise[i] = R - offset;
649 }
650
651 if (fixed <= 0) return;
652
653 for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
654 hi = i + fixed / 2;
655 lo = hi - fixed;
656 if ( hi>=n ) break;
657 if ( lo>=0 ) break;
658
659 tN = N[hi] + N[-lo];
660 tX = X[hi] - X[-lo];
661 tXX = XX[hi] + XX[-lo];
662 tY = Y[hi] + Y[-lo];
663 tXY = XY[hi] - XY[-lo];
664
665
666 A = tY * tXX - tX * tXY;
667 B = tN * tXY - tX * tY;
668 D = tN * tXX - tX * tX;
669 R = (A + x * B) / D;
670
671 if (R - offset < noise[i]) noise[i] = R - offset;
672 }
673 for ( ; i < n; i++, x += 1.f) {
674
675 hi = i + fixed / 2;
676 lo = hi - fixed;
677 if ( hi>=n ) break;
678 if ( lo<0 ) break;
679
680 tN = N[hi] - N[lo];
681 tX = X[hi] - X[lo];
682 tXX = XX[hi] - XX[lo];
683 tY = Y[hi] - Y[lo];
684 tXY = XY[hi] - XY[lo];
685
686 A = tY * tXX - tX * tXY;
687 B = tN * tXY - tX * tY;
688 D = tN * tXX - tX * tX;
689 R = (A + x * B) / D;
690
691 if (R - offset < noise[i]) noise[i] = R - offset;
692 }
693 for ( ; i < n; i++, x += 1.f) {
694 R = (A + x * B) / D;
695 if (R - offset < noise[i]) noise[i] = R - offset;
696 }
697 }
698
_vp_noisemask(vorbis_look_psy * p,float * logmdct,float * logmask)699 void _vp_noisemask(vorbis_look_psy *p,
700 float *logmdct,
701 float *logmask){
702
703 int i,n=p->n;
704 float *work = (float *) alloca(n*sizeof(*work));
705
706 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
707 140.,-1);
708
709 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
710
711 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
712 p->vi->noisewindowfixed);
713
714 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 for(i=0;i<n;i++){
739 int dB=logmask[i]+.5;
740 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
741 if(dB<0)dB=0;
742 logmask[i]= work[i]+p->vi->noisecompand[dB];
743 }
744
745 }
746
_vp_tonemask(vorbis_look_psy * p,float * logfft,float * logmask,float global_specmax,float local_specmax)747 void _vp_tonemask(vorbis_look_psy *p,
748 float *logfft,
749 float *logmask,
750 float global_specmax,
751 float local_specmax){
752
753 int i,n=p->n;
754
755 float *seed = (float *) alloca(sizeof(*seed)*p->total_octave_lines);
756 float att=local_specmax+p->vi->ath_adjatt;
757 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 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
762
763 for(i=0;i<n;i++)
764 logmask[i]=p->ath[i]+att;
765
766 /* tone masking */
767 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
768 max_seeds(p,seed,logmask);
769
770 }
771
_vp_offset_and_mix(vorbis_look_psy * p,float * noise,float * tone,int offset_select,float * logmask,float * mdct,float * logmdct)772 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 int i,n=p->n;
780 float de, coeffi, cx;/* AoTuV */
781 float toneatt=p->vi->tone_masteratt[offset_select];
782
783 cx = p->m_val;
784
785 for(i=0;i<n;i++){
786 float val= noise[i]+p->noiseoffset[offset_select][i];
787 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
788 logmask[i] = std::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 if(offset_select == 1) {
801 coeffi = -17.2; /* coeffi is a -17.2dB threshold */
802 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
803
804 if(val > coeffi){
805 /* mdct value is > -17.2 dB below floor */
806
807 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 if(de < 0) de = 0.0001;
815 }else
816 /* mdct value is <= -17.2 dB below floor */
817
818 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 mdct[i] *= de;
825
826 }
827 }
828 }
829
_vp_ampmax_decay(float amp,vorbis_dsp_state * vd)830 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
831 vorbis_info *vi=vd->vi;
832 codec_setup_info *ci = (codec_setup_info *) vi->codec_setup;
833 vorbis_info_psy_global *gi=&ci->psy_g_param;
834
835 int n=ci->blocksizes[vd->W]/2;
836 float secs=(float)n/vi->rate;
837
838 amp+=secs*gi->ampmax_att_per_sec;
839 if(amp<-9999)amp=-9999;
840 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 */
apsort(const void * a,const void * b)911 static int apsort(const void *a, const void *b){
912 float f1=**(float**)a;
913 float f2=**(float**)b;
914 return (f1<f2)-(f1>f2);
915 }
916
flag_lossless(int limit,float prepoint,float postpoint,float * mdct,float * floor,int * flag,int i,int jn)917 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 for(j=0;j<jn;j++){
921 float point = j>=limit-i ? postpoint : prepoint;
922 float r = fabs(mdct[j])/floor[j];
923 if(r<point)
924 flag[j]=0;
925 else
926 flag[j]=1;
927 }
928 }
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 */
noise_normalize(vorbis_look_psy * p,int limit,float * r,float * q,float * f,int * flags,float acc,int i,int n,int * out)934 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 vorbis_info_psy *vi=p->vi;
937 float **sort = (float **) alloca(n*sizeof(*sort));
938 int j,count=0;
939 int start = (vi->normal_p ? vi->normal_start-i : n);
940 if(start>n)start=n;
941
942 /* force classic behavior where only energy in the current band is considered */
943 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 for(j=0;j<start;j++){
948 if(!flags || !flags[j]){ /* lossless coupling already quantized.
949 Don't touch; requantizing based on
950 energy would be incorrect. */
951 float ve = q[j]/f[j];
952 if(r[j]<0)
953 out[j] = -rint(sqrt(ve));
954 else
955 out[j] = rint(sqrt(ve));
956 }
957 }
958
959 /* sort magnitudes for noise norm portion of partition */
960 for(;j<n;j++){
961 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 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 if(ve<.25f && (!flags || j>=limit-i)){
971 acc += ve;
972 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 if(r[j]<0)
976 out[j] = -rint(sqrt(ve));
977 else
978 out[j] = rint(sqrt(ve));
979 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 if(count){
987 /* noise norm to do */
988 qsort(sort,count,sizeof(*sort),apsort);
989 for(j=0;j<count;j++){
990 int k=sort[j]-q;
991 if(acc>=vi->normal_thresh){
992 out[k]=unitnorm(r[k]);
993 acc-=1.f;
994 q[k]=f[k];
995 }else{
996 out[k]=0;
997 q[k]=0.f;
998 }
999 }
1000 }
1001
1002 return acc;
1003 }
1004
1005 /* Noise normalization, quantization and coupling are not wholly
1006 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)1007 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 int n = p->n;
1019 int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1020 int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1021 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1022 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 float **raw = (float **) alloca(ch*sizeof(*raw));
1032
1033 /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1034 float **quant = (float **) alloca(ch*sizeof(*quant));
1035
1036 /* floor energy */
1037 float **floor = (float **) alloca(ch*sizeof(*floor));
1038
1039 /* flags indicating raw/quantized status of elements in raw vector */
1040 int **flag = (int **) alloca(ch*sizeof(*flag));
1041
1042 /* non-zero flag working vector */
1043 int *nz = (int *) alloca(ch*sizeof(*nz));
1044
1045 /* energy surplus/defecit tracking */
1046 float *acc = (float *) alloca((ch+vi->coupling_steps)*sizeof(*acc));
1047
1048 /* The threshold of a stereo is changed with the size of n */
1049 if(n > 1000)
1050 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1051
1052 raw[0] = (float *) alloca(ch*partition*sizeof(**raw));
1053 quant[0] = (float *) alloca(ch*partition*sizeof(**quant));
1054 floor[0] = (float *) alloca(ch*partition*sizeof(**floor));
1055 flag[0] = (int *) alloca(ch*partition*sizeof(**flag));
1056
1057 for(i=1;i<ch;i++){
1058 raw[i] = &raw[0][partition*i];
1059 quant[i] = &quant[0][partition*i];
1060 floor[i] = &floor[0][partition*i];
1061 flag[i] = &flag[0][partition*i];
1062 }
1063 for(i=0;i<ch+vi->coupling_steps;i++)
1064 acc[i]=0.f;
1065
1066 for(i=0;i<n;i+=partition){
1067 int k,j,jn = partition > n-i ? n-i : partition;
1068 int step,track = 0;
1069
1070 memcpy(nz,nonzero,sizeof(*nz)*ch);
1071
1072 /* prefill */
1073 memset(flag[0],0,ch*partition*sizeof(**flag));
1074 for(k=0;k<ch;k++){
1075 int *iout = &iwork[k][i];
1076 if(nz[k]){
1077
1078 for(j=0;j<jn;j++)
1079 floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1080
1081 flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1082
1083 for(j=0;j<jn;j++){
1084 quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1085 if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1086 floor[k][j]*=floor[k][j];
1087 }
1088
1089 acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],nullptr,acc[track],i,jn,iout);
1090
1091 }else{
1092 for(j=0;j<jn;j++){
1093 floor[k][j] = 1e-10f;
1094 raw[k][j] = 0.f;
1095 quant[k][j] = 0.f;
1096 flag[k][j] = 0;
1097 iout[j]=0;
1098 }
1099 acc[track]=0.f;
1100 }
1101 track++;
1102 }
1103
1104 /* coupling */
1105 for(step=0;step<vi->coupling_steps;step++){
1106 int Mi = vi->coupling_mag[step];
1107 int Ai = vi->coupling_ang[step];
1108 int *iM = &iwork[Mi][i];
1109 int *iA = &iwork[Ai][i];
1110 float *reM = raw[Mi];
1111 float *reA = raw[Ai];
1112 float *qeM = quant[Mi];
1113 float *qeA = quant[Ai];
1114 float *floorM = floor[Mi];
1115 float *floorA = floor[Ai];
1116 int *fM = flag[Mi];
1117 int *fA = flag[Ai];
1118
1119 if(nz[Mi] || nz[Ai]){
1120 nz[Mi] = nz[Ai] = 1;
1121
1122 for(j=0;j<jn;j++){
1123
1124 if(j<sliding_lowpass-i){
1125 if(fM[j] || fA[j]){
1126 /* lossless coupling */
1127
1128 reM[j] = fabs(reM[j])+fabs(reA[j]);
1129 qeM[j] = qeM[j]+qeA[j];
1130 fM[j]=fA[j]=1;
1131
1132 /* couple iM/iA */
1133 {
1134 int A = iM[j];
1135 int B = iA[j];
1136
1137 if(abs(A)>abs(B)){
1138 iA[j]=(A>0?A-B:B-A);
1139 }else{
1140 iA[j]=(B>0?A-B:B-A);
1141 iM[j]=B;
1142 }
1143
1144 /* collapse two equivalent tuples to one */
1145 if(iA[j]>=abs(iM[j])*2){
1146 iA[j]= -iA[j];
1147 iM[j]= -iM[j];
1148 }
1149
1150 }
1151
1152 }else{
1153 /* lossy (point) coupling */
1154 if(j<limit-i){
1155 /* dipole */
1156 reM[j] += reA[j];
1157 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 if(reM[j]+reA[j]<0){
1176 reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1177 }else{
1178 reM[j] = (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1179 }
1180 #endif
1181
1182 }
1183 reA[j]=qeA[j]=0.f;
1184 fA[j]=1;
1185 iA[j]=0;
1186 }
1187 }
1188 floorM[j]=floorA[j]=floorM[j]+floorA[j];
1189 }
1190 /* normalize the resulting mag vector */
1191 acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1192 track++;
1193 }
1194 }
1195 }
1196
1197 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 if(nonzero[vi->coupling_mag[i]] ||
1201 nonzero[vi->coupling_ang[i]]){
1202 nonzero[vi->coupling_mag[i]]=1;
1203 nonzero[vi->coupling_ang[i]]=1;
1204 }
1205 }
1206 }
1207