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