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