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