1 /************************************************************************
2 * Copyright (C) 2002-2009, Xiph.org Foundation
3 * Copyright (C) 2010, Robin Watts for Pinknoise Productions Ltd
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * * Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above
13 * copyright notice, this list of conditions and the following disclaimer
14 * in the documentation and/or other materials provided with the
15 * distribution.
16 * * Neither the names of the Xiph.org Foundation nor Pinknoise
17 * Productions Ltd nor the names of its contributors may be used to
18 * endorse or promote products derived from this software without
19 * specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 ************************************************************************
33
34 function: floor backend 0 implementation
35
36 ************************************************************************/
37
38 #include <stdlib.h>
39 #include <string.h>
40 #include <math.h>
41 #include <log/log.h>
42 #include "ogg.h"
43 #include "ivorbiscodec.h"
44 #include "codec_internal.h"
45 #include "codebook.h"
46 #include "misc.h"
47 #include "os.h"
48
49 #define LSP_FRACBITS 14
50 extern const ogg_int32_t FLOOR_fromdB_LOOKUP[];
51
52 /*************** LSP decode ********************/
53
54 #include "lsp_lookup.h"
55
56 /* interpolated 1./sqrt(p) where .5 <= a < 1. (.100000... to .111111...) in
57 16.16 format
58 returns in m.8 format */
59
60 static long ADJUST_SQRT2[2]={8192,5792};
vorbis_invsqlook_i(long a,long e)61 static inline ogg_int32_t vorbis_invsqlook_i(long a,long e){
62 long i=(a&0x7fff)>>(INVSQ_LOOKUP_I_SHIFT-1);
63 long d=a&INVSQ_LOOKUP_I_MASK; /* 0.10 */
64 long val=INVSQ_LOOKUP_I[i]- /* 1.16 */
65 ((INVSQ_LOOKUP_IDel[i]*d)>>INVSQ_LOOKUP_I_SHIFT); /* result 1.16 */
66 val*=ADJUST_SQRT2[e&1];
67 e=(e>>1)+21;
68 return(val>>e);
69 }
70
71 /* interpolated lookup based fromdB function, domain -140dB to 0dB only */
72 /* a is in n.12 format */
73 #ifdef _LOW_ACCURACY_
vorbis_fromdBlook_i(long a)74 static inline ogg_int32_t vorbis_fromdBlook_i(long a){
75 if(a>0) return 0x7fffffff;
76 if(a<(int)(((unsigned)-140)<<12)) return 0;
77 return FLOOR_fromdB_LOOKUP[((a+140)*467)>>20]<<9;
78 }
79 #else
vorbis_fromdBlook_i(long a)80 static inline ogg_int32_t vorbis_fromdBlook_i(long a){
81 if(a>0) return 0x7fffffff;
82 if(a<(int)(((unsigned)-140)<<12)) return 0;
83 return FLOOR_fromdB_LOOKUP[((a+(140<<12))*467)>>20];
84 }
85 #endif
86
87 /* interpolated lookup based cos function, domain 0 to PI only */
88 /* a is in 0.16 format, where 0==0, 2^^16-1==PI, return 0.14 */
vorbis_coslook_i(long a)89 static inline ogg_int32_t vorbis_coslook_i(long a){
90 int i=a>>COS_LOOKUP_I_SHIFT;
91 int d=a&COS_LOOKUP_I_MASK;
92 return COS_LOOKUP_I[i]- ((d*(COS_LOOKUP_I[i]-COS_LOOKUP_I[i+1]))>>
93 COS_LOOKUP_I_SHIFT);
94 }
95
96 /* interpolated half-wave lookup based cos function */
97 /* a is in 0.16 format, where 0==0, 2^^16==PI, return .LSP_FRACBITS */
vorbis_coslook2_i(long a)98 static inline ogg_int32_t vorbis_coslook2_i(long a){
99 int i=a>>COS_LOOKUP_I_SHIFT;
100 int d=a&COS_LOOKUP_I_MASK;
101 return ((COS_LOOKUP_I[i]<<COS_LOOKUP_I_SHIFT)-
102 d*(COS_LOOKUP_I[i]-COS_LOOKUP_I[i+1]))>>
103 (COS_LOOKUP_I_SHIFT-LSP_FRACBITS+14);
104 }
105
106 /* Values in barklook are defined such that toBARK(x) is an approximation to
107 POW(2, 15) * ((13.1*ATAN(0.00074*(x)))+(2.24*ATAN((x)*(x)*0.0000000185))+(0.0001*(x))) */
108 static const ogg_uint16_t barklook[]={
109 0,51,102,154, 206,258,311,365,
110 420,477,535,594, 656,719,785,854,
111 926,1002,1082,1166, 1256,1352,1454,1564,
112 1683,1812,1953,2107, 2276,2463,2670,2900,
113 3155,3440,3756,4106, 4493,4919,5387,5901,
114 6466,7094,7798,8599, 9528,10623,11935,13524,
115 15453,17775,20517,23667, 27183,31004,35069
116 };
117
118 /* used in init only; interpolate the long way */
toBARK(ogg_uint16_t n)119 static inline ogg_int32_t toBARK(ogg_uint16_t n){
120 int i;
121 int barklook_size = (sizeof(barklook) / sizeof(barklook[0]));
122 for(i=1;i<barklook_size;i++){
123 if(n<barklook[i]){
124 i--;
125 return (i<<14)+(((n-barklook[i])*
126 ((1UL<<31)/(barklook[i+1]-barklook[i])))>>17);
127 }
128 }
129 /* for a valid input n, which is half of info->rate (i.e. max 32767
130 as info->rate is 16 bit unsigned value), loop above will return
131 an output. So the following return will be used only when toBARK()
132 is called with invalid value */
133 return (barklook_size-1)<<14;
134 }
135
136 static const unsigned char MLOOP_1[64]={
137 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
138 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
139 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
140 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
141 };
142
143 static const unsigned char MLOOP_2[64]={
144 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
145 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
146 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
147 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
148 };
149
150 static const unsigned char MLOOP_3[8]={0,1,2,2,3,3,3,3};
151
vorbis_lsp_to_curve(ogg_int32_t * curve,int n,int ln,ogg_int32_t * lsp,int m,ogg_int32_t amp,ogg_int32_t ampoffset,ogg_int32_t nyq)152 void vorbis_lsp_to_curve(ogg_int32_t *curve,int n,int ln,
153 ogg_int32_t *lsp,int m,
154 ogg_int32_t amp,
155 ogg_int32_t ampoffset,
156 ogg_int32_t nyq){
157
158 /* 0 <= m < 256 */
159
160 /* set up for using all int later */
161 int i;
162 int ampoffseti=ampoffset*4096;
163 int ampi=amp;
164 ogg_int32_t *ilsp=(ogg_int32_t *)alloca(m*sizeof(*ilsp));
165
166 ogg_uint32_t inyq= (1UL<<31) / toBARK(nyq);
167 ogg_uint32_t imap= (1UL<<31) / ln;
168 ogg_uint32_t tBnyq1 = toBARK(nyq)<<1;
169
170 /* Besenham for frequency scale to avoid a division */
171 int f=0;
172 int fdx=n;
173 int fbase=nyq/fdx;
174 int ferr=0;
175 int fdy=nyq-fbase*fdx;
176 int map=0;
177
178 #ifdef _LOW_ACCURACY_
179 ogg_uint32_t nextbark=((tBnyq1<<11)/ln)>>12;
180 #else
181 ogg_uint32_t nextbark=MULT31(imap>>1,tBnyq1);
182 #endif
183 /* nextbark is guaranteed to be less than 54 << 14 here and that ensures index
184 to barklook can at max be 53 and 54 here */
185 int nextf=barklook[nextbark>>14]+(((nextbark&0x3fff)*
186 (barklook[(nextbark>>14)+1]-barklook[nextbark>>14]))>>14);
187
188 /* lsp is in 8.24, range 0 to PI; coslook wants it in .16 0 to 1*/
189 for(i=0;i<m;i++){
190 #ifndef _LOW_ACCURACY_
191 ogg_int32_t val=MULT32(lsp[i],0x517cc2);
192 #else
193 ogg_int32_t val=((lsp[i]>>10)*0x517d)>>14;
194 #endif
195
196 /* safeguard against a malicious stream */
197 if(val<0 || (val>>COS_LOOKUP_I_SHIFT)>=COS_LOOKUP_I_SZ){
198 memset(curve,0,sizeof(*curve)*n);
199 return;
200 }
201
202 ilsp[i]=vorbis_coslook_i(val);
203 }
204
205 i=0;
206 while(i<n){
207 int j;
208 ogg_uint32_t pi=46341; /* 2**-.5 in 0.16 */
209 ogg_uint32_t qi=46341;
210 ogg_int32_t qexp=0,shift;
211 ogg_int32_t wi;
212
213 wi=vorbis_coslook2_i((map*imap)>>15);
214
215
216 #ifdef _V_LSP_MATH_ASM
217 lsp_loop_asm(&qi,&pi,&qexp,ilsp,wi,m);
218
219 pi=((pi*pi)>>16);
220 qi=((qi*qi)>>16);
221
222 if(m&1){
223 qexp= qexp*2-28*((m+1)>>1)+m;
224 pi*=(1<<14)-((wi*wi)>>14);
225 qi+=pi>>14;
226 }else{
227 qexp= qexp*2-13*m;
228
229 pi*=(1<<14)-wi;
230 qi*=(1<<14)+wi;
231
232 qi=(qi+pi)>>14;
233 }
234
235 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
236 qi>>=1; qexp++;
237 }else
238 lsp_norm_asm(&qi,&qexp);
239
240 #else
241
242 qi*=labs(ilsp[0]-wi);
243 pi*=labs(ilsp[1]-wi);
244
245 for(j=3;j<m;j+=2){
246 if(!(shift=MLOOP_1[(pi|qi)>>25]))
247 if(!(shift=MLOOP_2[(pi|qi)>>19]))
248 shift=MLOOP_3[(pi|qi)>>16];
249
250 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
251 pi=(pi>>shift)*labs(ilsp[j]-wi);
252 qexp+=shift;
253 }
254 if(!(shift=MLOOP_1[(pi|qi)>>25]))
255 if(!(shift=MLOOP_2[(pi|qi)>>19]))
256 shift=MLOOP_3[(pi|qi)>>16];
257
258 /* pi,qi normalized collectively, both tracked using qexp */
259
260 if(m&1){
261 /* odd order filter; slightly assymetric */
262 /* the last coefficient */
263 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
264 pi=(pi>>shift)<<14;
265 qexp+=shift;
266
267 if(!(shift=MLOOP_1[(pi|qi)>>25]))
268 if(!(shift=MLOOP_2[(pi|qi)>>19]))
269 shift=MLOOP_3[(pi|qi)>>16];
270
271 pi>>=shift;
272 qi>>=shift;
273 qexp+=shift-14*((m+1)>>1);
274
275 pi=((pi*pi)>>16);
276 qi=((qi*qi)>>16);
277 qexp=qexp*2+m;
278
279 pi*=(1<<14)-((wi*wi)>>14);
280 qi+=pi>>14;
281
282 }else{
283 /* even order filter; still symmetric */
284
285 /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
286 worth tracking step by step */
287
288 pi>>=shift;
289 qi>>=shift;
290 qexp+=shift-7*m;
291
292 pi=((pi*pi)>>16);
293 qi=((qi*qi)>>16);
294 qexp=qexp*2+m;
295
296 pi*=(1<<14)-wi;
297 qi*=(1<<14)+wi;
298 qi=(qi+pi)>>14;
299
300 }
301
302
303 /* we've let the normalization drift because it wasn't important;
304 however, for the lookup, things must be normalized again. We
305 need at most one right shift or a number of left shifts */
306
307 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
308 qi>>=1; qexp++;
309 }else
310 while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
311 qi<<=1; qexp--;
312 }
313
314 #endif
315
316 amp=vorbis_fromdBlook_i(ampi* /* n.4 */
317 vorbis_invsqlook_i(qi,qexp)-
318 /* m.8, m+n<=8 */
319 ampoffseti); /* 8.12[0] */
320
321 #ifdef _LOW_ACCURACY_
322 amp>>=9;
323 #endif
324 curve[i]= MULT31_SHIFT15(curve[i],amp);
325
326 while(++i<n){
327
328 /* line plot to get new f */
329 ferr+=fdy;
330 if(ferr>=fdx){
331 ferr-=fdx;
332 f++;
333 }
334 f+=fbase;
335
336 if(f>=nextf)break;
337
338 curve[i]= MULT31_SHIFT15(curve[i],amp);
339 }
340
341 while(1){
342 map++;
343
344 if(map+1<ln){
345
346 #ifdef _LOW_ACCURACY_
347 nextbark=((tBnyq1<<11)/ln*(map+1))>>12;
348 #else
349 nextbark=MULT31((map+1)*(imap>>1),tBnyq1);
350 #endif
351 /* nextbark is guaranteed to be less than 54 << 14 here and that ensures index
352 to barklook can at max be 53 and 54 here */
353 nextf=barklook[nextbark>>14]+
354 (((nextbark&0x3fff)*
355 (barklook[(nextbark>>14)+1]-barklook[nextbark>>14]))>>14);
356 if(f<=nextf)break;
357
358 }else{
359 nextf=9999999;
360 break;
361 }
362 }
363 if(map>=ln){
364 map=ln-1; /* guard against the approximation */
365 nextf=9999999;
366 }
367 }
368 }
369
370 /*************** vorbis decode glue ************/
371
floor0_free_info(vorbis_info_floor * i)372 void floor0_free_info(vorbis_info_floor *i){
373 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
374 if(info)_ogg_free(info);
375 }
376
floor0_info_unpack(vorbis_info * vi,oggpack_buffer * opb)377 vorbis_info_floor *floor0_info_unpack (vorbis_info *vi,oggpack_buffer *opb){
378 codec_setup_info *ci=(codec_setup_info *)vi->codec_setup;
379 int j;
380
381 vorbis_info_floor0 *info=(vorbis_info_floor0 *)_ogg_malloc(sizeof(*info));
382 info->order=oggpack_read(opb,8);
383 info->rate=oggpack_read(opb,16);
384 info->barkmap=oggpack_read(opb,16);
385 info->ampbits=oggpack_read(opb,6);
386 info->ampdB=oggpack_read(opb,8);
387 info->numbooks=oggpack_read(opb,4)+1;
388
389 /* order must be greater than 1 to calculate p and q parameters for the linear floor value */
390 if(info->order<=1)goto err_out;
391 if(info->rate<1)goto err_out;
392 if(info->barkmap<1)goto err_out;
393
394 for(j=0;j<info->numbooks;j++){
395 info->books[j]=(char)oggpack_read(opb,8);
396 if(info->books[j]>=ci->books)goto err_out;
397 }
398
399 if(oggpack_eop(opb))goto err_out;
400 return(info);
401
402 err_out:
403 floor0_free_info(info);
404 return(NULL);
405 }
406
floor0_memosize(vorbis_info_floor * i)407 int floor0_memosize(vorbis_info_floor *i){
408 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
409 return info->order+1;
410 }
411
floor0_inverse1(vorbis_dsp_state * vd,vorbis_info_floor * i,ogg_int32_t * lsp)412 ogg_int32_t *floor0_inverse1(vorbis_dsp_state *vd,vorbis_info_floor *i,
413 ogg_int32_t *lsp){
414 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
415 int j,k;
416 uint64_t ampraw;
417 if(info->ampbits<=32){
418 ampraw=oggpack_read(&vd->opb,info->ampbits);
419 }else{
420 //max value possible for info->ampbits is 63
421 ampraw=oggpack_read(&vd->opb,info->ampbits-32);
422 ampraw<<=32;
423 ampraw|=oggpack_read(&vd->opb,32);
424 }
425 if(ampraw>0){ /* also handles the -1 out of data case */
426 uint64_t maxval=(1<<info->ampbits)-1;
427 float ampRatio=(float)ampraw/maxval;
428 int amp=ampRatio*(info->ampdB<<4);
429 int booknum=oggpack_read(&vd->opb,_ilog(info->numbooks));
430 if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
431 codec_setup_info *ci=(codec_setup_info *)vd->vi->codec_setup;
432 codebook *b=ci->book_param+info->books[booknum];
433 int sz = floor0_memosize(i);
434 if (sz < b->dim) {
435 ALOGE("lsp too small: %d < %ld", sz, b->dim);
436 return NULL;
437 }
438 ogg_int32_t last=0;
439
440 if(vorbis_book_decodev_set(b,lsp,&vd->opb,info->order,-24)==-1)goto eop;
441 for(j=0;j<info->order;){
442 for(k=0;k<b->dim && j<info->order;k++,j++)lsp[j]+=last;
443 last=lsp[j-1];
444 }
445
446 lsp[info->order]=amp;
447 return(lsp);
448 }
449 }
450 eop:
451 return(NULL);
452 }
453
floor0_inverse2(vorbis_dsp_state * vd,vorbis_info_floor * i,ogg_int32_t * lsp,ogg_int32_t * out)454 int floor0_inverse2(vorbis_dsp_state *vd,vorbis_info_floor *i,
455 ogg_int32_t *lsp,ogg_int32_t *out){
456 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
457 codec_setup_info *ci=(codec_setup_info *)vd->vi->codec_setup;
458
459 if(lsp){
460 ogg_int32_t amp=lsp[info->order];
461
462 /* take the coefficients back to a spectral envelope curve */
463 vorbis_lsp_to_curve(out,ci->blocksizes[vd->W]/2,info->barkmap,
464 lsp,info->order,amp,info->ampdB,
465 info->rate>>1);
466 return(1);
467 }
468 memset(out,0,sizeof(*out)*ci->blocksizes[vd->W]/2);
469 return(0);
470 }
471
472