Working on some stability issues (appears to be solved by making the pitch
[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 expo
39 nential */
40 static void exp_rotation(float *X, int len, float theta, int dir)
41 {
42    int i;
43    float c, s;
44    c = cos(theta);
45    s = sin(theta);
46    if (dir > 0)
47    {
48       for (i=0;i<(len/2)-1;i++)
49       {
50          float x1, x2;
51          x1 = X[2*i];
52          x2 = X[2*i+2];
53          X[2*i] = c*x1 - s*x2;
54          X[2*i+2] = c*x2 + s*x1;
55          
56          x1 = X[2*i+1];
57          x2 = X[2*i+3];
58          X[2*i+1] = c*x1 - s*x2;
59          X[2*i+3] = c*x2 + s*x1;
60       }
61       for (i=(len/2)-3;i>=0;i--)
62       {
63          float x1, x2;
64          x1 = X[2*i];
65          x2 = X[2*i+2];
66          X[2*i] = c*x1 - s*x2;
67          X[2*i+2] = c*x2 + s*x1;
68          
69          x1 = X[2*i+1];
70          x2 = X[2*i+3];
71          X[2*i+1] = c*x1 - s*x2;
72          X[2*i+3] = c*x2 + s*x1;
73       }
74
75    } else {
76       for (i=0;i<(len/2)-2;i++)
77       {
78          float x1, x2;
79          x1 = X[2*i];
80          x2 = X[2*i+2];
81          X[2*i] = c*x1 + s*x2;
82          X[2*i+2] = c*x2 - s*x1;
83          
84          x1 = X[2*i+1];
85          x2 = X[2*i+3];
86          X[2*i+1] = c*x1 + s*x2;
87          X[2*i+3] = c*x2 - s*x1;
88       }
89       
90       for (i=(len/2)-2;i>=0;i--)
91       {
92          float x1, x2;
93          x1 = X[2*i];
94          x2 = X[2*i+2];
95          X[2*i] = c*x1 + s*x2;
96          X[2*i+2] = c*x2 - s*x1;
97          
98          x1 = X[2*i+1];
99          x2 = X[2*i+3];
100          X[2*i+1] = c*x1 + s*x2;
101          X[2*i+3] = c*x2 - s*x1;
102       }
103    }
104 }
105
106
107 /* Compute the energy in each of the bands */
108 void compute_band_energies(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       bank[i] = 1e-10;
117       for (j=B*eBands[i];j<B*eBands[i+1];j++)
118          bank[i] += X[j]*X[j];
119       bank[i] = sqrt(bank[i]);
120    }
121 }
122
123 /* Normalise each band such that the energy is one. */
124 void normalise_bands(const CELTMode *m, float *X, float *bank)
125 {
126    int i, B;
127    const int *eBands = m->eBands;
128    B = m->nbMdctBlocks*m->nbChannels;
129    for (i=0;i<m->nbEBands;i++)
130    {
131       int j;
132       float x = 1.f/(1e-10+bank[i]);
133       for (j=B*eBands[i];j<B*eBands[i+1];j++)
134          X[j] *= x;
135    }
136    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
137       X[i] = 0;
138 }
139
140 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
141 void denormalise_bands(const CELTMode *m, float *X, float *bank)
142 {
143    int i, B;
144    const int *eBands = m->eBands;
145    B = m->nbMdctBlocks*m->nbChannels;
146    for (i=0;i<m->nbEBands;i++)
147    {
148       int j;
149       float x = bank[i];
150       for (j=B*eBands[i];j<B*eBands[i+1];j++)
151          X[j] *= x;
152    }
153    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
154       X[i] = 0;
155 }
156
157
158 /* Compute the best gain for each "pitch band" */
159 void compute_pitch_gain(const CELTMode *m, float *X, float *P, float *gains, float *bank)
160 {
161    int i, B;
162    const int *eBands = m->eBands;
163    const int *pBands = m->pBands;
164    B = m->nbMdctBlocks*m->nbChannels;
165    float w[B*eBands[m->nbEBands]];
166    for (i=0;i<m->nbEBands;i++)
167    {
168       int j;
169       for (j=B*eBands[i];j<B*eBands[i+1];j++)
170          w[j] = bank[i];
171    }
172
173    
174    for (i=0;i<m->nbPBands;i++)
175    {
176       float Sxy=0;
177       float Sxx = 0;
178       int j;
179       float gain;
180       for (j=B*pBands[i];j<B*pBands[i+1];j++)
181       {
182          Sxy += X[j]*P[j]*w[j];
183          Sxx += X[j]*X[j]*w[j];
184       }
185       gain = Sxy/(1e-10+Sxx);
186       //gain = Sxy/(2*(pbank[i+1]-pbank[i]));
187       //if (i<3)
188       //gain *= 1+.02*gain;
189       if (gain > .90)
190          gain = .90;
191       if (gain < 0.0)
192          gain = 0.0;
193       
194       gains[i] = gain;
195    }
196    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
197       P[i] = 0;
198 }
199
200 /* Apply the (quantised) gain to each "pitch band" */
201 void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains)
202 {
203    int i, B;
204    const int *pBands = m->pBands;
205    B = m->nbMdctBlocks*m->nbChannels;
206    for (i=0;i<m->nbPBands;i++)
207    {
208       int j;
209       for (j=B*pBands[i];j<B*pBands[i+1];j++)
210          P[j] *= gains[i];
211       //printf ("%f ", gain);
212    }
213    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
214       P[i] = 0;
215 }
216
217 void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc)
218 {
219    int i, j, B;
220    const int *eBands = m->eBands;
221    B = m->nbMdctBlocks*m->nbChannels;
222    float norm[B*eBands[m->nbEBands+1]];
223    
224    for (i=0;i<m->nbEBands;i++)
225    {
226       int q;
227       q = m->nbPulses[i];
228       if (q>0) {
229          float n = sqrt(B*(eBands[i+1]-eBands[i]));
230          alg_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
231          for (j=B*eBands[i];j<B*eBands[i+1];j++)
232             norm[j] = X[j] * n;
233          //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
234          //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
235       } else {
236          float n = sqrt(B*(eBands[i+1]-eBands[i]));
237          copy_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), -q, norm, B, eBands[i], enc);
238          for (j=B*eBands[i];j<B*eBands[i+1];j++)
239             norm[j] = X[j] * n;
240          //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs64(B*(eBands[i+1]-eBands[i]), -q)))/(B*(eBands[i+1]-eBands[i])));
241          //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs64(B*(eBands[i+1]-eBands[i]), -q))));
242       }
243    }
244    //printf ("\n");
245    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
246       X[i] = 0;
247 }
248
249 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
250 {
251    int i, j, B;
252    const int *eBands = m->eBands;
253    B = m->nbMdctBlocks*m->nbChannels;
254    float norm[B*eBands[m->nbEBands+1]];
255    
256    for (i=0;i<m->nbEBands;i++)
257    {
258       int q;
259       q = m->nbPulses[i];
260       if (q>0) {
261          float n = sqrt(B*(eBands[i+1]-eBands[i]));
262          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, dec);
263          for (j=B*eBands[i];j<B*eBands[i+1];j++)
264             norm[j] = X[j] * n;
265       } else {
266          float n = sqrt(B*(eBands[i+1]-eBands[i]));
267          for (j=B*eBands[i];j<B*eBands[i+1];j++)
268             X[j] = 0;
269          copy_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), -q, norm, B, eBands[i], dec);
270          for (j=B*eBands[i];j<B*eBands[i+1];j++)
271             norm[j] = X[j] * n;
272       }
273    }
274    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
275       X[i] = 0;
276 }
277
278 void band_rotation(const CELTMode *m, float *X, int dir)
279 {
280    int i, B;
281    const int *eBands = m->eBands;
282    B = m->nbMdctBlocks*m->nbChannels;
283    for (i=0;i<m->nbEBands;i++)
284    {
285       float theta;
286       theta = pow(.1f,1.f*abs(m->nbPulses[i])/(B*(eBands[i+1]-eBands[i])));
287       exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, dir);
288    }
289    //printf ("\n");
290 }