Comments/cleanup, no code change
[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 > .90)
174          gain = .90;
175       if (gain < 0.0)
176          gain = 0.0;
177       
178       gains[i] = gain;
179    }
180    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
181       P[i] = 0;
182 }
183
184 /* Apply the (quantised) gain to each "pitch band" */
185 void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains)
186 {
187    int i, B;
188    const int *pBands = m->pBands;
189    B = m->nbMdctBlocks*m->nbChannels;
190    for (i=0;i<m->nbPBands;i++)
191    {
192       int j;
193       for (j=B*pBands[i];j<B*pBands[i+1];j++)
194          P[j] *= gains[i];
195       //printf ("%f ", gain);
196    }
197    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
198       P[i] = 0;
199 }
200
201 /* Quantisation of the residual */
202 void quant_bands(const CELTMode *m, float *X, float *P, float *W, ec_enc *enc)
203 {
204    int i, j, B;
205    const int *eBands = m->eBands;
206    B = m->nbMdctBlocks*m->nbChannels;
207    float norm[B*eBands[m->nbEBands+1]];
208    
209    for (i=0;i<m->nbEBands;i++)
210    {
211       int q;
212       float theta, n;
213       q = m->nbPulses[i];
214       n = sqrt(B*(eBands[i+1]-eBands[i]));
215       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
216          
217       if (q<=0) {
218          q = -q;
219          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);
220       }
221          
222       if (q != 0)
223       {
224          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
225          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
226          alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
227          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
228       }
229       for (j=B*eBands[i];j<B*eBands[i+1];j++)
230          norm[j] = X[j] * n;
231       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
232       //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
233    }
234    //printf ("\n");
235    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
236       X[i] = 0;
237 }
238
239 /* Decoding of the residual */
240 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec)
241 {
242    int i, j, B;
243    const int *eBands = m->eBands;
244    B = m->nbMdctBlocks*m->nbChannels;
245    float norm[B*eBands[m->nbEBands+1]];
246    
247    for (i=0;i<m->nbEBands;i++)
248    {
249       int q;
250       float theta, n;
251       q = m->nbPulses[i];
252       n = sqrt(B*(eBands[i+1]-eBands[i]));
253       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+abs(m->nbPulses[i]));
254
255       if (q<=0) {
256          q = -q;
257          intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
258       }
259
260       if (q != 0)
261       {
262          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
263          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, dec);
264          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
265       }
266       for (j=B*eBands[i];j<B*eBands[i+1];j++)
267          norm[j] = X[j] * n;
268    }
269    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
270       X[i] = 0;
271 }