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