codecs/ilbc/enhancer.c

   1
   2 /******************************************************************
   3
   4     iLBC Speech Coder ANSI-C Source Code
   5
   6     enhancer.c
   7
   8     Copyright (C) The Internet Society (2004).
   9     All Rights Reserved.
  10
  11 ******************************************************************/
  12
  13 #include <math.h>
  14 #include <string.h>
  15 #include "iLBC_define.h"
  16 #include "enhancer.h"
  17 #include "constants.h"
  18 #include "filter.h"
  19
  20 /*----------------------------------------------------------------*
  21
  22
  23  * Find index in array such that the array element with said
  24  * index is the element of said array closest to "value"
  25  * according to the squared-error criterion
  26  *---------------------------------------------------------------*/
  27
  28 static void NearestNeighbor(
  29     int   *index,   /* (o) index of array element closest
  30                            to value */
  31     float *array,   /* (i) data array */
  32     float value,/* (i) value */
  33     int arlength/* (i) dimension of data array */
  34 ){
  35     int i;
  36     float bestcrit,crit;
  37
  38     crit=array[0]-value;
  39     bestcrit=crit*crit;
  40     *index=0;
  41     for (i=1; i<arlength; i++) {
  42         crit=array[i]-value;
  43         crit=crit*crit;
  44
  45         if (crit<bestcrit) {
  46             bestcrit=crit;
  47             *index=i;
  48         }
  49     }
  50 }
  51
  52 /*----------------------------------------------------------------*
  53  * compute cross correlation between sequences
  54  *---------------------------------------------------------------*/
  55
  56 static void mycorr1(
  57     float* corr,    /* (o) correlation of seq1 and seq2 */
  58     float* seq1,    /* (i) first sequence */
  59     int dim1,           /* (i) dimension first seq1 */
  60     const float *seq2,  /* (i) second sequence */
  61     int dim2        /* (i) dimension seq2 */
  62 ){
  63     int i,j;
  64
  65     for (i=0; i<=dim1-dim2; i++) {
  66         corr[i]=0.0;
  67         for (j=0; j<dim2; j++) {
  68             corr[i] += seq1[i+j] * seq2[j];
  69         }
  70     }
  71 }
  72
  73 /*----------------------------------------------------------------*
  74  * upsample finite array assuming zeros outside bounds
  75  *---------------------------------------------------------------*/
  76
  77
  78
  79 static void enh_upsample(
  80     float* useq1,   /* (o) upsampled output sequence */
  81     float* seq1,/* (i) unupsampled sequence */
  82     int dim1,       /* (i) dimension seq1 */
  83     int hfl         /* (i) polyphase filter length=2*hfl+1 */
  84 ){
  85     float *pu,*ps;
  86     int i,j,k,q,filterlength,hfl2;
  87     const float *polyp[ENH_UPS0]; /* pointers to
  88                                      polyphase columns */
  89     const float *pp;
  90
  91     /* define pointers for filter */
  92
  93     filterlength=2*hfl+1;
  94
  95     if ( filterlength > dim1 ) {
  96         hfl2=(int) (dim1/2);
  97         for (j=0; j<ENH_UPS0; j++) {
  98             polyp[j]=polyphaserTbl+j*filterlength+hfl-hfl2;
  99         }
 100         hfl=hfl2;
 101         filterlength=2*hfl+1;
 102     }
 103     else {
 104         for (j=0; j<ENH_UPS0; j++) {
 105             polyp[j]=polyphaserTbl+j*filterlength;
 106         }
 107     }
 108
 109     /* filtering: filter overhangs left side of sequence */
 110
 111     pu=useq1;
 112     for (i=hfl; i<filterlength; i++) {
 113         for (j=0; j<ENH_UPS0; j++) {
 114             *pu=0.0;
 115             pp = polyp[j];
 116             ps = seq1+i;
 117             for (k=0; k<=i; k++) {
 118                 *pu += *ps-- * *pp++;
 119             }
 120             pu++;
 121         }
 122     }
 123
 124     /* filtering: simple convolution=inner products */
 125
 126     for (i=filterlength; i<dim1; i++) {
 127         for (j=0;j<ENH_UPS0; j++){
 128             *pu=0.0;
 129             pp = polyp[j];
 130             ps = seq1+i;
 131             for (k=0; k<filterlength; k++) {
 132                 *pu += *ps-- * *pp++;
 133
 134
 135             }
 136             pu++;
 137         }
 138     }
 139
 140     /* filtering: filter overhangs right side of sequence */
 141
 142     for (q=1; q<=hfl; q++) {
 143         for (j=0; j<ENH_UPS0; j++) {
 144             *pu=0.0;
 145             pp = polyp[j]+q;
 146             ps = seq1+dim1-1;
 147             for (k=0; k<filterlength-q; k++) {
 148                 *pu += *ps-- * *pp++;
 149             }
 150             pu++;
 151         }
 152     }
 153 }
 154
 155
 156 /*----------------------------------------------------------------*
 157  * find segment starting near idata+estSegPos that has highest
 158  * correlation with idata+centerStartPos through
 159  * idata+centerStartPos+ENH_BLOCKL-1 segment is found at a
 160  * resolution of ENH_UPSO times the original of the original
 161  * sampling rate
 162  *---------------------------------------------------------------*/
 163
 164 static void refiner(
 165     float *seg,         /* (o) segment array */
 166     float *updStartPos, /* (o) updated start point */
 167     float* idata,       /* (i) original data buffer */
 168     int idatal,         /* (i) dimension of idata */
 169     int centerStartPos, /* (i) beginning center segment */
 170     float estSegPos,/* (i) estimated beginning other segment */
 171     float period    /* (i) estimated pitch period */
 172 ){
 173     int estSegPosRounded,searchSegStartPos,searchSegEndPos,corrdim;
 174     int tloc,tloc2,i,st,en,fraction;
 175     float vect[ENH_VECTL],corrVec[ENH_CORRDIM],maxv;
 176     float corrVecUps[ENH_CORRDIM*ENH_UPS0];
 177
 178     /* defining array bounds */
 179
 180     estSegPosRounded=(int)(estSegPos - 0.5);
 181
 182     searchSegStartPos=estSegPosRounded-ENH_SLOP;
 183
 184     if (searchSegStartPos<0) {
 185         searchSegStartPos=0;
 186     }
 187     searchSegEndPos=estSegPosRounded+ENH_SLOP;
 188
 189
 190
 191     if (searchSegEndPos+ENH_BLOCKL >= idatal) {
 192         searchSegEndPos=idatal-ENH_BLOCKL-1;
 193     }
 194     corrdim=searchSegEndPos-searchSegStartPos+1;
 195
 196     /* compute upsampled correlation (corr33) and find
 197        location of max */
 198
 199     mycorr1(corrVec,idata+searchSegStartPos,
 200         corrdim+ENH_BLOCKL-1,idata+centerStartPos,ENH_BLOCKL);
 201     enh_upsample(corrVecUps,corrVec,corrdim,ENH_FL0);
 202     tloc=0; maxv=corrVecUps[0];
 203     for (i=1; i<ENH_UPS0*corrdim; i++) {
 204
 205         if (corrVecUps[i]>maxv) {
 206             tloc=i;
 207             maxv=corrVecUps[i];
 208         }
 209     }
 210
 211     /* make vector can be upsampled without ever running outside
 212        bounds */
 213
 214     *updStartPos= (float)searchSegStartPos +
 215         (float)tloc/(float)ENH_UPS0+(float)1.0;
 216     tloc2=(int)(tloc/ENH_UPS0);
 217
 218     if (tloc>tloc2*ENH_UPS0) {
 219         tloc2++;
 220     }
 221     st=searchSegStartPos+tloc2-ENH_FL0;
 222
 223     if (st<0) {
 224         memset(vect,0,-st*sizeof(float));
 225         memcpy(&vect[-st],idata, (ENH_VECTL+st)*sizeof(float));
 226     }
 227     else {
 228         en=st+ENH_VECTL;
 229
 230         if (en>idatal) {
 231             memcpy(vect, &idata[st],
 232                 (ENH_VECTL-(en-idatal))*sizeof(float));
 233             memset(&vect[ENH_VECTL-(en-idatal)], 0,
 234                 (en-idatal)*sizeof(float));
 235         }
 236         else {
 237             memcpy(vect, &idata[st], ENH_VECTL*sizeof(float));
 238         }
 239     }
 240     fraction=tloc2*ENH_UPS0-tloc;
 241
 242     /* compute the segment (this is actually a convolution) */
 243
 244     mycorr1(seg,vect,ENH_VECTL,polyphaserTbl+(2*ENH_FL0+1)*fraction,
 245
 246
 247         2*ENH_FL0+1);
 248 }
 249
 250 /*----------------------------------------------------------------*
 251  * find the smoothed output data
 252  *---------------------------------------------------------------*/
 253
 254 static void smath(
 255     float *odata,   /* (o) smoothed output */
 256     float *sseq,/* (i) said second sequence of waveforms */
 257     int hl,         /* (i) 2*hl+1 is sseq dimension */
 258     float alpha0/* (i) max smoothing energy fraction */
 259 ){
 260     int i,k;
 261     float w00,w10,w11,A,B,C,*psseq,err,errs;
 262     float surround[BLOCKL_MAX]; /* shape contributed by other than
 263                                    current */
 264     float wt[2*ENH_HL+1];       /* waveform weighting to get
 265                                    surround shape */
 266     float denom;
 267
 268     /* create shape of contribution from all waveforms except the
 269        current one */
 270
 271     for (i=1; i<=2*hl+1; i++) {
 272         wt[i-1] = (float)0.5*(1 - (float)cos(2*PI*i/(2*hl+2)));
 273     }
 274     wt[hl]=0.0; /* for clarity, not used */
 275     for (i=0; i<ENH_BLOCKL; i++) {
 276         surround[i]=sseq[i]*wt[0];
 277     }
 278     for (k=1; k<hl; k++) {
 279         psseq=sseq+k*ENH_BLOCKL;
 280         for(i=0;i<ENH_BLOCKL; i++) {
 281             surround[i]+=psseq[i]*wt[k];
 282         }
 283     }
 284     for (k=hl+1; k<=2*hl; k++) {
 285         psseq=sseq+k*ENH_BLOCKL;
 286         for(i=0;i<ENH_BLOCKL; i++) {
 287             surround[i]+=psseq[i]*wt[k];
 288         }
 289     }
 290
 291     /* compute some inner products */
 292
 293     w00 = w10 = w11 = 0.0;
 294     psseq=sseq+hl*ENH_BLOCKL; /* current block  */
 295     for (i=0; i<ENH_BLOCKL;i++) {
 296         w00+=psseq[i]*psseq[i];
 297         w11+=surround[i]*surround[i];
 298         w10+=surround[i]*psseq[i];
 299     }
 300
 301
 302
 303     if (fabs(w11) < 1.0) {
 304         w11=1.0;
 305     }
 306     C = (float)sqrt( w00/w11);
 307
 308     /* first try enhancement without power-constraint */
 309
 310     errs=0.0;
 311     psseq=sseq+hl*ENH_BLOCKL;
 312     for (i=0; i<ENH_BLOCKL; i++) {
 313         odata[i]=C*surround[i];
 314         err=psseq[i]-odata[i];
 315         errs+=err*err;
 316     }
 317
 318     /* if constraint violated by first try, add constraint */
 319
 320     if (errs > alpha0 * w00) {
 321         if ( w00 < 1) {
 322             w00=1;
 323         }
 324         denom = (w11*w00-w10*w10)/(w00*w00);
 325
 326         if (denom > 0.0001) { /* eliminates numerical problems
 327                                  for if smooth */
 328             A = (float)sqrt( (alpha0- alpha0*alpha0/4)/denom);
 329             B = -alpha0/2 - A * w10/w00;
 330             B = B+1;
 331         }
 332         else { /* essentially no difference between cycles;
 333                   smoothing not needed */
 334             A= 0.0;
 335             B= 1.0;
 336         }
 337
 338         /* create smoothed sequence */
 339
 340         psseq=sseq+hl*ENH_BLOCKL;
 341         for (i=0; i<ENH_BLOCKL; i++) {
 342             odata[i]=A*surround[i]+B*psseq[i];
 343         }
 344     }
 345 }
 346
 347 /*----------------------------------------------------------------*
 348  * get the pitch-synchronous sample sequence
 349  *---------------------------------------------------------------*/
 350
 351 static void getsseq(
 352     float *sseq,    /* (o) the pitch-synchronous sequence */
 353     float *idata,       /* (i) original data */
 354     int idatal,         /* (i) dimension of data */
 355     int centerStartPos, /* (i) where current block starts */
 356     float *period,      /* (i) rough-pitch-period array */
 357
 358
 359     float *plocs,       /* (i) where periods of period array
 360                                are taken */
 361     int periodl,    /* (i) dimension period array */
 362     int hl              /* (i) 2*hl+1 is the number of sequences */
 363 ){
 364     int i,centerEndPos,q;
 365     float blockStartPos[2*ENH_HL+1];
 366     int lagBlock[2*ENH_HL+1];
 367     float plocs2[ENH_PLOCSL];
 368     float *psseq;
 369
 370     centerEndPos=centerStartPos+ENH_BLOCKL-1;
 371
 372     /* present */
 373
 374     NearestNeighbor(lagBlock+hl,plocs,
 375         (float)0.5*(centerStartPos+centerEndPos),periodl);
 376
 377     blockStartPos[hl]=(float)centerStartPos;
 378     psseq=sseq+ENH_BLOCKL*hl;
 379     memcpy(psseq, idata+centerStartPos, ENH_BLOCKL*sizeof(float));
 380
 381     /* past */
 382
 383     for (q=hl-1; q>=0; q--) {
 384         blockStartPos[q]=blockStartPos[q+1]-period[lagBlock[q+1]];
 385         NearestNeighbor(lagBlock+q,plocs,
 386             blockStartPos[q]+
 387             ENH_BLOCKL_HALF-period[lagBlock[q+1]], periodl);
 388
 389
 390         if (blockStartPos[q]-ENH_OVERHANG>=0) {
 391             refiner(sseq+q*ENH_BLOCKL, blockStartPos+q, idata,
 392                 idatal, centerStartPos, blockStartPos[q],
 393                 period[lagBlock[q+1]]);
 394         } else {
 395             psseq=sseq+q*ENH_BLOCKL;
 396             memset(psseq, 0, ENH_BLOCKL*sizeof(float));
 397         }
 398     }
 399
 400     /* future */
 401
 402     for (i=0; i<periodl; i++) {
 403         plocs2[i]=plocs[i]-period[i];
 404     }
 405     for (q=hl+1; q<=2*hl; q++) {
 406         NearestNeighbor(lagBlock+q,plocs2,
 407             blockStartPos[q-1]+ENH_BLOCKL_HALF,periodl);
 408
 409         blockStartPos[q]=blockStartPos[q-1]+period[lagBlock[q]];
 410         if (blockStartPos[q]+ENH_BLOCKL+ENH_OVERHANG<idatal) {
 411             refiner(sseq+ENH_BLOCKL*q, blockStartPos+q, idata,
 412                 idatal, centerStartPos, blockStartPos[q],
 413
 414
 415                 period[lagBlock[q]]);
 416         }
 417         else {
 418             psseq=sseq+q*ENH_BLOCKL;
 419             memset(psseq, 0, ENH_BLOCKL*sizeof(float));
 420         }
 421     }
 422 }
 423
 424 /*----------------------------------------------------------------*
 425  * perform enhancement on idata+centerStartPos through
 426  * idata+centerStartPos+ENH_BLOCKL-1
 427  *---------------------------------------------------------------*/
 428
 429 static void enhancer(
 430     float *odata,       /* (o) smoothed block, dimension blockl */
 431     float *idata,       /* (i) data buffer used for enhancing */
 432     int idatal,         /* (i) dimension idata */
 433     int centerStartPos, /* (i) first sample current block
 434                                within idata */
 435     float alpha0,       /* (i) max correction-energy-fraction
 436                               (in [0,1]) */
 437     float *period,      /* (i) pitch period array */
 438     float *plocs,       /* (i) locations where period array
 439                                values valid */
 440     int periodl         /* (i) dimension of period and plocs */
 441 ){
 442     float sseq[(2*ENH_HL+1)*ENH_BLOCKL];
 443
 444     /* get said second sequence of segments */
 445
 446     getsseq(sseq,idata,idatal,centerStartPos,period,
 447         plocs,periodl,ENH_HL);
 448
 449     /* compute the smoothed output from said second sequence */
 450
 451     smath(odata,sseq,ENH_HL,alpha0);
 452
 453 }
 454
 455 /*----------------------------------------------------------------*
 456  * cross correlation
 457  *---------------------------------------------------------------*/
 458
 459 float xCorrCoef(
 460     float *target,      /* (i) first array */
 461     float *regressor,   /* (i) second array */
 462     int subl        /* (i) dimension arrays */
 463 ){
 464     int i;
 465     float ftmp1, ftmp2;
 466
 467     ftmp1 = 0.0;
 468     ftmp2 = 0.0;
 469
 470
 471     for (i=0; i<subl; i++) {
 472         ftmp1 += target[i]*regressor[i];
 473         ftmp2 += regressor[i]*regressor[i];
 474     }
 475
 476     if (ftmp1 > 0.0) {
 477         return (float)(ftmp1*ftmp1/ftmp2);
 478     }
 479     else {
 480         return (float)0.0;
 481     }
 482 }
 483
 484 /*----------------------------------------------------------------*
 485  * interface for enhancer
 486  *---------------------------------------------------------------*/
 487
 488 int enhancerInterface(
 489     float *out,                     /* (o) enhanced signal */
 490     float *in,                      /* (i) unenhanced signal */
 491     iLBC_Dec_Inst_t *iLBCdec_inst   /* (i) buffers etc */
 492 ){
 493     float *enh_buf, *enh_period;
 494     int iblock, isample;
 495     int lag=0, ilag, i, ioffset;
 496     float cc, maxcc;
 497     float ftmp1, ftmp2;
 498     float *inPtr, *enh_bufPtr1, *enh_bufPtr2;
 499     float plc_pred[ENH_BLOCKL];
 500
 501     float lpState[6], downsampled[(ENH_NBLOCKS*ENH_BLOCKL+120)/2];
 502     int inLen=ENH_NBLOCKS*ENH_BLOCKL+120;
 503     int start, plc_blockl, inlag;
 504
 505     enh_buf=iLBCdec_inst->enh_buf;
 506     enh_period=iLBCdec_inst->enh_period;
 507
 508     memmove(enh_buf, &enh_buf[iLBCdec_inst->blockl],
 509         (ENH_BUFL-iLBCdec_inst->blockl)*sizeof(float));
 510
 511     memcpy(&enh_buf[ENH_BUFL-iLBCdec_inst->blockl], in,
 512         iLBCdec_inst->blockl*sizeof(float));
 513
 514     if (iLBCdec_inst->mode==30)
 515         plc_blockl=ENH_BLOCKL;
 516     else
 517         plc_blockl=40;
 518
 519     /* when 20 ms frame, move processing one block */
 520     ioffset=0;
 521     if (iLBCdec_inst->mode==20) ioffset=1;
 522
 523     i=3-ioffset;
 524     memmove(enh_period, &enh_period[i],
 525
 526
 527         (ENH_NBLOCKS_TOT-i)*sizeof(float));
 528
 529     /* Set state information to the 6 samples right before
 530        the samples to be downsampled. */
 531
 532     memcpy(lpState,
 533         enh_buf+(ENH_NBLOCKS_EXTRA+ioffset)*ENH_BLOCKL-126,
 534         6*sizeof(float));
 535
 536     /* Down sample a factor 2 to save computations */
 537
 538     DownSample(enh_buf+(ENH_NBLOCKS_EXTRA+ioffset)*ENH_BLOCKL-120,
 539                 lpFilt_coefsTbl, inLen-ioffset*ENH_BLOCKL,
 540                 lpState, downsampled);
 541
 542     /* Estimate the pitch in the down sampled domain. */
 543     for (iblock = 0; iblock<ENH_NBLOCKS-ioffset; iblock++) {
 544
 545         lag = 10;
 546         maxcc = xCorrCoef(downsampled+60+iblock*
 547             ENH_BLOCKL_HALF, downsampled+60+iblock*
 548             ENH_BLOCKL_HALF-lag, ENH_BLOCKL_HALF);
 549         for (ilag=11; ilag<60; ilag++) {
 550             cc = xCorrCoef(downsampled+60+iblock*
 551                 ENH_BLOCKL_HALF, downsampled+60+iblock*
 552                 ENH_BLOCKL_HALF-ilag, ENH_BLOCKL_HALF);
 553
 554             if (cc > maxcc) {
 555                 maxcc = cc;
 556                 lag = ilag;
 557             }
 558         }
 559
 560         /* Store the estimated lag in the non-downsampled domain */
 561         enh_period[iblock+ENH_NBLOCKS_EXTRA+ioffset] = (float)lag*2;
 562
 563
 564     }
 565
 566
 567     /* PLC was performed on the previous packet */
 568     if (iLBCdec_inst->prev_enh_pl==1) {
 569
 570         inlag=(int)enh_period[ENH_NBLOCKS_EXTRA+ioffset];
 571
 572         lag = inlag-1;
 573         maxcc = xCorrCoef(in, in+lag, plc_blockl);
 574         for (ilag=inlag; ilag<=inlag+1; ilag++) {
 575             cc = xCorrCoef(in, in+ilag, plc_blockl);
 576
 577             if (cc > maxcc) {
 578                 maxcc = cc;
 579                 lag = ilag;
 580             }
 581         }
 582
 583
 584
 585         enh_period[ENH_NBLOCKS_EXTRA+ioffset-1]=(float)lag;
 586
 587         /* compute new concealed residual for the old lookahead,
 588            mix the forward PLC with a backward PLC from
 589            the new frame */
 590
 591         inPtr=&in[lag-1];
 592
 593         enh_bufPtr1=&plc_pred[plc_blockl-1];
 594
 595         if (lag>plc_blockl) {
 596             start=plc_blockl;
 597         } else {
 598             start=lag;
 599         }
 600
 601         for (isample = start; isample>0; isample--) {
 602             *enh_bufPtr1-- = *inPtr--;
 603         }
 604
 605         enh_bufPtr2=&enh_buf[ENH_BUFL-1-iLBCdec_inst->blockl];
 606         for (isample = (plc_blockl-1-lag); isample>=0; isample--)
 607 {
 608             *enh_bufPtr1-- = *enh_bufPtr2--;
 609         }
 610
 611         /* limit energy change */
 612         ftmp2=0.0;
 613         ftmp1=0.0;
 614         for (i=0;i<plc_blockl;i++) {
 615             ftmp2+=enh_buf[ENH_BUFL-1-iLBCdec_inst->blockl-i]*
 616                 enh_buf[ENH_BUFL-1-iLBCdec_inst->blockl-i];
 617             ftmp1+=plc_pred[i]*plc_pred[i];
 618         }
 619         ftmp1=(float)sqrt(ftmp1/(float)plc_blockl);
 620         ftmp2=(float)sqrt(ftmp2/(float)plc_blockl);
 621         if (ftmp1>(float)2.0*ftmp2 && ftmp1>0.0) {
 622             for (i=0;i<plc_blockl-10;i++) {
 623                 plc_pred[i]*=(float)2.0*ftmp2/ftmp1;
 624             }
 625             for (i=plc_blockl-10;i<plc_blockl;i++) {
 626                 plc_pred[i]*=(float)(i-plc_blockl+10)*
 627                     ((float)1.0-(float)2.0*ftmp2/ftmp1)/(float)(10)+
 628                     (float)2.0*ftmp2/ftmp1;
 629             }
 630         }
 631
 632         enh_bufPtr1=&enh_buf[ENH_BUFL-1-iLBCdec_inst->blockl];
 633         for (i=0; i<plc_blockl; i++) {
 634             ftmp1 = (float) (i+1) / (float) (plc_blockl+1);
 635             *enh_bufPtr1 *= ftmp1;
 636             *enh_bufPtr1 += ((float)1.0-ftmp1)*
 637                                 plc_pred[plc_blockl-1-i];
 638             enh_bufPtr1--;
 639         }
 640
 641
 642     }
 643
 644     if (iLBCdec_inst->mode==20) {
 645         /* Enhancer with 40 samples delay */
 646         for (iblock = 0; iblock<2; iblock++) {
 647             enhancer(out+iblock*ENH_BLOCKL, enh_buf,
 648                 ENH_BUFL, (5+iblock)*ENH_BLOCKL+40,
 649                 ENH_ALPHA0, enh_period, enh_plocsTbl,
 650                     ENH_NBLOCKS_TOT);
 651         }
 652     } else if (iLBCdec_inst->mode==30) {
 653         /* Enhancer with 80 samples delay */
 654         for (iblock = 0; iblock<3; iblock++) {
 655             enhancer(out+iblock*ENH_BLOCKL, enh_buf,
 656                 ENH_BUFL, (4+iblock)*ENH_BLOCKL,
 657                 ENH_ALPHA0, enh_period, enh_plocsTbl,
 658                     ENH_NBLOCKS_TOT);
 659         }
 660     }
 661
 662     return (lag*2);
 663 }
 664
 665