Some work on index packing
[opus.git] / libcelt / celt.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 "os_support.h"
33 #include "mdct.h"
34 #include <math.h>
35 #include "celt.h"
36 #include "pitch.h"
37 #include "fftwrap.h"
38 #include "bands.h"
39 #include "modes.h"
40 #include "probenc.h"
41
42 #define MAX_PERIOD 1024
43
44 struct CELTState_ {
45    const CELTMode *mode;
46    int frame_size;
47    int block_size;
48    int nb_blocks;
49       
50    ec_byte_buffer buf;
51    ec_enc         enc;
52
53    float preemph;
54    float preemph_memE;
55    float preemph_memD;
56    
57    mdct_lookup mdct_lookup;
58    void *fft;
59    
60    float *window;
61    float *in_mem;
62    float *mdct_overlap;
63    float *out_mem;
64
65    float *oldBandE;
66 };
67
68
69
70 CELTState *celt_encoder_new(const CELTMode *mode)
71 {
72    int i, N, B;
73    N = mode->mdctSize;
74    B = mode->nbMdctBlocks;
75    CELTState *st = celt_alloc(sizeof(CELTState));
76    
77    st->mode = mode;
78    st->frame_size = B*N;
79    st->block_size = N;
80    st->nb_blocks  = B;
81    
82    ec_byte_writeinit(&st->buf);
83    ec_enc_init(&st->enc,&st->buf);
84
85    mdct_init(&st->mdct_lookup, 2*N);
86    st->fft = spx_fft_init(MAX_PERIOD);
87    
88    st->window = celt_alloc(2*N*sizeof(float));
89    st->in_mem = celt_alloc(N*sizeof(float));
90    st->mdct_overlap = celt_alloc(N*sizeof(float));
91    st->out_mem = celt_alloc(MAX_PERIOD*sizeof(float));
92    for (i=0;i<N;i++)
93       st->window[i] = st->window[2*N-i-1] = sin(.5*M_PI* sin(.5*M_PI*(i+.5)/N) * sin(.5*M_PI*(i+.5)/N));
94    
95    st->oldBandE = celt_alloc(mode->nbEBands*sizeof(float));
96
97    st->preemph = 0.8;
98    return st;
99 }
100
101 void celt_encoder_destroy(CELTState *st)
102 {
103    if (st == NULL)
104    {
105       celt_warning("NULL passed to celt_encoder_destroy");
106       return;
107    }
108    mdct_clear(&st->mdct_lookup);
109    spx_fft_destroy(st->fft);
110
111    celt_free(st->window);
112    celt_free(st->in_mem);
113    celt_free(st->mdct_overlap);
114    celt_free(st->out_mem);
115    
116    celt_free(st);
117 }
118
119 static void haar1(float *X, int N)
120 {
121    int i;
122    for (i=0;i<N;i+=2)
123    {
124       float a, b;
125       a = X[i];
126       b = X[i+1];
127       X[i] = .707107f*(a+b);
128       X[i+1] = .707107f*(a-b);
129    }
130 }
131
132 static void inv_haar1(float *X, int N)
133 {
134    int i;
135    for (i=0;i<N;i+=2)
136    {
137       float a, b;
138       a = X[i];
139       b = X[i+1];
140       X[i] = .707107f*(a+b);
141       X[i+1] = .707107f*(a-b);
142    }
143 }
144
145 static void compute_mdcts(mdct_lookup *mdct_lookup, float *window, float *in, float *out, int N, int B)
146 {
147    int i;
148    for (i=0;i<B;i++)
149    {
150       int j;
151       float x[2*N];
152       float tmp[N];
153       for (j=0;j<2*N;j++)
154          x[j] = window[j]*in[i*N+j];
155       mdct_forward(mdct_lookup, x, tmp);
156       /* Interleaving the sub-frames */
157       for (j=0;j<N;j++)
158          out[B*j+i] = tmp[j];
159    }
160
161 }
162
163 int celt_encode(CELTState *st, short *pcm)
164 {
165    int i, N, B;
166    N = st->block_size;
167    B = st->nb_blocks;
168    float in[(B+1)*N];
169    
170    float X[B*N];         /**< Interleaved signal MDCTs */
171    float P[B*N];         /**< Interleaved pitch MDCTs*/
172    float bandE[st->mode->nbEBands];
173    float gains[st->mode->nbPBands];
174    int pitch_index;
175    
176    for (i=0;i<N;i++)
177       in[i] = st->in_mem[i];
178    for (;i<(B+1)*N;i++)
179    {
180       float tmp = pcm[i-N];
181       in[i] = tmp - st->preemph*st->preemph_memE;
182       st->preemph_memE = tmp;
183    }
184    for (i=0;i<N;i++)
185       st->in_mem[i] = in[B*N+i];
186
187    /* Compute MDCTs */
188    compute_mdcts(&st->mdct_lookup, st->window, in, X, N, B);
189    
190    /* Pitch analysis */
191    for (i=0;i<N;i++)
192    {
193       in[i] *= st->window[i];
194       in[B*N+i] *= st->window[N+i];
195    }
196    find_spectral_pitch(st->fft, in, st->out_mem, MAX_PERIOD, (B+1)*N, &pitch_index);
197    
198    /* Compute MDCTs of the pitch part */
199    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index, P, N, B);
200    
201    /*int j;
202    for (j=0;j<B*N;j++)
203       printf ("%f ", X[j]);
204    for (j=0;j<B*N;j++)
205       printf ("%f ", P[j]);
206    printf ("\n");*/
207    //haar1(X, B*N);
208    //haar1(P, B*N);
209    
210    /* Band normalisation */
211    compute_band_energies(st->mode, X, bandE);
212    normalise_bands(st->mode, X, bandE);
213    //for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");
214    
215    {
216       float bandEp[st->mode->nbEBands];
217       compute_band_energies(st->mode, P, bandEp);
218       normalise_bands(st->mode, P, bandEp);
219    }
220    
221    quant_energy(st->mode, bandE, st->oldBandE);
222    
223    /* Pitch prediction */
224    compute_pitch_gain(st->mode, X, P, gains, bandE);
225    //quantise_pitch(gains, PBANDS);
226    pitch_quant_bands(st->mode, X, P, gains);
227    
228    //for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");
229    /* Subtract the pitch prediction from the signal to encode */
230    for (i=0;i<B*N;i++)
231       X[i] -= P[i];
232
233    /*float sum=0;
234    for (i=0;i<B*N;i++)
235       sum += X[i]*X[i];
236    printf ("%f\n", sum);*/
237    /* Residual quantisation */
238    quant_bands(st->mode, X, P, &st->enc);
239    
240    /* Synthesis */
241    denormalise_bands(st->mode, X, bandE);
242
243    //inv_haar1(X, B*N);
244
245    CELT_MOVE(st->out_mem, st->out_mem+B*N, MAX_PERIOD-B*N);
246    /* Compute inverse MDCTs */
247    for (i=0;i<B;i++)
248    {
249       int j;
250       float x[2*N];
251       float tmp[N];
252       /* De-interleaving the sub-frames */
253       for (j=0;j<N;j++)
254          tmp[j] = X[B*j+i];
255       mdct_backward(&st->mdct_lookup, tmp, x);
256       for (j=0;j<2*N;j++)
257          x[j] = st->window[j]*x[j];
258       for (j=0;j<N;j++)
259          st->out_mem[MAX_PERIOD+(i-B)*N+j] = x[j]+st->mdct_overlap[j];
260       for (j=0;j<N;j++)
261          st->mdct_overlap[j] = x[N+j];
262       
263       for (j=0;j<N;j++)
264       {
265          float tmp = st->out_mem[MAX_PERIOD+(i-B)*N+j] + st->preemph*st->preemph_memD;
266          st->preemph_memD = tmp;
267          pcm[i*N+j] = (short)floor(.5+tmp);
268       }
269    }
270    printf ("%d\n", ec_byte_bytes(&st->buf));
271    return 0;
272 }
273