Some work on index packing
[opus.git] / libcelt / vq.c
1 /* (C) 2007 Jean-Marc Valin, CSIRO
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #include <math.h>
33 #include <stdlib.h>
34 #include "cwrs.h"
35
36 /* Algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
37    a combination of pulses such that its norm is still equal to 1 */
38 void alg_quant(float *x, int N, int K, float *p)
39 {
40    float y[N];
41    int i,j;
42    float xy = 0;
43    float yy = 0;
44    float yp = 0;
45    float Rpp=0;
46    float gain=0;
47    for (j=0;j<N;j++)
48       Rpp += p[j]*p[j];
49    for (i=0;i<N;i++)
50       y[i] = 0;
51       
52    for (i=0;i<K;i++)
53    {
54       int best_id=0;
55       float max_val=-1e10;
56       float best_xy=0, best_yy=0, best_yp = 0;
57       for (j=0;j<N;j++)
58       {
59          float tmp_xy, tmp_yy, tmp_yp;
60          float score;
61          float g;
62          tmp_xy = xy + fabs(x[j]);
63          tmp_yy = yy + 2*fabs(y[j]) + 1;
64          if (x[j]>0)
65             tmp_yp = yp + p[j];
66          else
67             tmp_yp = yp - p[j];
68          g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
69          score = 2*g*tmp_xy - g*g*tmp_yy;
70          if (score>max_val)
71          {
72             max_val = score;
73             best_id = j;
74             best_xy = tmp_xy;
75             best_yy = tmp_yy;
76             best_yp = tmp_yp;
77             gain = g;
78          }
79       }
80
81       xy = best_xy;
82       yy = best_yy;
83       yp = best_yp;
84       if (x[best_id]>0)
85          y[best_id] += 1;
86       else
87          y[best_id] -= 1;
88    }
89    
90    for (i=0;i<N;i++)
91       x[i] = p[i]+gain*y[i];
92    
93 }
94
95 /* Improved algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
96    a combination of pulses such that its norm is still equal to 1. The only difference with 
97    the quantiser above is that the search is more complete. */
98 int alg_quant2(float *x, int N, int K, float *p)
99 {
100    int L = 5;
101    //float tata[200];
102    float y[L][N];
103    int iy[L][N];
104    //float tata2[200];
105    float ny[L][N];
106    int iny[L][N];
107    int i, j, m;
108    float xy[L], nxy[L];
109    float yy[L], nyy[L];
110    float yp[L], nyp[L];
111    float best_scores[L];
112    float Rpp=0, Rxp=0;
113    float gain[L];
114    int maxL = 1;
115    float alpha = .9;
116    
117    for (j=0;j<N;j++)
118       Rpp += p[j]*p[j];
119    for (j=0;j<N;j++)
120       Rxp += x[j]*p[j];
121    for (m=0;m<L;m++)
122       for (i=0;i<N;i++)
123          y[m][i] = 0;
124       
125    for (m=0;m<L;m++)
126       for (i=0;i<N;i++)
127          ny[m][i] = 0;
128
129    for (m=0;m<L;m++)
130       for (i=0;i<N;i++)
131          iy[m][i] = iny[m][i] = 0;
132
133    for (m=0;m<L;m++)
134       xy[m] = yy[m] = yp[m] = gain[m] = 0;
135    
136    for (i=0;i<K;i++)
137    {
138       int L2 = L;
139       if (L>maxL)
140       {
141          L2 = maxL;
142          maxL *= N;
143       }
144       for (m=0;m<L;m++)
145          best_scores[m] = -1e10;
146
147       for (m=0;m<L2;m++)
148       {
149          for (j=0;j<N;j++)
150          {
151             int sign;
152             for (sign=-1;sign<=1;sign+=2)
153             {
154                if (iy[m][j]*sign < 0)
155                   continue;
156                //fprintf (stderr, "%d/%d %d/%d %d/%d\n", i, K, m, L2, j, N);
157                float tmp_xy, tmp_yy, tmp_yp;
158                float score;
159                float g;
160                float s = sign;
161                tmp_xy = xy[m] + s*x[j]               - alpha*s*p[j]*Rxp;
162                tmp_yy = yy[m] + 2*s*y[m][j] + 1      +alpha*alpha*p[j]*p[j]*Rpp - 2*alpha*s*p[j]*yp[m] - 2*alpha*p[j]*p[j];
163                tmp_yp = yp[m] + s*p[j]               *(1-alpha*Rpp);
164                g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
165                score = 2*g*tmp_xy - g*g*tmp_yy;
166
167                if (score>best_scores[L-1])
168                {
169                   int k, n;
170                   int id = L-1;
171                   while (id > 0 && score > best_scores[id-1])
172                      id--;
173                
174                   for (k=L-1;k>id;k--)
175                   {
176                      nxy[k] = nxy[k-1];
177                      nyy[k] = nyy[k-1];
178                      nyp[k] = nyp[k-1];
179                      //fprintf(stderr, "%d %d \n", N, k);
180                      for (n=0;n<N;n++)
181                         ny[k][n] = ny[k-1][n];
182                      for (n=0;n<N;n++)
183                         iny[k][n] = iny[k-1][n];
184                      gain[k] = gain[k-1];
185                      best_scores[k] = best_scores[k-1];
186                   }
187
188                   nxy[id] = tmp_xy;
189                   nyy[id] = tmp_yy;
190                   nyp[id] = tmp_yp;
191                   gain[id] = g;
192                   for (n=0;n<N;n++)
193                      ny[id][n] = y[m][n];
194                   ny[id][j] += s;
195                   for (n=0;n<N;n++)
196                      ny[id][n] -= alpha*s*p[j]*p[n];
197                
198                   for (n=0;n<N;n++)
199                      iny[id][n] = iy[m][n];
200                   if (s>0)
201                      iny[id][j] += 1;
202                   else
203                      iny[id][j] -= 1;
204                   best_scores[id] = score;
205                }
206             }   
207          }
208          
209       }
210       int k,n;
211       for (k=0;k<L;k++)
212       {
213          xy[k] = nxy[k];
214          yy[k] = nyy[k];
215          yp[k] = nyp[k];
216          for (n=0;n<N;n++)
217             y[k][n] = ny[k][n];
218          for (n=0;n<N;n++)
219             iy[k][n] = iny[k][n];
220       }
221
222    }
223    
224    for (i=0;i<N;i++)
225       x[i] = p[i]+gain[0]*y[0][i];
226    if (0) {
227       float E=1e-15;
228       int ABS = 0;
229       for (i=0;i<N;i++)
230          ABS += abs(iy[0][i]);
231       //if (K != ABS)
232       //   printf ("%d %d\n", K, ABS);
233       for (i=0;i<N;i++)
234          E += x[i]*x[i];
235       //printf ("%f\n", E);
236       E = 1/sqrt(E);
237       for (i=0;i<N;i++)
238          x[i] *= E;
239    }
240    int comb[K];
241    int signs[K];
242    pulse2comb(N, K, comb, signs, iy[0]); 
243    return icwrs(N, K, comb, signs);
244 }
245
246 /* Just replace the band with noise of unit energy */
247 void noise_quant(float *x, int N, int K, float *p)
248 {
249    int i;
250    float E = 1e-10;
251    for (i=0;i<N;i++)
252    {
253       x[i] = (rand()%1000)/500.-1;
254       E += x[i]*x[i];
255    }
256    E = 1./sqrt(E);
257    for (i=0;i<N;i++)
258    {
259       x[i] *= E;
260    }
261 }
262
263 static const float pg[5] = {1.f, .82f, .75f, 0.7f, 0.6f};
264
265 /* Finds the right offset into Y and copy it */
266 void copy_quant(float *x, int N, int K, float *Y, int B, int N0)
267 {
268    int i,j;
269    int best=0;
270    float best_score=0;
271    float s = 1;
272    float E;
273    for (i=0;i<N0*B-N;i+=B)
274    {
275       int j;
276       float xy=0, yy=0;
277       float score;
278       for (j=0;j<N;j++)
279       {
280          xy += x[j]*Y[i+j];
281          yy += Y[i+j]*Y[i+j];
282       }
283       score = xy*xy/(.001+yy);
284       if (score > best_score)
285       {
286          best_score = score;
287          best = i;
288          if (xy>0)
289             s = 1;
290          else
291             s = -1;
292       }
293    }
294    //printf ("%d %f\n", best, best_score);
295    if (K==0)
296    {
297       E = 1e-10;
298       for (j=0;j<N;j++)
299       {
300          x[j] = s*Y[best+j];
301          E += x[j]*x[j];
302       }
303       E = 1/sqrt(E);
304       for (j=0;j<N;j++)
305          x[j] *= E;
306    } else {
307       float P[N];
308       float pred_gain;
309       if (K>4)
310          pred_gain = .5;
311       else
312          pred_gain = pg[K];
313       E = 1e-10;
314       for (j=0;j<N;j++)
315       {
316          P[j] = s*Y[best+j];
317          E += P[j]*P[j];
318       }
319       E = .8/sqrt(E);
320       for (j=0;j<N;j++)
321          P[j] *= E;
322       alg_quant2(x, N, K, P);
323    }
324 }