xref: /aosp_15_r20/external/libopus/celt/laplace.c (revision a58d3d2adb790c104798cd88c8a3aff4fa8b82cc)
1*a58d3d2aSXin Li /* Copyright (c) 2007 CSIRO
2*a58d3d2aSXin Li    Copyright (c) 2007-2009 Xiph.Org Foundation
3*a58d3d2aSXin Li    Written by Jean-Marc Valin */
4*a58d3d2aSXin Li /*
5*a58d3d2aSXin Li    Redistribution and use in source and binary forms, with or without
6*a58d3d2aSXin Li    modification, are permitted provided that the following conditions
7*a58d3d2aSXin Li    are met:
8*a58d3d2aSXin Li 
9*a58d3d2aSXin Li    - Redistributions of source code must retain the above copyright
10*a58d3d2aSXin Li    notice, this list of conditions and the following disclaimer.
11*a58d3d2aSXin Li 
12*a58d3d2aSXin Li    - Redistributions in binary form must reproduce the above copyright
13*a58d3d2aSXin Li    notice, this list of conditions and the following disclaimer in the
14*a58d3d2aSXin Li    documentation and/or other materials provided with the distribution.
15*a58d3d2aSXin Li 
16*a58d3d2aSXin Li    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17*a58d3d2aSXin Li    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18*a58d3d2aSXin Li    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19*a58d3d2aSXin Li    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
20*a58d3d2aSXin Li    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21*a58d3d2aSXin Li    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22*a58d3d2aSXin Li    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23*a58d3d2aSXin Li    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24*a58d3d2aSXin Li    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25*a58d3d2aSXin Li    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26*a58d3d2aSXin Li    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27*a58d3d2aSXin Li */
28*a58d3d2aSXin Li 
29*a58d3d2aSXin Li #ifdef HAVE_CONFIG_H
30*a58d3d2aSXin Li #include "config.h"
31*a58d3d2aSXin Li #endif
32*a58d3d2aSXin Li 
33*a58d3d2aSXin Li #include "laplace.h"
34*a58d3d2aSXin Li #include "mathops.h"
35*a58d3d2aSXin Li 
36*a58d3d2aSXin Li /* The minimum probability of an energy delta (out of 32768). */
37*a58d3d2aSXin Li #define LAPLACE_LOG_MINP (0)
38*a58d3d2aSXin Li #define LAPLACE_MINP (1<<LAPLACE_LOG_MINP)
39*a58d3d2aSXin Li /* The minimum number of guaranteed representable energy deltas (in one
40*a58d3d2aSXin Li     direction). */
41*a58d3d2aSXin Li #define LAPLACE_NMIN (16)
42*a58d3d2aSXin Li 
43*a58d3d2aSXin Li /* When called, decay is positive and at most 11456. */
ec_laplace_get_freq1(unsigned fs0,int decay)44*a58d3d2aSXin Li static unsigned ec_laplace_get_freq1(unsigned fs0, int decay)
45*a58d3d2aSXin Li {
46*a58d3d2aSXin Li    unsigned ft;
47*a58d3d2aSXin Li    ft = 32768 - LAPLACE_MINP*(2*LAPLACE_NMIN) - fs0;
48*a58d3d2aSXin Li    return ft*(opus_int32)(16384-decay)>>15;
49*a58d3d2aSXin Li }
50*a58d3d2aSXin Li 
ec_laplace_encode(ec_enc * enc,int * value,unsigned fs,int decay)51*a58d3d2aSXin Li void ec_laplace_encode(ec_enc *enc, int *value, unsigned fs, int decay)
52*a58d3d2aSXin Li {
53*a58d3d2aSXin Li    unsigned fl;
54*a58d3d2aSXin Li    int val = *value;
55*a58d3d2aSXin Li    fl = 0;
56*a58d3d2aSXin Li    if (val)
57*a58d3d2aSXin Li    {
58*a58d3d2aSXin Li       int s;
59*a58d3d2aSXin Li       int i;
60*a58d3d2aSXin Li       s = -(val<0);
61*a58d3d2aSXin Li       val = (val+s)^s;
62*a58d3d2aSXin Li       fl = fs;
63*a58d3d2aSXin Li       fs = ec_laplace_get_freq1(fs, decay);
64*a58d3d2aSXin Li       /* Search the decaying part of the PDF.*/
65*a58d3d2aSXin Li       for (i=1; fs > 0 && i < val; i++)
66*a58d3d2aSXin Li       {
67*a58d3d2aSXin Li          fs *= 2;
68*a58d3d2aSXin Li          fl += fs+2*LAPLACE_MINP;
69*a58d3d2aSXin Li          fs = (fs*(opus_int32)decay)>>15;
70*a58d3d2aSXin Li       }
71*a58d3d2aSXin Li       /* Everything beyond that has probability LAPLACE_MINP. */
72*a58d3d2aSXin Li       if (!fs)
73*a58d3d2aSXin Li       {
74*a58d3d2aSXin Li          int di;
75*a58d3d2aSXin Li          int ndi_max;
76*a58d3d2aSXin Li          ndi_max = (32768-fl+LAPLACE_MINP-1)>>LAPLACE_LOG_MINP;
77*a58d3d2aSXin Li          ndi_max = (ndi_max-s)>>1;
78*a58d3d2aSXin Li          di = IMIN(val - i, ndi_max - 1);
79*a58d3d2aSXin Li          fl += (2*di+1+s)*LAPLACE_MINP;
80*a58d3d2aSXin Li          fs = IMIN(LAPLACE_MINP, 32768-fl);
81*a58d3d2aSXin Li          *value = (i+di+s)^s;
82*a58d3d2aSXin Li       }
83*a58d3d2aSXin Li       else
84*a58d3d2aSXin Li       {
85*a58d3d2aSXin Li          fs += LAPLACE_MINP;
86*a58d3d2aSXin Li          fl += fs&~s;
87*a58d3d2aSXin Li       }
88*a58d3d2aSXin Li       celt_assert(fl+fs<=32768);
89*a58d3d2aSXin Li       celt_assert(fs>0);
90*a58d3d2aSXin Li    }
91*a58d3d2aSXin Li    ec_encode_bin(enc, fl, fl+fs, 15);
92*a58d3d2aSXin Li }
93*a58d3d2aSXin Li 
ec_laplace_decode(ec_dec * dec,unsigned fs,int decay)94*a58d3d2aSXin Li int ec_laplace_decode(ec_dec *dec, unsigned fs, int decay)
95*a58d3d2aSXin Li {
96*a58d3d2aSXin Li    int val=0;
97*a58d3d2aSXin Li    unsigned fl;
98*a58d3d2aSXin Li    unsigned fm;
99*a58d3d2aSXin Li    fm = ec_decode_bin(dec, 15);
100*a58d3d2aSXin Li    fl = 0;
101*a58d3d2aSXin Li    if (fm >= fs)
102*a58d3d2aSXin Li    {
103*a58d3d2aSXin Li       val++;
104*a58d3d2aSXin Li       fl = fs;
105*a58d3d2aSXin Li       fs = ec_laplace_get_freq1(fs, decay)+LAPLACE_MINP;
106*a58d3d2aSXin Li       /* Search the decaying part of the PDF.*/
107*a58d3d2aSXin Li       while(fs > LAPLACE_MINP && fm >= fl+2*fs)
108*a58d3d2aSXin Li       {
109*a58d3d2aSXin Li          fs *= 2;
110*a58d3d2aSXin Li          fl += fs;
111*a58d3d2aSXin Li          fs = ((fs-2*LAPLACE_MINP)*(opus_int32)decay)>>15;
112*a58d3d2aSXin Li          fs += LAPLACE_MINP;
113*a58d3d2aSXin Li          val++;
114*a58d3d2aSXin Li       }
115*a58d3d2aSXin Li       /* Everything beyond that has probability LAPLACE_MINP. */
116*a58d3d2aSXin Li       if (fs <= LAPLACE_MINP)
117*a58d3d2aSXin Li       {
118*a58d3d2aSXin Li          int di;
119*a58d3d2aSXin Li          di = (fm-fl)>>(LAPLACE_LOG_MINP+1);
120*a58d3d2aSXin Li          val += di;
121*a58d3d2aSXin Li          fl += 2*di*LAPLACE_MINP;
122*a58d3d2aSXin Li       }
123*a58d3d2aSXin Li       if (fm < fl+fs)
124*a58d3d2aSXin Li          val = -val;
125*a58d3d2aSXin Li       else
126*a58d3d2aSXin Li          fl += fs;
127*a58d3d2aSXin Li    }
128*a58d3d2aSXin Li    celt_assert(fl<32768);
129*a58d3d2aSXin Li    celt_assert(fs>0);
130*a58d3d2aSXin Li    celt_assert(fl<=fm);
131*a58d3d2aSXin Li    celt_assert(fm<IMIN(fl+fs,32768));
132*a58d3d2aSXin Li    ec_dec_update(dec, fl, IMIN(fl+fs,32768), 32768);
133*a58d3d2aSXin Li    return val;
134*a58d3d2aSXin Li }
135*a58d3d2aSXin Li 
ec_laplace_encode_p0(ec_enc * enc,int value,opus_uint16 p0,opus_uint16 decay)136*a58d3d2aSXin Li void ec_laplace_encode_p0(ec_enc *enc, int value, opus_uint16 p0, opus_uint16 decay)
137*a58d3d2aSXin Li {
138*a58d3d2aSXin Li    int s;
139*a58d3d2aSXin Li    opus_uint16 sign_icdf[3];
140*a58d3d2aSXin Li    sign_icdf[0] = 32768-p0;
141*a58d3d2aSXin Li    sign_icdf[1] = sign_icdf[0]/2;
142*a58d3d2aSXin Li    sign_icdf[2] = 0;
143*a58d3d2aSXin Li    s = value == 0 ? 0 : (value > 0 ? 1 : 2);
144*a58d3d2aSXin Li    ec_enc_icdf16(enc, s, sign_icdf, 15);
145*a58d3d2aSXin Li    value = abs(value);
146*a58d3d2aSXin Li    if (value)
147*a58d3d2aSXin Li    {
148*a58d3d2aSXin Li       int i;
149*a58d3d2aSXin Li       opus_uint16 icdf[8];
150*a58d3d2aSXin Li       icdf[0] = IMAX(7, decay);
151*a58d3d2aSXin Li       for (i=1;i<7;i++)
152*a58d3d2aSXin Li       {
153*a58d3d2aSXin Li          icdf[i] = IMAX(7-i, (icdf[i-1] * (opus_int32)decay) >> 15);
154*a58d3d2aSXin Li       }
155*a58d3d2aSXin Li       icdf[7] = 0;
156*a58d3d2aSXin Li       value--;
157*a58d3d2aSXin Li       do {
158*a58d3d2aSXin Li          ec_enc_icdf16(enc, IMIN(value, 7), icdf, 15);
159*a58d3d2aSXin Li          value -= 7;
160*a58d3d2aSXin Li       } while (value >= 0);
161*a58d3d2aSXin Li    }
162*a58d3d2aSXin Li }
163*a58d3d2aSXin Li 
ec_laplace_decode_p0(ec_dec * dec,opus_uint16 p0,opus_uint16 decay)164*a58d3d2aSXin Li int ec_laplace_decode_p0(ec_dec *dec, opus_uint16 p0, opus_uint16 decay)
165*a58d3d2aSXin Li {
166*a58d3d2aSXin Li    int s;
167*a58d3d2aSXin Li    int value;
168*a58d3d2aSXin Li    opus_uint16 sign_icdf[3];
169*a58d3d2aSXin Li    sign_icdf[0] = 32768-p0;
170*a58d3d2aSXin Li    sign_icdf[1] = sign_icdf[0]/2;
171*a58d3d2aSXin Li    sign_icdf[2] = 0;
172*a58d3d2aSXin Li    s = ec_dec_icdf16(dec, sign_icdf, 15);
173*a58d3d2aSXin Li    if (s==2) s = -1;
174*a58d3d2aSXin Li    if (s != 0)
175*a58d3d2aSXin Li    {
176*a58d3d2aSXin Li       int i;
177*a58d3d2aSXin Li       int v;
178*a58d3d2aSXin Li       opus_uint16 icdf[8];
179*a58d3d2aSXin Li       icdf[0] = IMAX(7, decay);
180*a58d3d2aSXin Li       for (i=1;i<7;i++)
181*a58d3d2aSXin Li       {
182*a58d3d2aSXin Li          icdf[i] = IMAX(7-i, (icdf[i-1] * (opus_int32)decay) >> 15);
183*a58d3d2aSXin Li       }
184*a58d3d2aSXin Li       icdf[7] = 0;
185*a58d3d2aSXin Li       value = 1;
186*a58d3d2aSXin Li       do {
187*a58d3d2aSXin Li          v = ec_dec_icdf16(dec, icdf, 15);
188*a58d3d2aSXin Li          value += v;
189*a58d3d2aSXin Li       } while (v == 7);
190*a58d3d2aSXin Li       return s*value;
191*a58d3d2aSXin Li    } else return 0;
192*a58d3d2aSXin Li }
193*a58d3d2aSXin Li 
194*a58d3d2aSXin Li #if 0
195*a58d3d2aSXin Li 
196*a58d3d2aSXin Li #include <stdio.h>
197*a58d3d2aSXin Li #define NB_VALS 10
198*a58d3d2aSXin Li #define DATA_SIZE 10000
199*a58d3d2aSXin Li int main() {
200*a58d3d2aSXin Li    ec_enc enc;
201*a58d3d2aSXin Li    ec_dec dec;
202*a58d3d2aSXin Li    unsigned char *ptr;
203*a58d3d2aSXin Li    int i;
204*a58d3d2aSXin Li    int decay, p0;
205*a58d3d2aSXin Li    int val[NB_VALS] = {6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
206*a58d3d2aSXin Li    /*for (i=0;i<NB_VALS;i++) {
207*a58d3d2aSXin Li       val[i] = -log(rand()/(float)RAND_MAX);
208*a58d3d2aSXin Li       if (rand()%2) val[i] = -val[i];
209*a58d3d2aSXin Li    }*/
210*a58d3d2aSXin Li    p0 = 16000;
211*a58d3d2aSXin Li    decay = 16000;
212*a58d3d2aSXin Li    ptr = (unsigned char *)malloc(DATA_SIZE);
213*a58d3d2aSXin Li    ec_enc_init(&enc,ptr,DATA_SIZE);
214*a58d3d2aSXin Li    for (i=0;i<NB_VALS;i++) {
215*a58d3d2aSXin Li       printf("%d ", val[i]);
216*a58d3d2aSXin Li    }
217*a58d3d2aSXin Li    printf("\n");
218*a58d3d2aSXin Li    for (i=0;i<NB_VALS;i++) {
219*a58d3d2aSXin Li       ec_laplace_encode_p0(&enc, val[i], p0, decay);
220*a58d3d2aSXin Li    }
221*a58d3d2aSXin Li 
222*a58d3d2aSXin Li    ec_enc_done(&enc);
223*a58d3d2aSXin Li 
224*a58d3d2aSXin Li    ec_dec_init(&dec,ec_get_buffer(&enc),ec_range_bytes(&enc));
225*a58d3d2aSXin Li 
226*a58d3d2aSXin Li    for (i=0;i<NB_VALS;i++) {
227*a58d3d2aSXin Li       val[i] = ec_laplace_decode_p0(&dec, p0, decay);
228*a58d3d2aSXin Li    }
229*a58d3d2aSXin Li    for (i=0;i<NB_VALS;i++) {
230*a58d3d2aSXin Li       printf("%d ", val[i]);
231*a58d3d2aSXin Li    }
232*a58d3d2aSXin Li    printf("\n");
233*a58d3d2aSXin Li }
234*a58d3d2aSXin Li 
235*a58d3d2aSXin Li #endif
236