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: psy.c 17077 2010-03-26 06:22:19Z xiphmont $ 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 _vp_global_look(vorbis_info * vi)36 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){ 37 codec_setup_info *ci=vi->codec_setup; 38 vorbis_info_psy_global *gi=&ci->psy_g_param; 39 vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look)); 40 41 look->channels=vi->channels; 42 43 look->ampmax=-9999.; 44 look->gi=gi; 45 return(look); 46 } 47 _vp_global_free(vorbis_look_psy_global * look)48 void _vp_global_free(vorbis_look_psy_global *look){ 49 if(look){ 50 memset(look,0,sizeof(*look)); 51 _ogg_free(look); 52 } 53 } 54 _vi_gpsy_free(vorbis_info_psy_global * i)55 void _vi_gpsy_free(vorbis_info_psy_global *i){ 56 if(i){ 57 memset(i,0,sizeof(*i)); 58 _ogg_free(i); 59 } 60 } 61 _vi_psy_free(vorbis_info_psy * i)62 void _vi_psy_free(vorbis_info_psy *i){ 63 if(i){ 64 memset(i,0,sizeof(*i)); 65 _ogg_free(i); 66 } 67 } 68 min_curve(float * c,float * c2)69 static void min_curve(float *c, 70 float *c2){ 71 int i; 72 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i]; 73 } max_curve(float * c,float * c2)74 static void max_curve(float *c, 75 float *c2){ 76 int i; 77 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i]; 78 } 79 attenuate_curve(float * c,float att)80 static void attenuate_curve(float *c,float att){ 81 int i; 82 for(i=0;i<EHMER_MAX;i++) 83 c[i]+=att; 84 } 85 setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,float center_boost,float center_decay_rate)86 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 float *brute_buffer=alloca(n*sizeof(*brute_buffer)); 93 94 float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS); 95 96 memset(workc,0,sizeof(workc)); 97 98 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 int ath_offset=i*4; 106 for(j=0;j<EHMER_MAX;j++){ 107 float min=999.; 108 for(k=0;k<4;k++) 109 if(j+k+ath_offset<MAX_ATH){ 110 if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset]; 111 }else{ 112 if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1]; 113 } 114 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 for(j=0;j<6;j++) 120 memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j])); 121 memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); 122 memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); 123 124 /* apply centered curve boost/decay */ 125 for(j=0;j<P_LEVELS;j++){ 126 for(k=0;k<EHMER_MAX;k++){ 127 float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate; 128 if(adj<0. && center_boost>0)adj=0.; 129 if(adj>0. && center_boost<0)adj=0.; 130 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 for(j=0;j<P_LEVELS;j++){ 137 attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0); 138 memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc)); 139 attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0); 140 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 for(j=1;j<P_LEVELS;j++){ 154 min_curve(athc[j],athc[j-1]); 155 min_curve(workc[i][j],athc[j]); 156 } 157 } 158 159 for(i=0;i<P_BANDS;i++){ 160 int hi_curve,lo_curve,bin; 161 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 bin=floor(fromOC(i*.5)/binHz); 174 lo_curve= ceil(toOC(bin*binHz+1)*2); 175 hi_curve= floor(toOC((bin+1)*binHz)*2); 176 if(lo_curve>i)lo_curve=i; 177 if(lo_curve<0)lo_curve=0; 178 if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1; 179 180 for(m=0;m<P_LEVELS;m++){ 181 ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2)); 182 183 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 for(k=lo_curve;k<=hi_curve;k++){ 189 int l=0; 190 191 for(j=0;j<EHMER_MAX;j++){ 192 int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz; 193 int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1; 194 195 if(lo_bin<0)lo_bin=0; 196 if(lo_bin>n)lo_bin=n; 197 if(lo_bin<l)l=lo_bin; 198 if(hi_bin<0)hi_bin=0; 199 if(hi_bin>n)hi_bin=n; 200 201 for(;l<hi_bin && l<n;l++) 202 if(brute_buffer[l]>workc[k][m][j]) 203 brute_buffer[l]=workc[k][m][j]; 204 } 205 206 for(;l<n;l++) 207 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) 208 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 if(i+1<P_BANDS){ 214 int l=0; 215 k=i+1; 216 for(j=0;j<EHMER_MAX;j++){ 217 int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz; 218 int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1; 219 220 if(lo_bin<0)lo_bin=0; 221 if(lo_bin>n)lo_bin=n; 222 if(lo_bin<l)l=lo_bin; 223 if(hi_bin<0)hi_bin=0; 224 if(hi_bin>n)hi_bin=n; 225 226 for(;l<hi_bin && l<n;l++) 227 if(brute_buffer[l]>workc[k][m][j]) 228 brute_buffer[l]=workc[k][m][j]; 229 } 230 231 for(;l<n;l++) 232 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) 233 brute_buffer[l]=workc[k][m][EHMER_MAX-1]; 234 235 } 236 237 238 for(j=0;j<EHMER_MAX;j++){ 239 int bin=fromOC(j*.125+i*.5-2.)/binHz; 240 if(bin<0){ 241 ret[i][m][j+2]=-999.; 242 }else{ 243 if(bin>=n){ 244 ret[i][m][j+2]=-999.; 245 }else{ 246 ret[i][m][j+2]=brute_buffer[bin]; 247 } 248 } 249 } 250 251 /* add fenceposts */ 252 for(j=0;j<EHMER_OFFSET;j++) 253 if(ret[i][m][j+2]>-200.f)break; 254 ret[i][m][0]=j; 255 256 for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--) 257 if(ret[i][m][j+2]>-200.f) 258 break; 259 ret[i][m][1]=j; 260 261 } 262 } 263 264 return(ret); 265 } 266 _vp_psy_init(vorbis_look_psy * p,vorbis_info_psy * vi,vorbis_info_psy_global * gi,int n,long rate)267 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi, 268 vorbis_info_psy_global *gi,int n,long rate){ 269 long i,j,lo=-99,hi=1; 270 long maxoc; 271 memset(p,0,sizeof(*p)); 272 273 p->eighth_octave_lines=gi->eighth_octave_lines; 274 p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1; 275 276 p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines; 277 maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f; 278 p->total_octave_lines=maxoc-p->firstoc+1; 279 p->ath=_ogg_malloc(n*sizeof(*p->ath)); 280 281 p->octave=_ogg_malloc(n*sizeof(*p->octave)); 282 p->bark=_ogg_malloc(n*sizeof(*p->bark)); 283 p->vi=vi; 284 p->n=n; 285 p->rate=rate; 286 287 /* AoTuV HF weighting */ 288 p->m_val = 1.; 289 if(rate < 26000) p->m_val = 0; 290 else if(rate < 38000) p->m_val = .94; /* 32kHz */ 291 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 for(i=0,j=0;i<MAX_ATH-1;i++){ 296 int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate); 297 float base=ATH[i]; 298 if(j<endpos){ 299 float delta=(ATH[i+1]-base)/(endpos-j); 300 for(;j<endpos && j<n;j++){ 301 p->ath[j]=base+100.; 302 base+=delta; 303 } 304 } 305 } 306 307 for(;j<n;j++){ 308 p->ath[j]=p->ath[j-1]; 309 } 310 311 for(i=0;i<n;i++){ 312 float bark=toBARK(rate/(2*n)*i); 313 314 for(;lo+vi->noisewindowlomin<i && 315 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++); 316 317 for(;hi<=n && (hi<i+vi->noisewindowhimin || 318 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++); 319 320 p->bark[i]=((lo-1)<<16)+(hi-1); 321 322 } 323 324 for(i=0;i<n;i++) 325 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f; 326 327 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 p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset)); 332 for(i=0;i<P_NOISECURVES;i++) 333 p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset)); 334 335 for(i=0;i<n;i++){ 336 float halfoc=toOC((i+.5)*rate/(2.*n))*2.; 337 int inthalfoc; 338 float del; 339 340 if(halfoc<0)halfoc=0; 341 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1; 342 inthalfoc=(int)halfoc; 343 del=halfoc-inthalfoc; 344 345 for(j=0;j<P_NOISECURVES;j++) 346 p->noiseoffset[j][i]= 347 p->vi->noiseoff[j][inthalfoc]*(1.-del) + 348 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 } 360 _vp_psy_clear(vorbis_look_psy * p)361 void _vp_psy_clear(vorbis_look_psy *p){ 362 int i,j; 363 if(p){ 364 if(p->ath)_ogg_free(p->ath); 365 if(p->octave)_ogg_free(p->octave); 366 if(p->bark)_ogg_free(p->bark); 367 if(p->tonecurves){ 368 for(i=0;i<P_BANDS;i++){ 369 for(j=0;j<P_LEVELS;j++){ 370 _ogg_free(p->tonecurves[i][j]); 371 } 372 _ogg_free(p->tonecurves[i]); 373 } 374 _ogg_free(p->tonecurves); 375 } 376 if(p->noiseoffset){ 377 for(i=0;i<P_NOISECURVES;i++){ 378 _ogg_free(p->noiseoffset[i]); 379 } 380 _ogg_free(p->noiseoffset); 381 } 382 memset(p,0,sizeof(*p)); 383 } 384 } 385 386 /* 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)387 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 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f); 397 choice=max(choice,0); 398 choice=min(choice,P_LEVELS-1); 399 posts=curves[choice]; 400 curve=posts+2; 401 post1=(int)posts[1]; 402 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1); 403 404 for(i=posts[0];i<post1;i++){ 405 if(seedptr>0){ 406 float lin=amp+curve[i]; 407 if(seed[seedptr]<lin)seed[seedptr]=lin; 408 } 409 seedptr+=linesper; 410 if(seedptr>=n)break; 411 } 412 } 413 seed_loop(vorbis_look_psy * p,const float *** curves,const float * f,const float * flr,float * seed,float specmax)414 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 vorbis_info_psy *vi=p->vi; 421 long n=p->n,i; 422 float dBoffset=vi->max_curve_dB-specmax; 423 424 /* prime the working vector with peak values */ 425 426 for(i=0;i<n;i++){ 427 float max=f[i]; 428 long oc=p->octave[i]; 429 while(i+1<n && p->octave[i+1]==oc){ 430 i++; 431 if(f[i]>max)max=f[i]; 432 } 433 434 if(max+6.f>flr[i]){ 435 oc=oc>>p->shiftoc; 436 437 if(oc>=P_BANDS)oc=P_BANDS-1; 438 if(oc<0)oc=0; 439 440 seed_curve(seed, 441 curves[oc], 442 max, 443 p->octave[i]-p->firstoc, 444 p->total_octave_lines, 445 p->eighth_octave_lines, 446 dBoffset); 447 } 448 } 449 } 450 seed_chase(float * seeds,int linesper,long n)451 static void seed_chase(float *seeds, int linesper, long n){ 452 long *posstack=alloca(n*sizeof(*posstack)); 453 float *ampstack=alloca(n*sizeof(*ampstack)); 454 long stack=0; 455 long pos=0; 456 long i; 457 458 for(i=0;i<n;i++){ 459 if(stack<2){ 460 posstack[stack]=i; 461 ampstack[stack++]=seeds[i]; 462 }else{ 463 while(1){ 464 if(seeds[i]<ampstack[stack-1]){ 465 posstack[stack]=i; 466 ampstack[stack++]=seeds[i]; 467 break; 468 }else{ 469 if(i<posstack[stack-1]+linesper){ 470 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] && 471 i<posstack[stack-2]+linesper){ 472 /* we completely overlap, making stack-1 irrelevant. pop it */ 473 stack--; 474 continue; 475 } 476 } 477 posstack[stack]=i; 478 ampstack[stack++]=seeds[i]; 479 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 for(i=0;i<stack;i++){ 490 long endpos; 491 if(i<stack-1 && ampstack[i+1]>ampstack[i]){ 492 endpos=posstack[i+1]; 493 }else{ 494 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is 495 discarded in short frames */ 496 } 497 if(endpos>n)endpos=n; 498 for(;pos<endpos;pos++) 499 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 } 506 507 /* bleaugh, this is more complicated than it needs to be */ 508 #include<stdio.h> max_seeds(vorbis_look_psy * p,float * seed,float * flr)509 static void max_seeds(vorbis_look_psy *p, 510 float *seed, 511 float *flr){ 512 long n=p->total_octave_lines; 513 int linesper=p->eighth_octave_lines; 514 long linpos=0; 515 long pos; 516 517 seed_chase(seed,linesper,n); /* for masking */ 518 519 pos=p->octave[0]-p->firstoc-(linesper>>1); 520 521 while(linpos+1<p->n){ 522 float minV=seed[pos]; 523 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc; 524 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit; 525 while(pos+1<=end){ 526 pos++; 527 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF) 528 minV=seed[pos]; 529 } 530 531 end=pos+p->firstoc; 532 for(;linpos<p->n && p->octave[linpos]<=end;linpos++) 533 if(flr[linpos]<minV)flr[linpos]=minV; 534 } 535 536 { 537 float minV=seed[p->total_octave_lines-1]; 538 for(;linpos<p->n;linpos++) 539 if(flr[linpos]<minV)flr[linpos]=minV; 540 } 541 542 } 543 bark_noise_hybridmp(int n,const long * b,const float * f,float * noise,const float offset,const int fixed)544 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 float *N=alloca(n*sizeof(*N)); 551 float *X=alloca(n*sizeof(*N)); 552 float *XX=alloca(n*sizeof(*N)); 553 float *Y=alloca(n*sizeof(*N)); 554 float *XY=alloca(n*sizeof(*N)); 555 556 float tN, tX, tXX, tY, tXY; 557 int i; 558 559 int lo, hi; 560 float R=0.f; 561 float A=0.f; 562 float B=0.f; 563 float D=1.f; 564 float w, x, y; 565 566 tN = tX = tXX = tY = tXY = 0.f; 567 568 y = f[0] + offset; 569 if (y < 1.f) y = 1.f; 570 571 w = y * y * .5; 572 573 tN += w; 574 tX += w; 575 tY += w * y; 576 577 N[0] = tN; 578 X[0] = tX; 579 XX[0] = tXX; 580 Y[0] = tY; 581 XY[0] = tXY; 582 583 for (i = 1, x = 1.f; i < n; i++, x += 1.f) { 584 585 y = f[i] + offset; 586 if (y < 1.f) y = 1.f; 587 588 w = y * y; 589 590 tN += w; 591 tX += w * x; 592 tXX += w * x * x; 593 tY += w * y; 594 tXY += w * x * y; 595 596 N[i] = tN; 597 X[i] = tX; 598 XX[i] = tXX; 599 Y[i] = tY; 600 XY[i] = tXY; 601 } 602 603 for (i = 0, x = 0.f;; i++, x += 1.f) { 604 605 lo = b[i] >> 16; 606 if( lo>=0 ) break; 607 hi = b[i] & 0xffff; 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) 620 R = 0.f; 621 622 noise[i] = R - offset; 623 } 624 625 for ( ;; i++, x += 1.f) { 626 627 lo = b[i] >> 16; 628 hi = b[i] & 0xffff; 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 for ( ; i < n; i++, x += 1.f) { 646 647 R = (A + x * B) / D; 648 if (R < 0.f) R = 0.f; 649 650 noise[i] = R - offset; 651 } 652 653 if (fixed <= 0) return; 654 655 for (i = 0, x = 0.f;; i++, x += 1.f) { 656 hi = i + fixed / 2; 657 lo = hi - fixed; 658 if(lo>=0)break; 659 660 tN = N[hi] + N[-lo]; 661 tX = X[hi] - X[-lo]; 662 tXX = XX[hi] + XX[-lo]; 663 tY = Y[hi] + Y[-lo]; 664 tXY = XY[hi] - XY[-lo]; 665 666 667 A = tY * tXX - tX * tXY; 668 B = tN * tXY - tX * tY; 669 D = tN * tXX - tX * tX; 670 R = (A + x * B) / D; 671 672 if (R - offset < noise[i]) noise[i] = R - offset; 673 } 674 for ( ;; i++, x += 1.f) { 675 676 hi = i + fixed / 2; 677 lo = hi - fixed; 678 if(hi>=n)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=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=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]=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=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 = 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 float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */ 1024 1025 /* mdct is our raw mdct output, floor not removed. */ 1026 /* inout passes in the ifloor, passes back quantized result */ 1027 1028 /* unquantized energy (negative indicates amplitude has negative sign) */ 1029 float **raw = alloca(ch*sizeof(*raw)); 1030 1031 /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */ 1032 float **quant = alloca(ch*sizeof(*quant)); 1033 1034 /* floor energy */ 1035 float **floor = alloca(ch*sizeof(*floor)); 1036 1037 /* flags indicating raw/quantized status of elements in raw vector */ 1038 int **flag = alloca(ch*sizeof(*flag)); 1039 1040 /* non-zero flag working vector */ 1041 int *nz = alloca(ch*sizeof(*nz)); 1042 1043 /* energy surplus/defecit tracking */ 1044 float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc)); 1045 1046 /* The threshold of a stereo is changed with the size of n */ 1047 if(n > 1000) 1048 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]]; 1049 1050 raw[0] = alloca(ch*partition*sizeof(**raw)); 1051 quant[0] = alloca(ch*partition*sizeof(**quant)); 1052 floor[0] = alloca(ch*partition*sizeof(**floor)); 1053 flag[0] = alloca(ch*partition*sizeof(**flag)); 1054 1055 for(i=1;i<ch;i++){ 1056 raw[i] = &raw[0][partition*i]; 1057 quant[i] = &quant[0][partition*i]; 1058 floor[i] = &floor[0][partition*i]; 1059 flag[i] = &flag[0][partition*i]; 1060 } 1061 for(i=0;i<ch+vi->coupling_steps;i++) 1062 acc[i]=0.f; 1063 1064 for(i=0;i<n;i+=partition){ 1065 int k,j,jn = partition > n-i ? n-i : partition; 1066 int step,track = 0; 1067 1068 memcpy(nz,nonzero,sizeof(*nz)*ch); 1069 1070 /* prefill */ 1071 memset(flag[0],0,ch*partition*sizeof(**flag)); 1072 for(k=0;k<ch;k++){ 1073 int *iout = &iwork[k][i]; 1074 if(nz[k]){ 1075 1076 for(j=0;j<jn;j++) 1077 floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]]; 1078 1079 flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn); 1080 1081 for(j=0;j<jn;j++){ 1082 quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j]; 1083 if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f; 1084 floor[k][j]*=floor[k][j]; 1085 } 1086 1087 acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout); 1088 1089 }else{ 1090 for(j=0;j<jn;j++){ 1091 floor[k][j] = 1e-10f; 1092 raw[k][j] = 0.f; 1093 quant[k][j] = 0.f; 1094 flag[k][j] = 0; 1095 iout[j]=0; 1096 } 1097 acc[track]=0.f; 1098 } 1099 track++; 1100 } 1101 1102 /* coupling */ 1103 for(step=0;step<vi->coupling_steps;step++){ 1104 int Mi = vi->coupling_mag[step]; 1105 int Ai = vi->coupling_ang[step]; 1106 int *iM = &iwork[Mi][i]; 1107 int *iA = &iwork[Ai][i]; 1108 float *reM = raw[Mi]; 1109 float *reA = raw[Ai]; 1110 float *qeM = quant[Mi]; 1111 float *qeA = quant[Ai]; 1112 float *floorM = floor[Mi]; 1113 float *floorA = floor[Ai]; 1114 int *fM = flag[Mi]; 1115 int *fA = flag[Ai]; 1116 1117 if(nz[Mi] || nz[Ai]){ 1118 nz[Mi] = nz[Ai] = 1; 1119 1120 for(j=0;j<jn;j++){ 1121 1122 if(j<sliding_lowpass-i){ 1123 if(fM[j] || fA[j]){ 1124 /* lossless coupling */ 1125 1126 reM[j] = fabs(reM[j])+fabs(reA[j]); 1127 qeM[j] = qeM[j]+qeA[j]; 1128 fM[j]=fA[j]=1; 1129 1130 /* couple iM/iA */ 1131 { 1132 int A = iM[j]; 1133 int B = iA[j]; 1134 1135 if(abs(A)>abs(B)){ 1136 iA[j]=(A>0?A-B:B-A); 1137 }else{ 1138 iA[j]=(B>0?A-B:B-A); 1139 iM[j]=B; 1140 } 1141 1142 /* collapse two equivalent tuples to one */ 1143 if(iA[j]>=abs(iM[j])*2){ 1144 iA[j]= -iA[j]; 1145 iM[j]= -iM[j]; 1146 } 1147 1148 } 1149 1150 }else{ 1151 /* lossy (point) coupling */ 1152 if(j<limit-i){ 1153 /* dipole */ 1154 reM[j] += reA[j]; 1155 qeM[j] = fabs(reM[j]); 1156 }else{ 1157 /* AoTuV */ 1158 /** @ M2 ** 1159 The boost problem by the combination of noise normalization and point stereo is eased. 1160 However, this is a temporary patch. 1161 by Aoyumi @ 2004/04/18 1162 */ 1163 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit))); 1164 1165 /* elliptical */ 1166 if(reM[j]+reA[j]<0){ 1167 reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); 1168 }else{ 1169 reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); 1170 } 1171 } 1172 reA[j]=qeA[j]=0.f; 1173 fA[j]=1; 1174 iA[j]=0; 1175 } 1176 } 1177 floorM[j]=floorA[j]=floorM[j]+floorA[j]; 1178 } 1179 /* normalize the resulting mag vector */ 1180 acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM); 1181 track++; 1182 } 1183 } 1184 } 1185 1186 for(i=0;i<vi->coupling_steps;i++){ 1187 /* make sure coupling a zero and a nonzero channel results in two 1188 nonzero channels. */ 1189 if(nonzero[vi->coupling_mag[i]] || 1190 nonzero[vi->coupling_ang[i]]){ 1191 nonzero[vi->coupling_mag[i]]=1; 1192 nonzero[vi->coupling_ang[i]]=1; 1193 } 1194 } 1195 } 1196