1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
16    - Neither the name of the Xiph.org Foundation nor the names of its
17    contributors may be used to endorse or promote products derived from
18    this software without specific prior written permission.
20    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
33 #ifdef HAVE_CONFIG_H
34 #include "config.h"
35 #endif
37 #include "mathops.h"
38 #include "cwrs.h"
39 #include "vq.h"
40 #include "arch.h"
41 #include "os_support.h"
42 #include "rate.h"
44 #ifndef M_PI
45 #define M_PI 3.141592653
46 #endif
48 static void exp_rotation(celt_norm *X, int len, int dir, int stride, int K)
49 {
50    int i, k, iter;
51    celt_word16 c, s;
52    celt_word16 gain, theta;
53    celt_norm *Xptr;
54    gain = celt_div((celt_word32)MULT16_16(Q15_ONE,len),(celt_word32)(3+len+6*K));
55    /* FIXME: Make that HALF16 instead of HALF32 */
56    theta = SUB16(Q15ONE, HALF32(MULT16_16_Q15(gain,gain)));
57    /*if (len==30)
58    {
59    for (i=0;i<len;i++)
60    X[i] = 0;
61    X = 1;
62 }*/
63    c = celt_cos_norm(EXTEND32(theta));
64    s = dir*celt_cos_norm(EXTEND32(SUB16(Q15ONE,theta))); /*  sin(theta) */
65    if (len > 8*stride)
66       stride *= len/(8*stride);
67    iter = 1;
68    for (k=0;k<iter;k++)
69    {
70       /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy,
71       but at this point, I really don't think it's necessary */
72       Xptr = X;
73       for (i=0;i<len-stride;i++)
74       {
75          celt_norm x1, x2;
76          x1 = Xptr;
77          x2 = Xptr[stride];
78          Xptr[stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
79          *Xptr++      = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
80       }
81       Xptr = &X[len-2*stride-1];
82       for (i=len-2*stride-1;i>=0;i--)
83       {
84          celt_norm x1, x2;
85          x1 = Xptr;
86          x2 = Xptr[stride];
87          Xptr[stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
88          *Xptr--      = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
89       }
90    }
91    /*if (len==30)
92    {
93    for (i=0;i<len;i++)
94    printf ("%f ", X[i]);
95    printf ("\n");
96    exit(0);
97 }*/
98 }
101 /** Takes the pitch vector and the decoded residual vector, computes the gain
102     that will give ||p+g*y||=1 and mixes the residual with the pitch. */
103 static void normalise_residual(int * restrict iy, celt_norm * restrict X, int N, int K, celt_word32 Ryy)
104 {
105    int i;
106    celt_word32 g;
108    g = celt_rsqrt(Ryy);
110    i=0;
111    do
112       X[i] = EXTRACT16(SHR32(MULT16_16(g, iy[i]),1));
113    while (++i < N);
114 }
116 void alg_quant(celt_norm *X, int N, int K, int spread, ec_enc *enc)
117 {
118    VARDECL(celt_norm, y);
119    VARDECL(int, iy);
120    VARDECL(celt_word16, signx);
121    int j, is;
122    celt_word16 s;
123    int pulsesLeft;
124    celt_word32 sum;
125    celt_word32 xy, yy;
126    int N_1; /* Inverse of N, in Q14 format (even for float) */
127 #ifdef FIXED_POINT
128    int yshift;
129 #endif
130    SAVE_STACK;
132    K = get_pulses(K);
133 #ifdef FIXED_POINT
134    yshift = 13-celt_ilog2(K);
135 #endif
137    ALLOC(y, N, celt_norm);
138    ALLOC(iy, N, int);
139    ALLOC(signx, N, celt_word16);
140    N_1 = 512/N;
143       exp_rotation(X, N, 1, spread, K);
145    sum = 0;
146    j=0; do {
147       if (X[j]>0)
148          signx[j]=1;
149       else {
150          signx[j]=-1;
151          X[j]=-X[j];
152       }
153       iy[j] = 0;
154       y[j] = 0;
155    } while (++j<N);
157    xy = yy = 0;
159    pulsesLeft = K;
161    /* Do a pre-search by projecting on the pyramid */
162    if (K > (N>>1))
163    {
164       celt_word16 rcp;
165       j=0; do {
166          sum += X[j];
167       }  while (++j<N);
169 #ifdef FIXED_POINT
170       if (sum <= K)
171 #else
172       if (sum <= EPSILON)
173 #endif
174       {
175          X = QCONST16(1.f,14);
176          j=1; do
177             X[j]=0;
178          while (++j<N);
179          sum = QCONST16(1.f,14);
180       }
181       /* Do we have sufficient accuracy here? */
182       rcp = EXTRACT16(MULT16_32_Q16(K-1, celt_rcp(sum)));
183       j=0; do {
184 #ifdef FIXED_POINT
185          /* It's really important to round *towards zero* here */
186          iy[j] = MULT16_16_Q15(X[j],rcp);
187 #else
188          iy[j] = floor(rcp*X[j]);
189 #endif
190          y[j] = SHL16(iy[j],yshift);
191          yy = MAC16_16(yy, y[j],y[j]);
192          xy = MAC16_16(xy, X[j],y[j]);
193          y[j] *= 2;
194          pulsesLeft -= iy[j];
195       }  while (++j<N);
196    }
197    celt_assert2(pulsesLeft>=1, "Allocated too many pulses in the quick pass");
199    while (pulsesLeft > 0)
200    {
201       int pulsesAtOnce=1;
202       int best_id;
203       celt_word16 magnitude;
204       celt_word32 best_num = -VERY_LARGE16;
205       celt_word16 best_den = 0;
206 #ifdef FIXED_POINT
207       int rshift;
208 #endif
209       /* Decide on how many pulses to find at once */
210       pulsesAtOnce = (pulsesLeft*N_1)>>9; /* pulsesLeft/N */
211       if (pulsesAtOnce<1)
212          pulsesAtOnce = 1;
213 #ifdef FIXED_POINT
214       rshift = yshift+1+celt_ilog2(K-pulsesLeft+pulsesAtOnce);
215 #endif
216       magnitude = SHL16(pulsesAtOnce, yshift);
218       best_id = 0;
219       /* The squared magnitude term gets added anyway, so we might as well
220          add it outside the loop */
221       yy = MAC16_16(yy, magnitude,magnitude);
222       /* Choose between fast and accurate strategy depending on where we are in the search */
223          /* This should ensure that anything we can process will have a better score */
224       j=0;
225       do {
226          celt_word16 Rxy, Ryy;
227          /* Select sign based on X[j] alone */
228          s = magnitude;
229          /* Temporary sums of the new pulse(s) */
230          Rxy = EXTRACT16(SHR32(MAC16_16(xy, s,X[j]),rshift));
231          /* We're multiplying y[j] by two so we don't have to do it here */
232          Ryy = EXTRACT16(SHR32(MAC16_16(yy, s,y[j]),rshift));
234             /* Approximate score: we maximise Rxy/sqrt(Ryy) (we're guaranteed that
235          Rxy is positive because the sign is pre-computed) */
236          Rxy = MULT16_16_Q15(Rxy,Rxy);
237             /* The idea is to check for num/den >= best_num/best_den, but that way
238          we can do it without any division */
239          /* OPT: Make sure to use conditional moves here */
240          if (MULT16_16(best_den, Rxy) > MULT16_16(Ryy, best_num))
241          {
242             best_den = Ryy;
243             best_num = Rxy;
244             best_id = j;
245          }
246       } while (++j<N);
248       j = best_id;
249       is = pulsesAtOnce;
250       s = SHL16(is, yshift);
252       /* Updating the sums of the new pulse(s) */
253       xy = xy + MULT16_16(s,X[j]);
254       /* We're multiplying y[j] by two so we don't have to do it here */
255       yy = yy + MULT16_16(s,y[j]);
257       /* Only now that we've made the final choice, update y/iy */
258       /* Multiplying y[j] by 2 so we don't have to do it everywhere else */
259       y[j] += 2*s;
260       iy[j] += is;
261       pulsesLeft -= pulsesAtOnce;
262    }
263    j=0;
264    do {
265       X[j] = MULT16_16(signx[j],X[j]);
266       if (signx[j] < 0)
267          iy[j] = -iy[j];
268    } while (++j<N);
269    encode_pulses(iy, N, K, enc);
271    /* Recompute the gain in one pass to reduce the encoder-decoder mismatch
272    due to the recursive computation used in quantisation. */
273    normalise_residual(iy, X, N, K, EXTRACT16(SHR32(yy,2*yshift)));
275       exp_rotation(X, N, -1, spread, K);
276    RESTORE_STACK;
277 }
280 /** Decode pulse vector and combine the result with the pitch vector to produce
281     the final normalised signal in the current band. */
282 void alg_unquant(celt_norm *X, int N, int K, int spread, ec_dec *dec)
283 {
284    int i;
285    celt_word32 Ryy;
286    VARDECL(int, iy);
287    SAVE_STACK;
288    K = get_pulses(K);
289    ALLOC(iy, N, int);
290    decode_pulses(iy, N, K, dec);
291    Ryy = 0;
292    i=0;
293    do {
294       Ryy = MAC16_16(Ryy, iy[i], iy[i]);
295    } while (++i < N);
296    normalise_residual(iy, X, N, K, Ryy);
298       exp_rotation(X, N, -1, spread, K);
299    RESTORE_STACK;
300 }
302 celt_word16 renormalise_vector(celt_norm *X, celt_word16 value, int N, int stride)
303 {
304    int i;
305    celt_word32 E = EPSILON;
306    celt_word16 rE;
307    celt_word16 g;
308    celt_norm *xptr = X;
309    for (i=0;i<N;i++)
310    {
311       E = MAC16_16(E, *xptr, *xptr);
312       xptr += stride;
313    }
315    rE = celt_sqrt(E);
316 #ifdef FIXED_POINT
317    if (rE <= 128)
318       g = Q15ONE;
319    else
320 #endif
321       g = MULT16_16_Q15(value,celt_rcp(SHL32(rE,9)));
322    xptr = X;
323    for (i=0;i<N;i++)
324    {
325       *xptr = PSHR32(MULT16_16(g, *xptr),8);
326       xptr += stride;
327    }
328    return rE;
329 }
331 static void fold(const CELTMode *m, int N, const celt_norm * restrict Y, celt_norm * restrict P, int N0, int B)
332 {
333    int j;
334    int id = N0 % B;
335    /* Here, we assume that id will never be greater than N0, i.e. that
336       no band is wider than N0. In the unlikely case it happens, we set
337       everything to zero */
338    /*{
339            int offset = (N0*C - (id+C*N))/2;
340            if (offset > C*N0/16)
341                    offset = C*N0/16;
342            offset -= offset % (C*B);
343            if (offset < 0)
344                    offset = 0;
345            //printf ("%d\n", offset);
346            id += offset;
347    }*/
348    if (id+N>N0)
349       for (j=0;j<N;j++)
350          P[j] = 0;
351    else
352       for (j=0;j<N;j++)
353          P[j] = Y[id++];
354 }
356 void intra_fold(const CELTMode *m, int N, const celt_norm * restrict Y, celt_norm * restrict P, int N0, int B)
357 {
358    fold(m, N, Y, P, N0, B);
359    renormalise_vector(P, Q15ONE, N, 1);
360 }
{"ft_lang_label":"__label__en","ft_lang_prob":0.6266636,"math_prob":0.954808,"size":6982,"snap":"2019-51-2020-05","text_gpt3_token_len":2093,"char_repetition_ratio":0.11149327,"word_repetition_ratio":0.067780875,"special_character_ratio":0.37711257,"punctuation_ratio":0.16938342,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9802923,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2019-12-10T05:53:42Z\",\"WARC-Record-ID\":\"<urn:uuid:ab17c59b-c359-41c3-bbbc-c320218a9482>\",\"Content-Length\":\"82113\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:659cf175-2e88-4c3f-ab47-7a1ab1f5dfa4>\",\"WARC-Concurrent-To\":\"<urn:uuid:59b595a3-d8c6-4ec6-8e16-afbfd4b5bd5a>\",\"WARC-IP-Address\":\"140.211.15.28\",\"WARC-Target-URI\":\"https://git.xiph.org/?p=opus.git;a=blob;f=libcelt/vq.c;hb=a3803a98afeb8bddfee2038ecdf0623124d06d46\",\"WARC-Payload-Digest\":\"sha1:6WP25VGWMHWMLP56I4SS3W25FHEBUW6D\",\"WARC-Block-Digest\":\"sha1:WNF4LURRBIF6BUNVJ6KBGFN7XFCM7MVW\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2019/CC-MAIN-2019-51/CC-MAIN-2019-51_segments_1575540525821.56_warc_CC-MAIN-20191210041836-20191210065836-00307.warc.gz\"}"}