Now using 8 bands for the pitch gain, with a 128-entry codebook.
[opus.git] / libcelt / bands.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 "bands.h"
34 #include "modes.h"
35 #include "vq.h"
36 #include "cwrs.h"
37
38 /* Applies a series of rotations so that pulses are spread like a two-sided
39 exponential. The effect of this is to reduce the tonal noise created by the
40 sparse spectrum resulting from the pulse codebook */
41 static void exp_rotation(float *X, int len, float theta, int dir, int stride, int iter)
42 {
43    int i, k;
44    float c, s;
45    c = cos(theta);
46    s = sin(theta);
47    if (dir > 0)
48    {
49       for (k=0;k<iter;k++)
50       {
51          for (i=0;i<len-stride;i++)
52          {
53             float x1, x2;
54             x1 = X[i];
55             x2 = X[i+stride];
56             X[i] = c*x1 - s*x2;
57             X[i+stride] = c*x2 + s*x1;
58          }
59          for (i=len-2*stride-1;i>=0;i--)
60          {
61             float x1, x2;
62             x1 = X[i];
63             x2 = X[i+stride];
64             X[i] = c*x1 - s*x2;
65             X[i+stride] = c*x2 + s*x1;
66          }
67       }
68    } else {
69       for (k=0;k<iter;k++)
70       {
71          for (i=0;i<len-2*stride;i++)
72          {
73             float x1, x2;
74             x1 = X[i];
75             x2 = X[i+stride];
76             X[i] = c*x1 + s*x2;
77             X[i+stride] = c*x2 - s*x1;
78          }
79          for (i=len-stride-1;i>=0;i--)
80          {
81             float x1, x2;
82             x1 = X[i];
83             x2 = X[i+stride];
84             X[i] = c*x1 + s*x2;
85             X[i+stride] = c*x2 - s*x1;
86          }
87       }
88    }
89 }
90
91 /* Compute the amplitude (sqrt energy) in each of the bands */
92 void compute_band_energies(const CELTMode *m, float *X, float *bank)
93 {
94    int i, B;
95    const int *eBands = m->eBands;
96    B = m->nbMdctBlocks*m->nbChannels;
97    for (i=0;i<m->nbEBands;i++)
98    {
99       int j;
100       bank[i] = 1e-10;
101       for (j=B*eBands[i];j<B*eBands[i+1];j++)
102          bank[i] += X[j]*X[j];
103       bank[i] = sqrt(bank[i]);
104    }
105 }
106
107 /* Normalise each band such that the energy is one. */
108 void normalise_bands(const CELTMode *m, float *X, float *bank)
109 {
110    int i, B;
111    const int *eBands = m->eBands;
112    B = m->nbMdctBlocks*m->nbChannels;
113    for (i=0;i<m->nbEBands;i++)
114    {
115       int j;
116       float x = 1.f/(1e-10+bank[i]);
117       for (j=B*eBands[i];j<B*eBands[i+1];j++)
118          X[j] *= x;
119    }
120    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
121       X[i] = 0;
122 }
123
124 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
125 void denormalise_bands(const CELTMode *m, float *X, float *bank)
126 {
127    int i, B;
128    const int *eBands = m->eBands;
129    B = m->nbMdctBlocks*m->nbChannels;
130    for (i=0;i<m->nbEBands;i++)
131    {
132       int j;
133       float x = bank[i];
134       for (j=B*eBands[i];j<B*eBands[i+1];j++)
135          X[j] *= x;
136    }
137    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
138       X[i] = 0;
139 }
140
141
142 /* Compute the best gain for each "pitch band" */
143 void compute_pitch_gain(const CELTMode *m, float *X, float *P, float *gains, float *bank)
144 {
145    int i, B;
146    const int *eBands = m->eBands;
147    const int *pBands = m->pBands;
148    B = m->nbMdctBlocks*m->nbChannels;
149    float w[B*eBands[m->nbEBands]];
150    for (i=0;i<m->nbEBands;i++)
151    {
152       int j;
153       for (j=B*eBands[i];j<B*eBands[i+1];j++)
154          w[j] = bank[i];
155    }
156
157    
158    for (i=0;i<m->nbPBands;i++)
159    {
160       float Sxy=0;
161       float Sxx = 0;
162       int j;
163       float gain;
164       for (j=B*pBands[i];j<B*pBands[i+1];j++)
165       {
166          Sxy += X[j]*P[j]*w[j];
167          Sxx += X[j]*X[j]*w[j];
168       }
169       gain = Sxy/(1e-10+Sxx);
170       //gain = Sxy/(2*(pbank[i+1]-pbank[i]));
171       //if (i<3)
172       //gain *= 1+.02*gain;
173       if (gain > 1.f)
174          gain = 1.f;
175       if (gain < 0.0f)
176          gain = 0.0f;
177       /* We need to be a bit conservative, otherwise residual doesn't quantise well */
178       gain *= .9f;
179       gains[i] = gain;
180       //printf ("%f ", 1-sqrt(1-gain*gain));
181    }
182    /*if(rand()%10==0)
183    {
184       for (i=0;i<m->nbPBands;i++)
185          printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
186       printf ("\n");
187    }*/
188    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
189       P[i] = 0;
190 }
191
192 /* Apply the (quantised) gain to each "pitch band" */
193 void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains)
194 {
195    int i, B;
196    const int *pBands = m->pBands;
197    B = m->nbMdctBlocks*m->nbChannels;
198    for (i=0;i<m->nbPBands;i++)
199    {
200       int j;
201       for (j=B*pBands[i];j<B*pBands[i+1];j++)
202          P[j] *= gains[i];
203       //printf ("%f ", gain);
204    }
205    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
206       P[i] = 0;
207 }
208
209 /* Quantisation of the residual */
210 void quant_bands(const CELTMode *m, float *X, float *P, float *W, ec_enc *enc)
211 {
212    int i, j, B;
213    const int *eBands = m->eBands;
214    B = m->nbMdctBlocks*m->nbChannels;
215    float norm[B*eBands[m->nbEBands+1]];
216    
217    for (i=0;i<m->nbEBands;i++)
218    {
219       int q;
220       float theta, n;
221       q = m->nbPulses[i];
222       n = sqrt(B*(eBands[i+1]-eBands[i]));
223       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
224          
225       if (q<=0) {
226          q = -q;
227          intra_prediction(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], enc);
228       }
229          
230       if (q != 0)
231       {
232          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
233          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
234          alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
235          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
236       }
237       for (j=B*eBands[i];j<B*eBands[i+1];j++)
238          norm[j] = X[j] * n;
239       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
240       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
241    }
242    //printf ("\n");
243    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
244       X[i] = 0;
245 }
246
247 /* Decoding of the residual */
248 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
249 {
250    int i, j, B;
251    const int *eBands = m->eBands;
252    B = m->nbMdctBlocks*m->nbChannels;
253    float norm[B*eBands[m->nbEBands+1]];
254    
255    for (i=0;i<m->nbEBands;i++)
256    {
257       int q;
258       float theta, n;
259       q = m->nbPulses[i];
260       n = sqrt(B*(eBands[i+1]-eBands[i]));
261       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
262
263       if (q<=0) {
264          q = -q;
265          intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
266       }
267
268       if (q != 0)
269       {
270          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
271          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, dec);
272          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
273       }
274       for (j=B*eBands[i];j<B*eBands[i+1];j++)
275          norm[j] = X[j] * n;
276    }
277    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
278       X[i] = 0;
279 }