Doing init/reset of the entropy coder properly
[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
84    mdct_init(&st->mdct_lookup, 2*N);
85    st->fft = spx_fft_init(MAX_PERIOD);
86    
87    st->window = celt_alloc(2*N*sizeof(float));
88    st->in_mem = celt_alloc(N*sizeof(float));
89    st->mdct_overlap = celt_alloc(N*sizeof(float));
90    st->out_mem = celt_alloc(MAX_PERIOD*sizeof(float));
91    for (i=0;i<N;i++)
92       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));
93    
94    st->oldBandE = celt_alloc(mode->nbEBands*sizeof(float));
95
96    st->preemph = 0.8;
97    return st;
98 }
99
100 void celt_encoder_destroy(CELTState *st)
101 {
102    if (st == NULL)
103    {
104       celt_warning("NULL passed to celt_encoder_destroy");
105       return;
106    }
107    mdct_clear(&st->mdct_lookup);
108    spx_fft_destroy(st->fft);
109
110    celt_free(st->window);
111    celt_free(st->in_mem);
112    celt_free(st->mdct_overlap);
113    celt_free(st->out_mem);
114    
115    celt_free(st);
116 }
117
118 static void haar1(float *X, int N)
119 {
120    int i;
121    for (i=0;i<N;i+=2)
122    {
123       float a, b;
124       a = X[i];
125       b = X[i+1];
126       X[i] = .707107f*(a+b);
127       X[i+1] = .707107f*(a-b);
128    }
129 }
130
131 static void inv_haar1(float *X, int N)
132 {
133    int i;
134    for (i=0;i<N;i+=2)
135    {
136       float a, b;
137       a = X[i];
138       b = X[i+1];
139       X[i] = .707107f*(a+b);
140       X[i+1] = .707107f*(a-b);
141    }
142 }
143
144 static void compute_mdcts(mdct_lookup *mdct_lookup, float *window, float *in, float *out, int N, int B)
145 {
146    int i;
147    for (i=0;i<B;i++)
148    {
149       int j;
150       float x[2*N];
151       float tmp[N];
152       for (j=0;j<2*N;j++)
153          x[j] = window[j]*in[i*N+j];
154       mdct_forward(mdct_lookup, x, tmp);
155       /* Interleaving the sub-frames */
156       for (j=0;j<N;j++)
157          out[B*j+i] = tmp[j];
158    }
159
160 }
161
162 int celt_encode(CELTState *st, short *pcm)
163 {
164    int i, N, B;
165    N = st->block_size;
166    B = st->nb_blocks;
167    float in[(B+1)*N];
168    
169    float X[B*N];         /**< Interleaved signal MDCTs */
170    float P[B*N];         /**< Interleaved pitch MDCTs*/
171    float bandE[st->mode->nbEBands];
172    float gains[st->mode->nbPBands];
173    int pitch_index;
174    
175
176    ec_byte_reset(&st->buf);
177    ec_enc_init(&st->enc,&st->buf);
178
179    for (i=0;i<N;i++)
180       in[i] = st->in_mem[i];
181    for (;i<(B+1)*N;i++)
182    {
183       float tmp = pcm[i-N];
184       in[i] = tmp - st->preemph*st->preemph_memE;
185       st->preemph_memE = tmp;
186    }
187    for (i=0;i<N;i++)
188       st->in_mem[i] = in[B*N+i];
189
190    /* Compute MDCTs */
191    compute_mdcts(&st->mdct_lookup, st->window, in, X, N, B);
192    
193    /* Pitch analysis */
194    for (i=0;i<N;i++)
195    {
196       in[i] *= st->window[i];
197       in[B*N+i] *= st->window[N+i];
198    }
199    find_spectral_pitch(st->fft, in, st->out_mem, MAX_PERIOD, (B+1)*N, &pitch_index);
200    
201    /* Compute MDCTs of the pitch part */
202    compute_mdcts(&st->mdct_lookup, st->window, st->out_mem+pitch_index, P, N, B);
203    
204    /*int j;
205    for (j=0;j<B*N;j++)
206       printf ("%f ", X[j]);
207    for (j=0;j<B*N;j++)
208       printf ("%f ", P[j]);
209    printf ("\n");*/
210    //haar1(X, B*N);
211    //haar1(P, B*N);
212    
213    /* Band normalisation */
214    compute_band_energies(st->mode, X, bandE);
215    normalise_bands(st->mode, X, bandE);
216    //for (i=0;i<st->mode->nbEBands;i++)printf("%f ", bandE[i]);printf("\n");
217    
218    {
219       float bandEp[st->mode->nbEBands];
220       compute_band_energies(st->mode, P, bandEp);
221       normalise_bands(st->mode, P, bandEp);
222    }
223    
224    quant_energy(st->mode, bandE, st->oldBandE);
225    
226    /* Pitch prediction */
227    compute_pitch_gain(st->mode, X, P, gains, bandE);
228    //quantise_pitch(gains, PBANDS);
229    pitch_quant_bands(st->mode, X, P, gains);
230    
231    //for (i=0;i<B*N;i++) printf("%f ",P[i]);printf("\n");
232    /* Subtract the pitch prediction from the signal to encode */
233    for (i=0;i<B*N;i++)
234       X[i] -= P[i];
235
236    /*float sum=0;
237    for (i=0;i<B*N;i++)
238       sum += X[i]*X[i];
239    printf ("%f\n", sum);*/
240    /* Residual quantisation */
241    quant_bands(st->mode, X, P, &st->enc);
242    
243    /* Synthesis */
244    denormalise_bands(st->mode, X, bandE);
245
246    //inv_haar1(X, B*N);
247
248    CELT_MOVE(st->out_mem, st->out_mem+B*N, MAX_PERIOD-B*N);
249    /* Compute inverse MDCTs */
250    for (i=0;i<B;i++)
251    {
252       int j;
253       float x[2*N];
254       float tmp[N];
255       /* De-interleaving the sub-frames */
256       for (j=0;j<N;j++)
257          tmp[j] = X[B*j+i];
258       mdct_backward(&st->mdct_lookup, tmp, x);
259       for (j=0;j<2*N;j++)
260          x[j] = st->window[j]*x[j];
261       for (j=0;j<N;j++)
262          st->out_mem[MAX_PERIOD+(i-B)*N+j] = x[j]+st->mdct_overlap[j];
263       for (j=0;j<N;j++)
264          st->mdct_overlap[j] = x[N+j];
265       
266       for (j=0;j<N;j++)
267       {
268          float tmp = st->out_mem[MAX_PERIOD+(i-B)*N+j] + st->preemph*st->preemph_memD;
269          st->preemph_memD = tmp;
270          pcm[i*N+j] = (short)floor(.5+tmp);
271       }
272    }
273    ec_enc_done(&st->enc);
274    //printf ("%d\n", ec_byte_bytes(&st->buf));
275    return 0;
276 }
277