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 static const ogg_uint16_t barklook[54]={
107 0,51,102,154, 206,258,311,365,
108 420,477,535,594, 656,719,785,854,
109 926,1002,1082,1166, 1256,1352,1454,1564,
110 1683,1812,1953,2107, 2276,2463,2670,2900,
111 3155,3440,3756,4106, 4493,4919,5387,5901,
112 6466,7094,7798,8599, 9528,10623,11935,13524,
113 15453,17775,20517,23667, 27183,31004
114 };
115
116 /* used in init only; interpolate the long way */
toBARK(int n)117 static inline ogg_int32_t toBARK(int n){
118 int i;
119 for(i=0;i<54;i++)
120 if(n>=barklook[i] && n<barklook[i+1])break;
121
122 if(i==54){
123 return 54<<14;
124 }else{
125 return (i<<14)+(((n-barklook[i])*
126 ((1UL<<31)/(barklook[i+1]-barklook[i])))>>17);
127 }
128 }
129
130 static const unsigned char MLOOP_1[64]={
131 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
132 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
133 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
134 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
135 };
136
137 static const unsigned char MLOOP_2[64]={
138 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
139 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
140 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
141 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
142 };
143
144 static const unsigned char MLOOP_3[8]={0,1,2,2,3,3,3,3};
145
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)146 void vorbis_lsp_to_curve(ogg_int32_t *curve,int n,int ln,
147 ogg_int32_t *lsp,int m,
148 ogg_int32_t amp,
149 ogg_int32_t ampoffset,
150 ogg_int32_t nyq){
151
152 /* 0 <= m < 256 */
153
154 /* set up for using all int later */
155 int i;
156 int ampoffseti=ampoffset*4096;
157 int ampi=amp;
158 ogg_int32_t *ilsp=(ogg_int32_t *)alloca(m*sizeof(*ilsp));
159
160 ogg_uint32_t inyq= (1UL<<31) / toBARK(nyq);
161 ogg_uint32_t imap= (1UL<<31) / ln;
162 ogg_uint32_t tBnyq1 = toBARK(nyq)<<1;
163
164 /* Besenham for frequency scale to avoid a division */
165 int f=0;
166 int fdx=n;
167 int fbase=nyq/fdx;
168 int ferr=0;
169 int fdy=nyq-fbase*fdx;
170 int map=0;
171
172 #ifdef _LOW_ACCURACY_
173 ogg_uint32_t nextbark=((tBnyq1<<11)/ln)>>12;
174 #else
175 ogg_uint32_t nextbark=MULT31(imap>>1,tBnyq1);
176 #endif
177 int nextf=barklook[nextbark>>14]+(((nextbark&0x3fff)*
178 (barklook[(nextbark>>14)+1]-barklook[nextbark>>14]))>>14);
179
180 /* lsp is in 8.24, range 0 to PI; coslook wants it in .16 0 to 1*/
181 for(i=0;i<m;i++){
182 #ifndef _LOW_ACCURACY_
183 ogg_int32_t val=MULT32(lsp[i],0x517cc2);
184 #else
185 ogg_int32_t val=((lsp[i]>>10)*0x517d)>>14;
186 #endif
187
188 /* safeguard against a malicious stream */
189 if(val<0 || (val>>COS_LOOKUP_I_SHIFT)>=COS_LOOKUP_I_SZ){
190 memset(curve,0,sizeof(*curve)*n);
191 return;
192 }
193
194 ilsp[i]=vorbis_coslook_i(val);
195 }
196
197 i=0;
198 while(i<n){
199 int j;
200 ogg_uint32_t pi=46341; /* 2**-.5 in 0.16 */
201 ogg_uint32_t qi=46341;
202 ogg_int32_t qexp=0,shift;
203 ogg_int32_t wi;
204
205 wi=vorbis_coslook2_i((map*imap)>>15);
206
207
208 #ifdef _V_LSP_MATH_ASM
209 lsp_loop_asm(&qi,&pi,&qexp,ilsp,wi,m);
210
211 pi=((pi*pi)>>16);
212 qi=((qi*qi)>>16);
213
214 if(m&1){
215 qexp= qexp*2-28*((m+1)>>1)+m;
216 pi*=(1<<14)-((wi*wi)>>14);
217 qi+=pi>>14;
218 }else{
219 qexp= qexp*2-13*m;
220
221 pi*=(1<<14)-wi;
222 qi*=(1<<14)+wi;
223
224 qi=(qi+pi)>>14;
225 }
226
227 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
228 qi>>=1; qexp++;
229 }else
230 lsp_norm_asm(&qi,&qexp);
231
232 #else
233
234 qi*=labs(ilsp[0]-wi);
235 pi*=labs(ilsp[1]-wi);
236
237 for(j=3;j<m;j+=2){
238 if(!(shift=MLOOP_1[(pi|qi)>>25]))
239 if(!(shift=MLOOP_2[(pi|qi)>>19]))
240 shift=MLOOP_3[(pi|qi)>>16];
241
242 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
243 pi=(pi>>shift)*labs(ilsp[j]-wi);
244 qexp+=shift;
245 }
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 /* pi,qi normalized collectively, both tracked using qexp */
251
252 if(m&1){
253 /* odd order filter; slightly assymetric */
254 /* the last coefficient */
255 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
256 pi=(pi>>shift)<<14;
257 qexp+=shift;
258
259 if(!(shift=MLOOP_1[(pi|qi)>>25]))
260 if(!(shift=MLOOP_2[(pi|qi)>>19]))
261 shift=MLOOP_3[(pi|qi)>>16];
262
263 pi>>=shift;
264 qi>>=shift;
265 qexp+=shift-14*((m+1)>>1);
266
267 pi=((pi*pi)>>16);
268 qi=((qi*qi)>>16);
269 qexp=qexp*2+m;
270
271 pi*=(1<<14)-((wi*wi)>>14);
272 qi+=pi>>14;
273
274 }else{
275 /* even order filter; still symmetric */
276
277 /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
278 worth tracking step by step */
279
280 pi>>=shift;
281 qi>>=shift;
282 qexp+=shift-7*m;
283
284 pi=((pi*pi)>>16);
285 qi=((qi*qi)>>16);
286 qexp=qexp*2+m;
287
288 pi*=(1<<14)-wi;
289 qi*=(1<<14)+wi;
290 qi=(qi+pi)>>14;
291
292 }
293
294
295 /* we've let the normalization drift because it wasn't important;
296 however, for the lookup, things must be normalized again. We
297 need at most one right shift or a number of left shifts */
298
299 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
300 qi>>=1; qexp++;
301 }else
302 while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
303 qi<<=1; qexp--;
304 }
305
306 #endif
307
308 amp=vorbis_fromdBlook_i(ampi* /* n.4 */
309 vorbis_invsqlook_i(qi,qexp)-
310 /* m.8, m+n<=8 */
311 ampoffseti); /* 8.12[0] */
312
313 #ifdef _LOW_ACCURACY_
314 amp>>=9;
315 #endif
316 curve[i]= MULT31_SHIFT15(curve[i],amp);
317
318 while(++i<n){
319
320 /* line plot to get new f */
321 ferr+=fdy;
322 if(ferr>=fdx){
323 ferr-=fdx;
324 f++;
325 }
326 f+=fbase;
327
328 if(f>=nextf)break;
329
330 curve[i]= MULT31_SHIFT15(curve[i],amp);
331 }
332
333 while(1){
334 map++;
335
336 if(map+1<ln){
337
338 #ifdef _LOW_ACCURACY_
339 nextbark=((tBnyq1<<11)/ln*(map+1))>>12;
340 #else
341 nextbark=MULT31((map+1)*(imap>>1),tBnyq1);
342 #endif
343 nextf=barklook[nextbark>>14]+
344 (((nextbark&0x3fff)*
345 (barklook[(nextbark>>14)+1]-barklook[nextbark>>14]))>>14);
346 if(f<=nextf)break;
347
348 }else{
349 nextf=9999999;
350 break;
351 }
352 }
353 if(map>=ln){
354 map=ln-1; /* guard against the approximation */
355 nextf=9999999;
356 }
357 }
358 }
359
360 /*************** vorbis decode glue ************/
361
floor0_free_info(vorbis_info_floor * i)362 void floor0_free_info(vorbis_info_floor *i){
363 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
364 if(info)_ogg_free(info);
365 }
366
floor0_info_unpack(vorbis_info * vi,oggpack_buffer * opb)367 vorbis_info_floor *floor0_info_unpack (vorbis_info *vi,oggpack_buffer *opb){
368 codec_setup_info *ci=(codec_setup_info *)vi->codec_setup;
369 int j;
370
371 vorbis_info_floor0 *info=(vorbis_info_floor0 *)_ogg_malloc(sizeof(*info));
372 info->order=oggpack_read(opb,8);
373 info->rate=oggpack_read(opb,16);
374 info->barkmap=oggpack_read(opb,16);
375 info->ampbits=oggpack_read(opb,6);
376 info->ampdB=oggpack_read(opb,8);
377 info->numbooks=oggpack_read(opb,4)+1;
378
379 if(info->order<1)goto err_out;
380 if(info->rate<1)goto err_out;
381 if(info->barkmap<1)goto err_out;
382
383 for(j=0;j<info->numbooks;j++){
384 info->books[j]=(char)oggpack_read(opb,8);
385 if(info->books[j]>=ci->books)goto err_out;
386 }
387
388 if(oggpack_eop(opb))goto err_out;
389 return(info);
390
391 err_out:
392 floor0_free_info(info);
393 return(NULL);
394 }
395
floor0_memosize(vorbis_info_floor * i)396 int floor0_memosize(vorbis_info_floor *i){
397 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
398 return info->order+1;
399 }
400
floor0_inverse1(vorbis_dsp_state * vd,vorbis_info_floor * i,ogg_int32_t * lsp)401 ogg_int32_t *floor0_inverse1(vorbis_dsp_state *vd,vorbis_info_floor *i,
402 ogg_int32_t *lsp){
403 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
404 int j,k;
405 uint64_t ampraw;
406 if(info->ampbits<=32){
407 ampraw=oggpack_read(&vd->opb,info->ampbits);
408 }else{
409 //max value possible for info->ampbits is 63
410 ampraw=oggpack_read(&vd->opb,info->ampbits-32);
411 ampraw<<=32;
412 ampraw|=oggpack_read(&vd->opb,32);
413 }
414 if(ampraw>0){ /* also handles the -1 out of data case */
415 uint64_t maxval=(1<<info->ampbits)-1;
416 float ampRatio=(float)ampraw/maxval;
417 int amp=ampRatio*(info->ampdB<<4);
418 int booknum=oggpack_read(&vd->opb,_ilog(info->numbooks));
419 if(booknum!=-1 && booknum<info->numbooks){ /* be paranoid */
420 codec_setup_info *ci=(codec_setup_info *)vd->vi->codec_setup;
421 codebook *b=ci->book_param+info->books[booknum];
422 int sz = floor0_memosize(i);
423 if (sz < b->dim) {
424 ALOGE("lsp too small: %d < %ld", sz, b->dim);
425 return NULL;
426 }
427 ogg_int32_t last=0;
428
429 for(j=0;j<info->order;j+=b->dim)
430 if(vorbis_book_decodev_set(b,lsp+j,&vd->opb,b->dim,-24)==-1)goto eop;
431 for(j=0;j<info->order;){
432 for(k=0;k<b->dim;k++,j++)lsp[j]+=last;
433 last=lsp[j-1];
434 }
435
436 lsp[info->order]=amp;
437 return(lsp);
438 }
439 }
440 eop:
441 return(NULL);
442 }
443
floor0_inverse2(vorbis_dsp_state * vd,vorbis_info_floor * i,ogg_int32_t * lsp,ogg_int32_t * out)444 int floor0_inverse2(vorbis_dsp_state *vd,vorbis_info_floor *i,
445 ogg_int32_t *lsp,ogg_int32_t *out){
446 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
447 codec_setup_info *ci=(codec_setup_info *)vd->vi->codec_setup;
448
449 if(lsp){
450 ogg_int32_t amp=lsp[info->order];
451
452 /* take the coefficients back to a spectral envelope curve */
453 vorbis_lsp_to_curve(out,ci->blocksizes[vd->W]/2,info->barkmap,
454 lsp,info->order,amp,info->ampdB,
455 info->rate>>1);
456 return(1);
457 }
458 memset(out,0,sizeof(*out)*ci->blocksizes[vd->W]/2);
459 return(0);
460 }
461
462