Nothing to see here.
[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
35 /* Algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
36    a combination of pulses such that its norm is still equal to 1 */
37 void alg_quant(float *x, int N, int K, float *p)
38 {
39    float y[N];
40    int i,j;
41    float xy = 0;
42    float yy = 0;
43    float yp = 0;
44    float Rpp=0;
45    float gain=0;
46    for (j=0;j<N;j++)
47       Rpp += p[j]*p[j];
48    for (i=0;i<N;i++)
49       y[i] = 0;
50       
51    for (i=0;i<K;i++)
52    {
53       int best_id=0;
54       float max_val=-1e10;
55       float best_xy=0, best_yy=0, best_yp = 0;
56       for (j=0;j<N;j++)
57       {
58          float tmp_xy, tmp_yy, tmp_yp;
59          float score;
60          float g;
61          tmp_xy = xy + fabs(x[j]);
62          tmp_yy = yy + 2*fabs(y[j]) + 1;
63          if (x[j]>0)
64             tmp_yp = yp + p[j];
65          else
66             tmp_yp = yp - p[j];
67          g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
68          score = 2*g*tmp_xy - g*g*tmp_yy;
69          if (score>max_val)
70          {
71             max_val = score;
72             best_id = j;
73             best_xy = tmp_xy;
74             best_yy = tmp_yy;
75             best_yp = tmp_yp;
76             gain = g;
77          }
78       }
79
80       xy = best_xy;
81       yy = best_yy;
82       yp = best_yp;
83       if (x[best_id]>0)
84          y[best_id] += 1;
85       else
86          y[best_id] -= 1;
87    }
88    
89    for (i=0;i<N;i++)
90       x[i] = p[i]+gain*y[i];
91    
92 }
93
94 /* Improved algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
95    a combination of pulses such that its norm is still equal to 1. The only difference with 
96    the quantiser above is that the search is more complete. */
97 void alg_quant2(float *x, int N, int K, float *p)
98 {
99    int L = 5;
100    //float tata[200];
101    float y[L][N];
102    //float tata2[200];
103    float ny[L][N];
104    int i, j, m;
105    float xy[L], nxy[L];
106    float yy[L], nyy[L];
107    float yp[L], nyp[L];
108    float best_scores[L];
109    float Rpp=0;
110    float gain[L];
111    int maxL = 1;
112    for (j=0;j<N;j++)
113       Rpp += p[j]*p[j];
114    for (m=0;m<L;m++)
115       for (i=0;i<N;i++)
116          y[m][i] = 0;
117       
118    for (m=0;m<L;m++)
119       for (i=0;i<N;i++)
120          ny[m][i] = 0;
121
122    for (m=0;m<L;m++)
123       xy[m] = yy[m] = yp[m] = gain[m] = 0;
124    
125    for (i=0;i<K;i++)
126    {
127       int L2 = L;
128       if (L>maxL)
129       {
130          L2 = maxL;
131          maxL *= N;
132       }
133       for (m=0;m<L;m++)
134          best_scores[m] = -1e10;
135
136       for (m=0;m<L2;m++)
137       {
138          for (j=0;j<N;j++)
139          {
140             //fprintf (stderr, "%d/%d %d/%d %d/%d\n", i, K, m, L2, j, N);
141             float tmp_xy, tmp_yy, tmp_yp;
142             float score;
143             float g;
144             tmp_xy = xy[m] + fabs(x[j]);
145             tmp_yy = yy[m] + 2*fabs(y[m][j]) + 1;
146             if (x[j]>0)
147                tmp_yp = yp[m] + p[j];
148             else
149                tmp_yp = yp[m] - p[j];
150             g = (sqrt(tmp_yp*tmp_yp + tmp_yy - tmp_yy*Rpp) - tmp_yp)/tmp_yy;
151             score = 2*g*tmp_xy - g*g*tmp_yy;
152
153             if (score>best_scores[L-1])
154             {
155                int k, n;
156                int id = L-1;
157                while (id > 0 && score > best_scores[id-1])
158                   id--;
159                
160                for (k=L-1;k>id;k--)
161                {
162                   nxy[k] = nxy[k-1];
163                   nyy[k] = nyy[k-1];
164                   nyp[k] = nyp[k-1];
165                   //fprintf(stderr, "%d %d \n", N, k);
166                   for (n=0;n<N;n++)
167                      ny[k][n] = ny[k-1][n];
168                   gain[k] = gain[k-1];
169                   best_scores[k] = best_scores[k-1];
170                }
171
172                nxy[id] = tmp_xy;
173                nyy[id] = tmp_yy;
174                nyp[id] = tmp_yp;
175                gain[id] = g;
176                for (n=0;n<N;n++)
177                   ny[id][n] = y[m][n];
178                if (x[j]>0)
179                   ny[id][j] += 1;
180                else
181                   ny[id][j] -= 1;
182                best_scores[id] = score;
183             }
184             
185          }
186          
187       }
188       int k,n;
189       for (k=0;k<L;k++)
190       {
191          xy[k] = nxy[k];
192          yy[k] = nyy[k];
193          yp[k] = nyp[k];
194          for (n=0;n<N;n++)
195             y[k][n] = ny[k][n];
196       }
197
198    }
199    
200    for (i=0;i<N;i++)
201       x[i] = p[i]+gain[0]*y[0][i];
202    
203 }
204
205 /* Just replace the band with noise of unit energy */
206 void noise_quant(float *x, int N, int K, float *p)
207 {
208    int i;
209    float E = 1e-10;
210    for (i=0;i<N;i++)
211    {
212       x[i] = (rand()%1000)/500.-1;
213       E += x[i]*x[i];
214    }
215    E = 1./sqrt(E);
216    for (i=0;i<N;i++)
217    {
218       x[i] *= E;
219    }
220 }