dbf-halloween2015

annotate libs/vorbis/lsp.c @ 1:c3f5c32cb210

find changesets by author, revision, files, or words in the commit message
barfed all the libraries in the source tree to make porting easier
author John Tsiombikas <nuclear@member.fsf.org>
date Sun, 01 Nov 2015 00:36:56 +0200
parents
children
rev   line source
nuclear@1 1 /********************************************************************
nuclear@1 2 * *
nuclear@1 3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
nuclear@1 4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
nuclear@1 5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
nuclear@1 6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
nuclear@1 7 * *
nuclear@1 8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
nuclear@1 9 * by the Xiph.Org Foundation http://www.xiph.org/ *
nuclear@1 10 * *
nuclear@1 11 ********************************************************************
nuclear@1 12
nuclear@1 13 function: LSP (also called LSF) conversion routines
nuclear@1 14 last mod: $Id: lsp.c 17538 2010-10-15 02:52:29Z tterribe $
nuclear@1 15
nuclear@1 16 The LSP generation code is taken (with minimal modification and a
nuclear@1 17 few bugfixes) from "On the Computation of the LSP Frequencies" by
nuclear@1 18 Joseph Rothweiler (see http://www.rothweiler.us for contact info).
nuclear@1 19 The paper is available at:
nuclear@1 20
nuclear@1 21 http://www.myown1.com/joe/lsf
nuclear@1 22
nuclear@1 23 ********************************************************************/
nuclear@1 24
nuclear@1 25 /* Note that the lpc-lsp conversion finds the roots of polynomial with
nuclear@1 26 an iterative root polisher (CACM algorithm 283). It *is* possible
nuclear@1 27 to confuse this algorithm into not converging; that should only
nuclear@1 28 happen with absurdly closely spaced roots (very sharp peaks in the
nuclear@1 29 LPC f response) which in turn should be impossible in our use of
nuclear@1 30 the code. If this *does* happen anyway, it's a bug in the floor
nuclear@1 31 finder; find the cause of the confusion (probably a single bin
nuclear@1 32 spike or accidental near-float-limit resolution problems) and
nuclear@1 33 correct it. */
nuclear@1 34
nuclear@1 35 #include <math.h>
nuclear@1 36 #include <string.h>
nuclear@1 37 #include <stdlib.h>
nuclear@1 38 #include "lsp.h"
nuclear@1 39 #include "os.h"
nuclear@1 40 #include "misc.h"
nuclear@1 41 #include "lookup.h"
nuclear@1 42 #include "scales.h"
nuclear@1 43
nuclear@1 44 /* three possible LSP to f curve functions; the exact computation
nuclear@1 45 (float), a lookup based float implementation, and an integer
nuclear@1 46 implementation. The float lookup is likely the optimal choice on
nuclear@1 47 any machine with an FPU. The integer implementation is *not* fixed
nuclear@1 48 point (due to the need for a large dynamic range and thus a
nuclear@1 49 separately tracked exponent) and thus much more complex than the
nuclear@1 50 relatively simple float implementations. It's mostly for future
nuclear@1 51 work on a fully fixed point implementation for processors like the
nuclear@1 52 ARM family. */
nuclear@1 53
nuclear@1 54 /* define either of these (preferably FLOAT_LOOKUP) to have faster
nuclear@1 55 but less precise implementation. */
nuclear@1 56 #undef FLOAT_LOOKUP
nuclear@1 57 #undef INT_LOOKUP
nuclear@1 58
nuclear@1 59 #ifdef FLOAT_LOOKUP
nuclear@1 60 #include "lookup.c" /* catch this in the build system; we #include for
nuclear@1 61 compilers (like gcc) that can't inline across
nuclear@1 62 modules */
nuclear@1 63
nuclear@1 64 /* side effect: changes *lsp to cosines of lsp */
nuclear@1 65 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
nuclear@1 66 float amp,float ampoffset){
nuclear@1 67 int i;
nuclear@1 68 float wdel=M_PI/ln;
nuclear@1 69 vorbis_fpu_control fpu;
nuclear@1 70
nuclear@1 71 vorbis_fpu_setround(&fpu);
nuclear@1 72 for(i=0;i<m;i++)lsp[i]=vorbis_coslook(lsp[i]);
nuclear@1 73
nuclear@1 74 i=0;
nuclear@1 75 while(i<n){
nuclear@1 76 int k=map[i];
nuclear@1 77 int qexp;
nuclear@1 78 float p=.7071067812f;
nuclear@1 79 float q=.7071067812f;
nuclear@1 80 float w=vorbis_coslook(wdel*k);
nuclear@1 81 float *ftmp=lsp;
nuclear@1 82 int c=m>>1;
nuclear@1 83
nuclear@1 84 while(c--){
nuclear@1 85 q*=ftmp[0]-w;
nuclear@1 86 p*=ftmp[1]-w;
nuclear@1 87 ftmp+=2;
nuclear@1 88 }
nuclear@1 89
nuclear@1 90 if(m&1){
nuclear@1 91 /* odd order filter; slightly assymetric */
nuclear@1 92 /* the last coefficient */
nuclear@1 93 q*=ftmp[0]-w;
nuclear@1 94 q*=q;
nuclear@1 95 p*=p*(1.f-w*w);
nuclear@1 96 }else{
nuclear@1 97 /* even order filter; still symmetric */
nuclear@1 98 q*=q*(1.f+w);
nuclear@1 99 p*=p*(1.f-w);
nuclear@1 100 }
nuclear@1 101
nuclear@1 102 q=frexp(p+q,&qexp);
nuclear@1 103 q=vorbis_fromdBlook(amp*
nuclear@1 104 vorbis_invsqlook(q)*
nuclear@1 105 vorbis_invsq2explook(qexp+m)-
nuclear@1 106 ampoffset);
nuclear@1 107
nuclear@1 108 do{
nuclear@1 109 curve[i++]*=q;
nuclear@1 110 }while(map[i]==k);
nuclear@1 111 }
nuclear@1 112 vorbis_fpu_restore(fpu);
nuclear@1 113 }
nuclear@1 114
nuclear@1 115 #else
nuclear@1 116
nuclear@1 117 #ifdef INT_LOOKUP
nuclear@1 118 #include "lookup.c" /* catch this in the build system; we #include for
nuclear@1 119 compilers (like gcc) that can't inline across
nuclear@1 120 modules */
nuclear@1 121
nuclear@1 122 static const int MLOOP_1[64]={
nuclear@1 123 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
nuclear@1 124 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
nuclear@1 125 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
nuclear@1 126 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
nuclear@1 127 };
nuclear@1 128
nuclear@1 129 static const int MLOOP_2[64]={
nuclear@1 130 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
nuclear@1 131 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
nuclear@1 132 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
nuclear@1 133 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
nuclear@1 134 };
nuclear@1 135
nuclear@1 136 static const int MLOOP_3[8]={0,1,2,2,3,3,3,3};
nuclear@1 137
nuclear@1 138
nuclear@1 139 /* side effect: changes *lsp to cosines of lsp */
nuclear@1 140 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
nuclear@1 141 float amp,float ampoffset){
nuclear@1 142
nuclear@1 143 /* 0 <= m < 256 */
nuclear@1 144
nuclear@1 145 /* set up for using all int later */
nuclear@1 146 int i;
nuclear@1 147 int ampoffseti=rint(ampoffset*4096.f);
nuclear@1 148 int ampi=rint(amp*16.f);
nuclear@1 149 long *ilsp=alloca(m*sizeof(*ilsp));
nuclear@1 150 for(i=0;i<m;i++)ilsp[i]=vorbis_coslook_i(lsp[i]/M_PI*65536.f+.5f);
nuclear@1 151
nuclear@1 152 i=0;
nuclear@1 153 while(i<n){
nuclear@1 154 int j,k=map[i];
nuclear@1 155 unsigned long pi=46341; /* 2**-.5 in 0.16 */
nuclear@1 156 unsigned long qi=46341;
nuclear@1 157 int qexp=0,shift;
nuclear@1 158 long wi=vorbis_coslook_i(k*65536/ln);
nuclear@1 159
nuclear@1 160 qi*=labs(ilsp[0]-wi);
nuclear@1 161 pi*=labs(ilsp[1]-wi);
nuclear@1 162
nuclear@1 163 for(j=3;j<m;j+=2){
nuclear@1 164 if(!(shift=MLOOP_1[(pi|qi)>>25]))
nuclear@1 165 if(!(shift=MLOOP_2[(pi|qi)>>19]))
nuclear@1 166 shift=MLOOP_3[(pi|qi)>>16];
nuclear@1 167 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
nuclear@1 168 pi=(pi>>shift)*labs(ilsp[j]-wi);
nuclear@1 169 qexp+=shift;
nuclear@1 170 }
nuclear@1 171 if(!(shift=MLOOP_1[(pi|qi)>>25]))
nuclear@1 172 if(!(shift=MLOOP_2[(pi|qi)>>19]))
nuclear@1 173 shift=MLOOP_3[(pi|qi)>>16];
nuclear@1 174
nuclear@1 175 /* pi,qi normalized collectively, both tracked using qexp */
nuclear@1 176
nuclear@1 177 if(m&1){
nuclear@1 178 /* odd order filter; slightly assymetric */
nuclear@1 179 /* the last coefficient */
nuclear@1 180 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
nuclear@1 181 pi=(pi>>shift)<<14;
nuclear@1 182 qexp+=shift;
nuclear@1 183
nuclear@1 184 if(!(shift=MLOOP_1[(pi|qi)>>25]))
nuclear@1 185 if(!(shift=MLOOP_2[(pi|qi)>>19]))
nuclear@1 186 shift=MLOOP_3[(pi|qi)>>16];
nuclear@1 187
nuclear@1 188 pi>>=shift;
nuclear@1 189 qi>>=shift;
nuclear@1 190 qexp+=shift-14*((m+1)>>1);
nuclear@1 191
nuclear@1 192 pi=((pi*pi)>>16);
nuclear@1 193 qi=((qi*qi)>>16);
nuclear@1 194 qexp=qexp*2+m;
nuclear@1 195
nuclear@1 196 pi*=(1<<14)-((wi*wi)>>14);
nuclear@1 197 qi+=pi>>14;
nuclear@1 198
nuclear@1 199 }else{
nuclear@1 200 /* even order filter; still symmetric */
nuclear@1 201
nuclear@1 202 /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
nuclear@1 203 worth tracking step by step */
nuclear@1 204
nuclear@1 205 pi>>=shift;
nuclear@1 206 qi>>=shift;
nuclear@1 207 qexp+=shift-7*m;
nuclear@1 208
nuclear@1 209 pi=((pi*pi)>>16);
nuclear@1 210 qi=((qi*qi)>>16);
nuclear@1 211 qexp=qexp*2+m;
nuclear@1 212
nuclear@1 213 pi*=(1<<14)-wi;
nuclear@1 214 qi*=(1<<14)+wi;
nuclear@1 215 qi=(qi+pi)>>14;
nuclear@1 216
nuclear@1 217 }
nuclear@1 218
nuclear@1 219
nuclear@1 220 /* we've let the normalization drift because it wasn't important;
nuclear@1 221 however, for the lookup, things must be normalized again. We
nuclear@1 222 need at most one right shift or a number of left shifts */
nuclear@1 223
nuclear@1 224 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
nuclear@1 225 qi>>=1; qexp++;
nuclear@1 226 }else
nuclear@1 227 while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
nuclear@1 228 qi<<=1; qexp--;
nuclear@1 229 }
nuclear@1 230
nuclear@1 231 amp=vorbis_fromdBlook_i(ampi* /* n.4 */
nuclear@1 232 vorbis_invsqlook_i(qi,qexp)-
nuclear@1 233 /* m.8, m+n<=8 */
nuclear@1 234 ampoffseti); /* 8.12[0] */
nuclear@1 235
nuclear@1 236 curve[i]*=amp;
nuclear@1 237 while(map[++i]==k)curve[i]*=amp;
nuclear@1 238 }
nuclear@1 239 }
nuclear@1 240
nuclear@1 241 #else
nuclear@1 242
nuclear@1 243 /* old, nonoptimized but simple version for any poor sap who needs to
nuclear@1 244 figure out what the hell this code does, or wants the other
nuclear@1 245 fraction of a dB precision */
nuclear@1 246
nuclear@1 247 /* side effect: changes *lsp to cosines of lsp */
nuclear@1 248 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
nuclear@1 249 float amp,float ampoffset){
nuclear@1 250 int i;
nuclear@1 251 float wdel=M_PI/ln;
nuclear@1 252 for(i=0;i<m;i++)lsp[i]=2.f*cos(lsp[i]);
nuclear@1 253
nuclear@1 254 i=0;
nuclear@1 255 while(i<n){
nuclear@1 256 int j,k=map[i];
nuclear@1 257 float p=.5f;
nuclear@1 258 float q=.5f;
nuclear@1 259 float w=2.f*cos(wdel*k);
nuclear@1 260 for(j=1;j<m;j+=2){
nuclear@1 261 q *= w-lsp[j-1];
nuclear@1 262 p *= w-lsp[j];
nuclear@1 263 }
nuclear@1 264 if(j==m){
nuclear@1 265 /* odd order filter; slightly assymetric */
nuclear@1 266 /* the last coefficient */
nuclear@1 267 q*=w-lsp[j-1];
nuclear@1 268 p*=p*(4.f-w*w);
nuclear@1 269 q*=q;
nuclear@1 270 }else{
nuclear@1 271 /* even order filter; still symmetric */
nuclear@1 272 p*=p*(2.f-w);
nuclear@1 273 q*=q*(2.f+w);
nuclear@1 274 }
nuclear@1 275
nuclear@1 276 q=fromdB(amp/sqrt(p+q)-ampoffset);
nuclear@1 277
nuclear@1 278 curve[i]*=q;
nuclear@1 279 while(map[++i]==k)curve[i]*=q;
nuclear@1 280 }
nuclear@1 281 }
nuclear@1 282
nuclear@1 283 #endif
nuclear@1 284 #endif
nuclear@1 285
nuclear@1 286 static void cheby(float *g, int ord) {
nuclear@1 287 int i, j;
nuclear@1 288
nuclear@1 289 g[0] *= .5f;
nuclear@1 290 for(i=2; i<= ord; i++) {
nuclear@1 291 for(j=ord; j >= i; j--) {
nuclear@1 292 g[j-2] -= g[j];
nuclear@1 293 g[j] += g[j];
nuclear@1 294 }
nuclear@1 295 }
nuclear@1 296 }
nuclear@1 297
nuclear@1 298 static int comp(const void *a,const void *b){
nuclear@1 299 return (*(float *)a<*(float *)b)-(*(float *)a>*(float *)b);
nuclear@1 300 }
nuclear@1 301
nuclear@1 302 /* Newton-Raphson-Maehly actually functioned as a decent root finder,
nuclear@1 303 but there are root sets for which it gets into limit cycles
nuclear@1 304 (exacerbated by zero suppression) and fails. We can't afford to
nuclear@1 305 fail, even if the failure is 1 in 100,000,000, so we now use
nuclear@1 306 Laguerre and later polish with Newton-Raphson (which can then
nuclear@1 307 afford to fail) */
nuclear@1 308
nuclear@1 309 #define EPSILON 10e-7
nuclear@1 310 static int Laguerre_With_Deflation(float *a,int ord,float *r){
nuclear@1 311 int i,m;
nuclear@1 312 double *defl=alloca(sizeof(*defl)*(ord+1));
nuclear@1 313 for(i=0;i<=ord;i++)defl[i]=a[i];
nuclear@1 314
nuclear@1 315 for(m=ord;m>0;m--){
nuclear@1 316 double new=0.f,delta;
nuclear@1 317
nuclear@1 318 /* iterate a root */
nuclear@1 319 while(1){
nuclear@1 320 double p=defl[m],pp=0.f,ppp=0.f,denom;
nuclear@1 321
nuclear@1 322 /* eval the polynomial and its first two derivatives */
nuclear@1 323 for(i=m;i>0;i--){
nuclear@1 324 ppp = new*ppp + pp;
nuclear@1 325 pp = new*pp + p;
nuclear@1 326 p = new*p + defl[i-1];
nuclear@1 327 }
nuclear@1 328
nuclear@1 329 /* Laguerre's method */
nuclear@1 330 denom=(m-1) * ((m-1)*pp*pp - m*p*ppp);
nuclear@1 331 if(denom<0)
nuclear@1 332 return(-1); /* complex root! The LPC generator handed us a bad filter */
nuclear@1 333
nuclear@1 334 if(pp>0){
nuclear@1 335 denom = pp + sqrt(denom);
nuclear@1 336 if(denom<EPSILON)denom=EPSILON;
nuclear@1 337 }else{
nuclear@1 338 denom = pp - sqrt(denom);
nuclear@1 339 if(denom>-(EPSILON))denom=-(EPSILON);
nuclear@1 340 }
nuclear@1 341
nuclear@1 342 delta = m*p/denom;
nuclear@1 343 new -= delta;
nuclear@1 344
nuclear@1 345 if(delta<0.f)delta*=-1;
nuclear@1 346
nuclear@1 347 if(fabs(delta/new)<10e-12)break;
nuclear@1 348 }
nuclear@1 349
nuclear@1 350 r[m-1]=new;
nuclear@1 351
nuclear@1 352 /* forward deflation */
nuclear@1 353
nuclear@1 354 for(i=m;i>0;i--)
nuclear@1 355 defl[i-1]+=new*defl[i];
nuclear@1 356 defl++;
nuclear@1 357
nuclear@1 358 }
nuclear@1 359 return(0);
nuclear@1 360 }
nuclear@1 361
nuclear@1 362
nuclear@1 363 /* for spit-and-polish only */
nuclear@1 364 static int Newton_Raphson(float *a,int ord,float *r){
nuclear@1 365 int i, k, count=0;
nuclear@1 366 double error=1.f;
nuclear@1 367 double *root=alloca(ord*sizeof(*root));
nuclear@1 368
nuclear@1 369 for(i=0; i<ord;i++) root[i] = r[i];
nuclear@1 370
nuclear@1 371 while(error>1e-20){
nuclear@1 372 error=0;
nuclear@1 373
nuclear@1 374 for(i=0; i<ord; i++) { /* Update each point. */
nuclear@1 375 double pp=0.,delta;
nuclear@1 376 double rooti=root[i];
nuclear@1 377 double p=a[ord];
nuclear@1 378 for(k=ord-1; k>= 0; k--) {
nuclear@1 379
nuclear@1 380 pp= pp* rooti + p;
nuclear@1 381 p = p * rooti + a[k];
nuclear@1 382 }
nuclear@1 383
nuclear@1 384 delta = p/pp;
nuclear@1 385 root[i] -= delta;
nuclear@1 386 error+= delta*delta;
nuclear@1 387 }
nuclear@1 388
nuclear@1 389 if(count>40)return(-1);
nuclear@1 390
nuclear@1 391 count++;
nuclear@1 392 }
nuclear@1 393
nuclear@1 394 /* Replaced the original bubble sort with a real sort. With your
nuclear@1 395 help, we can eliminate the bubble sort in our lifetime. --Monty */
nuclear@1 396
nuclear@1 397 for(i=0; i<ord;i++) r[i] = root[i];
nuclear@1 398 return(0);
nuclear@1 399 }
nuclear@1 400
nuclear@1 401
nuclear@1 402 /* Convert lpc coefficients to lsp coefficients */
nuclear@1 403 int vorbis_lpc_to_lsp(float *lpc,float *lsp,int m){
nuclear@1 404 int order2=(m+1)>>1;
nuclear@1 405 int g1_order,g2_order;
nuclear@1 406 float *g1=alloca(sizeof(*g1)*(order2+1));
nuclear@1 407 float *g2=alloca(sizeof(*g2)*(order2+1));
nuclear@1 408 float *g1r=alloca(sizeof(*g1r)*(order2+1));
nuclear@1 409 float *g2r=alloca(sizeof(*g2r)*(order2+1));
nuclear@1 410 int i;
nuclear@1 411
nuclear@1 412 /* even and odd are slightly different base cases */
nuclear@1 413 g1_order=(m+1)>>1;
nuclear@1 414 g2_order=(m) >>1;
nuclear@1 415
nuclear@1 416 /* Compute the lengths of the x polynomials. */
nuclear@1 417 /* Compute the first half of K & R F1 & F2 polynomials. */
nuclear@1 418 /* Compute half of the symmetric and antisymmetric polynomials. */
nuclear@1 419 /* Remove the roots at +1 and -1. */
nuclear@1 420
nuclear@1 421 g1[g1_order] = 1.f;
nuclear@1 422 for(i=1;i<=g1_order;i++) g1[g1_order-i] = lpc[i-1]+lpc[m-i];
nuclear@1 423 g2[g2_order] = 1.f;
nuclear@1 424 for(i=1;i<=g2_order;i++) g2[g2_order-i] = lpc[i-1]-lpc[m-i];
nuclear@1 425
nuclear@1 426 if(g1_order>g2_order){
nuclear@1 427 for(i=2; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+2];
nuclear@1 428 }else{
nuclear@1 429 for(i=1; i<=g1_order;i++) g1[g1_order-i] -= g1[g1_order-i+1];
nuclear@1 430 for(i=1; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+1];
nuclear@1 431 }
nuclear@1 432
nuclear@1 433 /* Convert into polynomials in cos(alpha) */
nuclear@1 434 cheby(g1,g1_order);
nuclear@1 435 cheby(g2,g2_order);
nuclear@1 436
nuclear@1 437 /* Find the roots of the 2 even polynomials.*/
nuclear@1 438 if(Laguerre_With_Deflation(g1,g1_order,g1r) ||
nuclear@1 439 Laguerre_With_Deflation(g2,g2_order,g2r))
nuclear@1 440 return(-1);
nuclear@1 441
nuclear@1 442 Newton_Raphson(g1,g1_order,g1r); /* if it fails, it leaves g1r alone */
nuclear@1 443 Newton_Raphson(g2,g2_order,g2r); /* if it fails, it leaves g2r alone */
nuclear@1 444
nuclear@1 445 qsort(g1r,g1_order,sizeof(*g1r),comp);
nuclear@1 446 qsort(g2r,g2_order,sizeof(*g2r),comp);
nuclear@1 447
nuclear@1 448 for(i=0;i<g1_order;i++)
nuclear@1 449 lsp[i*2] = acos(g1r[i]);
nuclear@1 450
nuclear@1 451 for(i=0;i<g2_order;i++)
nuclear@1 452 lsp[i*2+1] = acos(g2r[i]);
nuclear@1 453 return(0);
nuclear@1 454 }