libghost/ceft.c

   1 /* Copyright (C) 2007
   2
   3    Code-Excited Fourier Transform -- This is highly experimental and
   4    it's not clear at all it even has a remote chance of working
   5
   6
   7    This library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Lesser General Public
   9    License as published by the Free Software Foundation; either
  10    version 2.1 of the License, or (at your option) any later version.
  11
  12    This library is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    Lesser General Public License for more details.
  16
  17    You should have received a copy of the GNU Lesser General Public
  18    License along with this library; if not, write to the Free Software
  19    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 */
  21
  22 #include "ceft.h"
  23 #include <stdlib.h>
  24 #include <stdio.h>
  25 #include <math.h>
  26
  27 #include "fftwrap.h"
  28
  29 #define NBANDS 15
  30 int qbank[] =   {1, 2, 4, 6, 8, 12, 16, 20, 24, 28, 36, 44, 52, 68, 84, 116, 128};
  31 //int qpulses[] = {2, 2, 2, 2, 2,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1};
  32 //int qpulses[] = {2, 2, 2, 2, 2,  2,  1,  1,  1,  1,  0,  0,  0,  0,  0};
  33 //int qpulses[] = {2, 3, 2, 2, 3,  2,  2,  2,  1,  2,  1,  0,  0,  0,  0};
  34 //int qpulses[] = {3, 4, 3, 2, 3,  2,  2,  2,  1,  2,  1,  0,  0,  0,  0};
  35
  36 //32 kbps
  37 int qpulses[] = {3, 4, 4, 3, 3,  2,  2,  2,  2,  2,  2,  0,  0,  0,  0};
  38
  39 //44 kbps
  40 //int qpulses[] = {4, 7, 6, 4, 4,  3,  3,  3,  3,  3,  3,  3,  0,  0,  0};
  41
  42 //int qpulses[] = {5, 5, 5, 5, 5,  5,  2,  2,  1,  2,  1,  0,  0,  0,  0};
  43 //int qpulses[] = {5, 5, 2, 2, 3,  2,  5,  5,  5,  5,  5,  5,  0,  0,  0};
  44
  45
  46 //int qpulses[] = {4, 4, 3, 3, 3,  3,  3,  2,  2,  3,  2,  2,  0,  0,  0};
  47 //int qpulses[] = {5, 5, 5, 5, 5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5};
  48
  49
  50 #define PBANDS 5
  51 int pbank[] = {1, 4, 8, 12, 20, 44};
  52
  53 void alg_quant(float *x, int N, int K, float *p)
  54 {
  55    float y[N];
  56    int i,j;
  57    float xy = 0;
  58    float yy = 0;
  59    float yp = 0;
  60    float Rpp=0;
  61    float gain=0;
  62    for (j=0;j<N;j++)
  63       Rpp += p[j]*p[j];
  64    for (i=0;i<N;i++)
  65       y[i] = 0;
  66
  67    for (i=0;i<K;i++)
  68    {
  69       int best_id=0;
  70       float max_val=-1e10;
  71       float best_xy=0, best_yy=0, best_yp = 0;
  72       for (j=0;j<N;j++)
  73       {
  74          float tmp_xy, tmp_yy, tmp_yp;
  75          float score;
  76          float g;
  77          tmp_xy = xy + fabs(x[j]);
  78          tmp_yy = yy + 2*fabs(y[j]) + 1;
  79          if (x[j]>0)
  80             tmp_yp = yp + p[j];
  81          else
  82             tmp_yp = yp - p[j];
  83          g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
  84          score = 2*g*tmp_xy - g*g*tmp_yy;
  85          if (score>max_val)
  86          {
  87             max_val = score;
  88             best_id = j;
  89             best_xy = tmp_xy;
  90             best_yy = tmp_yy;
  91             best_yp = tmp_yp;
  92             gain = g;
  93          }
  94       }
  95
  96       xy = best_xy;
  97       yy = best_yy;
  98       yp = best_yp;
  99       if (x[best_id]>0)
 100          y[best_id] += 1;
 101       else
 102          y[best_id] -= 1;
 103    }
 104
 105    for (i=0;i<N;i++)
 106       x[i] = p[i]+gain*y[i];
 107
 108 }
 109
 110 void alg_quant2(float *x, int N, int K, float *p)
 111 {
 112    int L = 5;
 113    float tata[200];
 114    float y[L][N];
 115    float tata2[200];
 116    float ny[L][N];
 117    int i, j, m;
 118    float xy[L], nxy[L];
 119    float yy[L], nyy[L];
 120    float yp[L], nyp[L];
 121    float best_scores[L];
 122    float Rpp=0;
 123    float gain[L];
 124    int maxL = 1;
 125    for (j=0;j<N;j++)
 126       Rpp += p[j]*p[j];
 127    for (m=0;m<L;m++)
 128       for (i=0;i<N;i++)
 129          y[m][i] = 0;
 130
 131    for (m=0;m<L;m++)
 132       for (i=0;i<N;i++)
 133          ny[m][i] = 0;
 134
 135    for (m=0;m<L;m++)
 136       xy[m] = yy[m] = yp[m] = gain[m] = 0;
 137
 138    for (i=0;i<K;i++)
 139    {
 140       int L2 = L;
 141       if (L>maxL)
 142       {
 143          L2 = maxL;
 144          maxL *= N;
 145       }
 146       for (m=0;m<L;m++)
 147          best_scores[m] = -1e10;
 148
 149       for (m=0;m<L2;m++)
 150       {
 151          for (j=0;j<N;j++)
 152          {
 153             //fprintf (stderr, "%d/%d %d/%d %d/%d\n", i, K, m, L2, j, N);
 154             float tmp_xy, tmp_yy, tmp_yp;
 155             float score;
 156             float g;
 157             tmp_xy = xy[m] + fabs(x[j]);
 158             tmp_yy = yy[m] + 2*fabs(y[m][j]) + 1;
 159             if (x[j]>0)
 160                tmp_yp = yp[m] + p[j];
 161             else
 162                tmp_yp = yp[m] - p[j];
 163             g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
 164             score = 2*g*tmp_xy - g*g*tmp_yy;
 165
 166             if (score>best_scores[L-1])
 167             {
 168                int k, n;
 169                int id = L-1;
 170                while (id > 0 && score > best_scores[id-1])
 171                   id--;
 172
 173                for (k=L-1;k>id;k--)
 174                {
 175                   nxy[k] = nxy[k-1];
 176                   nyy[k] = nyy[k-1];
 177                   nyp[k] = nyp[k-1];
 178                   //fprintf(stderr, "%d %d \n", N, k);
 179                   for (n=0;n<N;n++)
 180                      ny[k][n] = ny[k-1][n];
 181                   gain[k] = gain[k-1];
 182                   best_scores[k] = best_scores[k-1];
 183                }
 184
 185                nxy[id] = tmp_xy;
 186                nyy[id] = tmp_yy;
 187                nyp[id] = tmp_yp;
 188                gain[id] = g;
 189                for (n=0;n<N;n++)
 190                   ny[id][n] = y[m][n];
 191                if (x[j]>0)
 192                   ny[id][j] += 1;
 193                else
 194                   ny[id][j] -= 1;
 195                best_scores[id] = score;
 196             }
 197
 198          }
 199
 200       }
 201       int k,n;
 202       for (k=0;k<L;k++)
 203       {
 204          xy[k] = nxy[k];
 205          yy[k] = nyy[k];
 206          yp[k] = nyp[k];
 207          for (n=0;n<N;n++)
 208             y[k][n] = ny[k][n];
 209       }
 210
 211    }
 212
 213    for (i=0;i<N;i++)
 214       x[i] = p[i]+gain[0]*y[0][i];
 215
 216 }
 217
 218 void noise_quant(float *x, int N, int K, float *p)
 219 {
 220    int i;
 221    float E = 1e-10;
 222    for (i=0;i<N;i++)
 223    {
 224       x[i] = (rand()%1000)/500.-1;
 225       E += x[i]*x[i];
 226    }
 227    E = 1./sqrt(E);
 228    for (i=0;i<N;i++)
 229    {
 230       x[i] *= E;
 231    }
 232 }
 233
 234 void compute_bank(float *X, float *bank)
 235 {
 236    int i;
 237    for (i=0;i<NBANDS;i++)
 238    {
 239       int j;
 240       bank[i] = 1e-10;
 241       for (j=qbank[i];j<qbank[i+1];j++)
 242       {
 243          bank[i] += X[j*2-1]*X[j*2-1];
 244          bank[i] += X[j*2]*X[j*2];
 245       }
 246       //bank[i] = sqrt(.5*bank[i]/(qbank[i+1]-qbank[i]));
 247       bank[i] = sqrt(bank[i]);
 248    }
 249 }
 250
 251 void normalise_bank(float *X, float *bank)
 252 {
 253    int i;
 254    for (i=0;i<NBANDS;i++)
 255    {
 256       int j;
 257       float x = 1.f/bank[i];
 258       for (j=qbank[i];j<qbank[i+1];j++)
 259       {
 260          X[j*2-1] *= x;
 261          X[j*2]   *= x;
 262       }
 263    }
 264    for (i=2*qbank[NBANDS]-1;i<256;i++)
 265       X[i] = 0;
 266 }
 267
 268 void denormalise_bank(float *X, float *bank)
 269 {
 270    int i;
 271    for (i=0;i<NBANDS;i++)
 272    {
 273       int j;
 274       float x = bank[i];
 275       for (j=qbank[i];j<qbank[i+1];j++)
 276       {
 277          X[j*2-1] *= x;
 278          X[j*2]   *= x;
 279       }
 280    }
 281    //FIXME: Kludge
 282    X[255] = 0;
 283 }
 284
 285 float norm2(float *x, int len)
 286 {
 287    float E=0;
 288    int i;
 289    for (i=0;i<len;i++)
 290       E += x[i]*x[i];
 291    return E;
 292 }
 293
 294 void crappy_fft(float *X, int len, int R, int dir)
 295 {
 296    int i, j;
 297    float fact = 2*M_PI/(len>>1);
 298    float out[len];
 299    for (i=0;i<len>>1;i++)
 300    {
 301       out[2*i] = 0;
 302       out[2*i+1] = 0;
 303       for (j=0;j<len>>1;j++)
 304       {
 305          float c, d;
 306          c = cos(fact*j*i);
 307          if (dir == 1)
 308             d = -sin(fact*j*i);
 309          else
 310             d = -sin(fact*j*i);
 311
 312          out[2*i  ] += c*X[2*j] + d*X[2*j+1];
 313          out[2*i+1] += d*X[2*j] + c*X[2*j+1];
 314       }
 315    }
 316    for (i=0;i<len;i++)
 317       X[i] = out[i]/sqrt(len/2);
 318 }
 319
 320 void crappy_rfft(float *X, int len, int R, int dir)
 321 {
 322    int N=len*2;
 323    float x[N];
 324    int i;
 325    /*for (i=0;i<len;i++)
 326       printf ("%f ", X[i]);*/
 327    if (dir>0)
 328    {
 329       for (i=0;i<len;i++)
 330       {
 331          x[2*i] = X[i];
 332          x[2*i+1] = 0;
 333       }
 334       crappy_fft(x, N, R, 1);
 335       for (i=2;i<len;i++)
 336          X[i] = sqrt(2)*x[i];
 337       X[1] = x[len];
 338       X[0] = x[0];
 339    } else {
 340       for (i=1;i<len>>1;i++)
 341       {
 342          x[2*i] = sqrt(.5)*X[2*i];
 343          x[2*i+1] = sqrt(.5)*X[2*i+1];
 344          x[2*(len-i)] = sqrt(.5)*X[2*i];
 345          x[2*(len-i)+1] = -sqrt(.5)*X[2*i+1];
 346       }
 347       x[len] = X[1];
 348       x[len+1] = 0;
 349       x[1] = 0;
 350       x[0] = X[0];
 351       crappy_fft(x, N, R, -1);
 352       for (i=0;i<len;i++)
 353          X[i] = x[2*i];
 354    }
 355    /*printf ("  ");
 356    for (i=0;i<len;i++)
 357       printf ("%f ", X[i]);
 358    printf ("\n");*/
 359 }
 360
 361 void exp_rotation(float *X, int len, float theta, int dir)
 362 {
 363    int i;
 364    float c, s;
 365    c = cos(theta);
 366    s = sin(theta);
 367    if (dir > 0)
 368    {
 369       for (i=0;i<(len/2)-1;i++)
 370       {
 371          float x1, x2;
 372          x1 = X[2*i];
 373          x2 = X[2*i+2];
 374          X[2*i] = c*x1 - s*x2;
 375          X[2*i+2] = c*x2 + s*x1;
 376
 377          x1 = X[2*i+1];
 378          x2 = X[2*i+3];
 379          X[2*i+1] = c*x1 - s*x2;
 380          X[2*i+3] = c*x2 + s*x1;
 381       }
 382       for (i=(len/2)-3;i>=0;i--)
 383       {
 384          float x1, x2;
 385          x1 = X[2*i];
 386          x2 = X[2*i+2];
 387          X[2*i] = c*x1 - s*x2;
 388          X[2*i+2] = c*x2 + s*x1;
 389
 390          x1 = X[2*i+1];
 391          x2 = X[2*i+3];
 392          X[2*i+1] = c*x1 - s*x2;
 393          X[2*i+3] = c*x2 + s*x1;
 394       }
 395
 396    } else {
 397       for (i=0;i<(len/2)-2;i++)
 398       {
 399          float x1, x2;
 400          x1 = X[2*i];
 401          x2 = X[2*i+2];
 402          X[2*i] = c*x1 + s*x2;
 403          X[2*i+2] = c*x2 - s*x1;
 404
 405          x1 = X[2*i+1];
 406          x2 = X[2*i+3];
 407          X[2*i+1] = c*x1 + s*x2;
 408          X[2*i+3] = c*x2 - s*x1;
 409       }
 410
 411       for (i=(len/2)-2;i>=0;i--)
 412       {
 413          float x1, x2;
 414          x1 = X[2*i];
 415          x2 = X[2*i+2];
 416          X[2*i] = c*x1 + s*x2;
 417          X[2*i+2] = c*x2 - s*x1;
 418
 419          x1 = X[2*i+1];
 420          x2 = X[2*i+3];
 421          X[2*i+1] = c*x1 + s*x2;
 422          X[2*i+3] = c*x2 - s*x1;
 423       }
 424    }
 425 }
 426
 427 void random_rotation(float *X, int R, int dir)
 428 {
 429    int i;
 430    for (i=0;i<NBANDS-4;i++)
 431    {
 432       //crappy_rfft(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i]), R+i, dir);
 433       //printf ("%f ", norm2(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i])));
 434       //rotate_vect(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i]), R+i, dir);
 435
 436       float theta;
 437       if (qbank[i+1]-qbank[i] < 12)
 438          theta = .25;
 439       else
 440          theta = .5;
 441       exp_rotation(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i]), theta, dir);
 442    }
 443    //printf ("\n");
 444 }
 445
 446
 447 void quant_bank(float *X, float *P, float centre)
 448 {
 449    int i;
 450    for (i=0;i<NBANDS;i++)
 451    {
 452       int q;
 453       /*if (centre < 5)
 454          q =qpulses3[i];
 455       else if (centre < 8)
 456          q =qpulses2[i];
 457       else
 458          q =qpulses[i];*/
 459       q =qpulses[i];
 460       if (q)
 461          alg_quant2(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i]), q, P+qbank[i]*2-1);
 462       else
 463          noise_quant(X+qbank[i]*2-1, 2*(qbank[i+1]-qbank[i]), q, P+qbank[i]*2-1);
 464    }
 465    //FIXME: This is a kludge, even though I don't think it really matters much
 466    X[255] = 0;
 467 }
 468
 469 void compute_pitch_gain(float *X, float *P, float *gains, float *bank)
 470 {
 471    int i;
 472    float w[256];
 473    for (i=0;i<NBANDS;i++)
 474    {
 475       int j;
 476       for (j=qbank[i];j<qbank[i+1];j++)
 477          w[j] = bank[i];
 478    }
 479
 480
 481    for (i=0;i<PBANDS;i++)
 482    {
 483       float Sxy=0;
 484       float Sxx = 0;
 485       int j;
 486       float gain;
 487       for (j=pbank[i];j<pbank[i+1];j++)
 488       {
 489          Sxy += X[j*2-1]*P[j*2-1]*w[j];
 490          Sxy += X[j*2]*P[j*2]*w[j];
 491          Sxx += X[j*2-1]*X[j*2-1]*w[j] + X[j*2]*X[j*2]*w[j];
 492       }
 493       gain = Sxy/(1e-10+Sxx);
 494       //gain = Sxy/(2*(pbank[i+1]-pbank[i]));
 495       //if (i<3)
 496       //gain *= 1+.02*gain;
 497       if (gain > .90)
 498          gain = .90;
 499       if (gain < 0.0)
 500          gain = 0.0;
 501
 502       gains[i] = gain;
 503    }
 504    for (i=pbank[PBANDS]*2-1;i<256;i++)
 505       P[i] = 0;
 506    /*if (rand()%10 == 0)
 507    {
 508       for (i=0;i<PBANDS;i++)
 509          printf ("%f ", gains[i]*gains[i]);
 510       printf ("\n");
 511    }*/
 512 }
 513
 514 void pitch_quant_bank(float *X, float *P, float *gains)
 515 {
 516    int i;
 517    for (i=0;i<PBANDS;i++)
 518    {
 519       int j;
 520       for (j=pbank[i];j<pbank[i+1];j++)
 521       {
 522          P[j*2-1] *= gains[i];
 523          P[j*2] *= gains[i];
 524       }
 525       //printf ("%f ", gain);
 526    }
 527    for (i=pbank[PBANDS]*2-1;i<256;i++)
 528       P[i] = 0;
 529 }
 530
 531
 532 void pitch_renormalise_bank(float *X, float *P)
 533 {
 534    int i;
 535    for (i=0;i<NBANDS;i++)
 536    {
 537       int j;
 538       float Rpp=0;
 539       float Rxp=0;
 540       float Rxx=0;
 541       float gain1;
 542       for (j=qbank[i];j<qbank[i+1];j++)
 543       {
 544          Rxp += X[j*2-1]*P[j*2-1];
 545          Rxp += X[j*2  ]*P[j*2  ];
 546          Rpp += P[j*2-1]*P[j*2-1];
 547          Rpp += P[j*2  ]*P[j*2  ];
 548          Rxx += X[j*2-1]*X[j*2-1];
 549          Rxx += X[j*2  ]*X[j*2  ];
 550       }
 551       //Rxx *= .5/(qbank[i+1]-qbank[i]);
 552       //Rxp *= .5/(qbank[i+1]-qbank[i]);
 553       //Rpp *= .5/(qbank[i+1]-qbank[i]);
 554       float arg = Rxp*Rxp + 1 - Rpp;
 555       if (arg < 0)
 556       {
 557          printf ("arg: %f %f %f %f\n", arg, Rxp, Rpp, Rxx);
 558          arg = 0;
 559       }
 560       gain1 = sqrt(arg)-Rxp;
 561       if (Rpp>.9999)
 562          Rpp = .9999;
 563       //gain1 = sqrt(1.-Rpp);
 564
 565       //gain2 = -sqrt(Rxp*Rxp + 1 - Rpp)-Rxp;
 566       //if (fabs(gain2)<fabs(gain1))
 567       //   gain1 = gain2;
 568       //printf ("%f ", Rxx, Rxp, Rpp, gain);
 569       //printf ("%f ", gain1);
 570       Rxx = 0;
 571       for (j=qbank[i];j<qbank[i+1];j++)
 572       {
 573          X[j*2-1] = P[j*2-1]+gain1*X[j*2-1];
 574          X[j*2  ] = P[j*2  ]+gain1*X[j*2  ];
 575          Rxx += X[j*2-1]*X[j*2-1];
 576          Rxx += X[j*2  ]*X[j*2  ];
 577       }
 578       //printf ("%f %f ", gain1, Rxx);
 579    }
 580    //printf ("\n");
 581    //FIXME: Kludge
 582    X[255] = 0;
 583 }
 584
 585
 586 struct CEFTState_ {
 587    void *frame_fft;
 588    int length;
 589 };
 590
 591 CEFTState *ceft_init(int len)
 592 {
 593    CEFTState *st = malloc(sizeof(CEFTState));
 594    st->length = len;
 595    st->frame_fft = spx_fft_init(st->length);
 596    return st;
 597 }
 598
 599
 600
 601 void ceft_encode(CEFTState *st, float *in, float *out, float *pitch, float *window)
 602 {
 603    //float bark[BARK_BANDS];
 604    float X[st->length];
 605    float Xbak[st->length];
 606    float Xp[st->length];
 607    int i;
 608    float bank[NBANDS];
 609    float pitch_bank[NBANDS];
 610    float p[st->length];
 611    float gains[PBANDS];
 612    float mask[NBANDS];
 613    static float obank[NBANDS+1];
 614    spx_fft_float(st->frame_fft, in, X);
 615
 616    /* Bands for the input signal */
 617    compute_bank(X, bank);
 618
 619    for (i=0;i<st->length;i++)
 620       p[i] = pitch[i]*window[i];
 621
 622    spx_fft_float(st->frame_fft, p, Xp);
 623
 624    /* Bands for the pitch signal */
 625    compute_bank(Xp, pitch_bank);
 626 #if 0
 627    if (rand()%10 ==0 && fabs(X[0]) > 2 && (fabs(X[0]) > 10 || rand() % 5 == 0))
 628    {
 629       /*printf ("%f ", 20*log10(1+fabs(X[0])));
 630       for (i=0;i<NBANDS;i++)
 631          printf ("%f ", 20*log10(1+bank[i]));
 632       for (i=0;i<NBANDS+1;i++)
 633          printf ("%f ", obank[i]);
 634       printf ("%f ", 20*log10(1+fabs(Xp[0])));
 635       for (i=0;i<NBANDS;i++)
 636          printf ("%f ", 20*log10(1+pitch_bank[i]));
 637       printf (" \n");*/
 638       for (i=0;i<NBANDS;i++)
 639          printf ("%f ", 20*log10(1+bank[i])-.9*obank[i+1]);
 640       /*printf ("%f ", 20*log10(1+fabs(X[0])));
 641       for (i=0;i<NBANDS;i++)
 642          printf ("%f ", 20*log10(1+bank[i]));*/
 643       printf ("\n");
 644
 645    }
 646    obank[0] = 20*log10(1+fabs(X[0]));
 647    for (i=0;i<NBANDS;i++)
 648       obank[i+1] = 20*log10(1+bank[i]);
 649    return;
 650 #endif
 651
 652 #if 0
 653    float tmp = 1.+X[0]*X[0];
 654    for (i=0;i<NBANDS;i++)
 655    {
 656       tmp = 1. + .1*tmp + bank[i]*bank[i];
 657       mask[i] = bank[i]*bank[i]/tmp;
 658       printf ("%f ", 10*log10(mask[i]));
 659    }
 660    printf ("\n");
 661 #endif
 662
 663    float centre=0;
 664 #if 0
 665    float Sxw = 0, Sw = 0;
 666    for (i=0;i<NBANDS;i++)
 667    {
 668       printf ("%f ", 20*log10(bank[i]));
 669       Sxw += i*sqrt(bank[i]);
 670       Sw  += sqrt(bank[i]);
 671    }
 672    centre = Sxw/Sw;
 673             //printf ("%f\n", Sxw/Sw);
 674 #endif
 675
 676
 677    normalise_bank(X, bank);
 678
 679    float in_bank[NBANDS];
 680    float qbank[NBANDS];
 681    static float last_err[NBANDS];
 682    static float last_bank[NBANDS];
 683
 684    for (i=0;i<NBANDS;i++)
 685    {
 686       in_bank[i] = 20*log10(bank[i]+1) + .0*last_err[i+1] - .9*last_bank[i];
 687    }
 688    for (i=0;i<NBANDS;i++)
 689       qbank[i] = in_bank[i];
 690
 691    quantise_bands(in_bank, qbank, NBANDS);
 692
 693 #if 0
 694    float q = .25f;
 695    for (i=0;i<NBANDS+1;i++)
 696    {
 697       if (i<1)
 698          q = 1.;
 699       else if (i<4)
 700          q = 2.;
 701       else if (i<5)
 702          q = 3.;
 703       else if (i<10)
 704          q = 4.;
 705       else if (i<14)
 706          q = 5.;
 707       else
 708          q = 6.;
 709       int sc = floor(.5 + (in_bank[i]-qbank[i])/q);
 710       printf ("%d ", sc);
 711       qbank[i] = q * sc;
 712    }
 713    printf ("\n");
 714 #endif
 715
 716    /*for (i=0;i<NBANDS+1;i++)
 717       printf ("%f ", in_bank[i]-qbank[i]);
 718    printf ("\n");*/
 719
 720
 721    for (i=0;i<NBANDS;i++)
 722       last_err[i] = qbank[i]-in_bank[i];
 723
 724    for (i=0;i<NBANDS;i++)
 725    {
 726       qbank[i] += .9*last_bank[i];
 727       if (qbank[i]<0)
 728          qbank[i] = 0;
 729       bank[i] = pow(10,(qbank[i])/20)-1;
 730       if (bank[i] < .1)
 731          bank[i] = .1;
 732    }
 733    for (i=0;i<NBANDS;i++)
 734       last_bank[i] = qbank[i];
 735
 736    {
 737       float sign;
 738       int id;
 739       float q = 2.;
 740       if (X[0]<0)
 741          sign = -1;
 742       else
 743          sign = 1;
 744       id = floor(.5+20/q*log10(1+fabs(X[0])));
 745       if (id < 0)
 746          id = 0;
 747       if (id > 32)
 748       {
 749          printf("%d %f\n", id, X[0]);
 750          id = 32;
 751       }
 752       //printf ("%d %f ", id, X[0]);
 753       X[0] = sign*pow(10,(q*id)/20)-1;
 754       //printf ("%f\n", X[0]);
 755    }
 756
 757    normalise_bank(Xp, pitch_bank);
 758
 759    /*
 760    for(i=0;i<st->length;i++)
 761       Xbak[i] = X[i];
 762    for(i=0;i<st->length;i++)
 763       printf ("%f ", X[i]);
 764    printf ("\n");
 765    */
 766
 767    random_rotation(X, 10, -1);
 768    random_rotation(Xp, 10, -1);
 769
 770    compute_pitch_gain(X, Xp, gains, bank);
 771    quantise_pitch(gains, PBANDS);
 772    pitch_quant_bank(X, Xp, gains);
 773
 774    for (i=1;i<st->length;i++)
 775       X[i] -= Xp[i];
 776
 777    //Quantise input
 778    quant_bank(X, Xp, centre);
 779
 780    random_rotation(X, 10, 1);
 781    //pitch_renormalise_bank(X, Xp);
 782
 783 #if 0
 784    float err = 0;
 785    for(i=1;i<19;i++)
 786       err += (Xbak[i] - X[i])*(Xbak[i] - X[i]);
 787    printf ("%f\n", err);
 788 #endif
 789    //printf ("\n");
 790
 791    /* Denormalise back to real power */
 792    denormalise_bank(X, bank);
 793
 794    //Synthesis
 795    spx_ifft_float(st->frame_fft, X, out);
 796
 797 }
 798