diff options
author | Monty <xiphmont@xiph.org> | 2001-04-29 22:21:06 +0000 |
---|---|---|
committer | Monty <xiphmont@xiph.org> | 2001-04-29 22:21:06 +0000 |
commit | 2b1605de070613fa21f435844d9c5ea8a3688f19 (patch) | |
tree | 26c9fdaad52bcd1c3e223128948753f711761f7d | |
parent | 78540c3db714ab4dda28e11659f2a69e9961bf69 (diff) | |
download | libvorbis-git-2b1605de070613fa21f435844d9c5ea8a3688f19.tar.gz |
...such that I don;t lose anything (floor 1 not yet tested/functional)
Monty
svn path=/branches/monty-branch-20010404/vorbis/; revision=1435
-rw-r--r-- | lib/backends.h | 32 | ||||
-rw-r--r-- | lib/barkmel.c | 64 | ||||
-rw-r--r-- | lib/floor0.c | 10 | ||||
-rw-r--r-- | lib/floor1.c | 949 | ||||
-rw-r--r-- | lib/mapping0.c | 401 | ||||
-rw-r--r-- | lib/modes/mode_A.h | 316 | ||||
-rw-r--r-- | lib/psy.c | 727 | ||||
-rw-r--r-- | lib/psy.h | 101 |
8 files changed, 2477 insertions, 123 deletions
diff --git a/lib/backends.h b/lib/backends.h index 72097000..8c64ca41 100644 --- a/lib/backends.h +++ b/lib/backends.h @@ -12,7 +12,7 @@ function: libvorbis backend and mapping structures; needed for static mode headers - last mod: $Id: backends.h,v 1.6.4.1 2001/04/05 00:22:48 xiphmont Exp $ + last mod: $Id: backends.h,v 1.6.4.2 2001/04/29 22:21:04 xiphmont Exp $ ********************************************************************/ @@ -61,7 +61,7 @@ typedef struct{ void (*free_info) (vorbis_info_floor *); void (*free_look) (vorbis_look_floor *); int (*forward) (struct vorbis_block *,vorbis_look_floor *, - float *, float *); + float *); int (*inverse) (struct vorbis_block *,vorbis_look_floor *, float *); } vorbis_func_floor; @@ -82,21 +82,27 @@ typedef struct{ } vorbis_info_floor0; +#define VIF_POSIT 63 +#define VIF_CLASS 16 +#define VIF_PARTS 31 typedef struct{ - long maxrange; - int rangeres; /* 6, 7 or 8 */ + int partitions; /* 0 to 31 */ + int partitionclass[VIF_PARTS]; /* 0 to 15 */ - int ranges; /* up to 8 */ - int rangesizes[8]; /* 1 to 8 */ - int rangelist[64]; - - int positionbook[8]; - int ampbook[8]; + int class_dim[VIF_CLASS]; /* 1 to 8 */ + int class_subs[VIF_CLASS]; /* 0,1,2,3 (bits: 1<<n poss) */ + int class_book[VIF_CLASS]; /* subs ^ dim entries */ + int class_subbook[VIF_CLASS][8]; /* [VIF_CLASS][subs] */ + + + int mult; /* 1 2 3 or 4 */ + int postlist[VIF_POSIT+2]; /* first two implicit */ - float maxover; /* encode side analysis parameter */ - float maxunder; /* encode side analysis parameter */ - float maxerr; /* encode side analysis parameter */ + float maxover; /* encode side analysis parameter */ + float maxunder; /* encode side analysis parameter */ + float maxerr; /* encode side analysis parameter */ + int searchstart; /* encode side analysis parameter */ } vorbis_info_floor1; /* Residue backend generic *****************************************/ diff --git a/lib/barkmel.c b/lib/barkmel.c new file mode 100644 index 00000000..1f3ba38c --- /dev/null +++ b/lib/barkmel.c @@ -0,0 +1,64 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * + * by the XIPHOPHORUS Company http://www.xiph.org/ * + + ******************************************************************** + + function: bark scale utility + last mod: $Id: barkmel.c,v 1.5.4.1 2001/04/29 22:21:04 xiphmont Exp $ + + ********************************************************************/ + +#include <stdio.h> +#include "scales.h" +int main(){ + int i; + double rate; + for(i=64;i<32000;i*=2){ + rate=48000.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=44100.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=32000.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=22050.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=16000.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=11025.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + rate=8000.f; + fprintf(stderr,"rate=%gHz, block=%d, f(1)=%.2gHz bark(1)=%.2g (of %.2g)\n\n", + rate,i,rate/2 / (i/2),toBARK(rate/2 /(i/2)),toBARK(rate/2)); + + + } + { + float i; + int j; + for(i=0.,j=0;i<28;i+=.816,j++){ + fprintf(stderr,"(%d) bark=%f %gHz (%d of 1024)\n", + j,i,fromBARK(i),(int)(fromBARK(i)/22050.*1024.)); + } + } + return(0); +} + diff --git a/lib/floor0.c b/lib/floor0.c index 6be7fa90..e2f443f1 100644 --- a/lib/floor0.c +++ b/lib/floor0.c @@ -11,7 +11,7 @@ ******************************************************************** function: floor backend 0 implementation - last mod: $Id: floor0.c,v 1.39.4.1 2001/04/05 00:22:48 xiphmont Exp $ + last mod: $Id: floor0.c,v 1.39.4.2 2001/04/29 22:21:04 xiphmont Exp $ ********************************************************************/ @@ -240,8 +240,8 @@ float _curve_to_lpc(float *curve,float *lpc, for(i=0;i<l->n;i++) curve[i]-=150; - _analysis_output_always("barkfloor",seq,work,bark,0,0); - _analysis_output_always("barkcurve",seq++,curve,l->n,1,0); + _analysis_output("barkfloor",seq,work,bark,0,0); + _analysis_output("barkcurve",seq++,curve,l->n,1,0); for(i=0;i<l->n;i++) curve[i]+=150; @@ -255,7 +255,7 @@ float _curve_to_lpc(float *curve,float *lpc, /* didn't need in->out seperation, modifies the flr[] vector; takes in a dB scale floor, puts out linear */ static int floor0_forward(vorbis_block *vb,vorbis_look_floor *i, - float *flr,float *dummy){ + float *flr){ long j; vorbis_look_floor0 *look=(vorbis_look_floor0 *)i; vorbis_info_floor0 *info=look->vi; @@ -385,7 +385,7 @@ static int floor0_forward(vorbis_block *vb,vorbis_look_floor *i, vorbis_lsp_to_curve(flr,look->linearmap,look->n,look->ln, lspwork,look->m,amp,info->ampdB); - _analysis_output_always("barklsp",seq-1,flr,look->n,1,1); + _analysis_output("barklsp",seq-1,flr,look->n,1,1); _analysis_output("lsp3",seq-1,flr,look->n,0,1); return(val); } diff --git a/lib/floor1.c b/lib/floor1.c index 6457f506..84128f80 100644 --- a/lib/floor1.c +++ b/lib/floor1.c @@ -11,7 +11,7 @@ ******************************************************************** function: floor backend 1 implementation - last mod: $Id: floor1.c,v 1.1.2.1 2001/04/05 00:22:48 xiphmont Exp $ + last mod: $Id: floor1.c,v 1.1.2.2 2001/04/29 22:21:04 xiphmont Exp $ ********************************************************************/ @@ -24,21 +24,43 @@ #include "registry.h" #include "codebook.h" #include "misc.h" +#include "scales.h" #include <stdio.h> +#define VORBIS_IEEE_FLOAT32 + #define floor1_rangedB 140 /* floor 0 fixed at -140dB to 0dB range */ typedef struct { - int lo_neighbor[66]; - int hi_neighbor[66]; + int sorted_index[VIF_POSIT+2]; /* en/de */ + int forward_index[VIF_POSIT+2]; /* en/de */ + int reverse_index[VIF_POSIT+2]; /* en */ + + int hineighbor[VIF_POSIT]; /* de */ + int loneighbor[VIF_POSIT]; /* de */ + int posts; + int n; + int quant_q; vorbis_info_floor1 *vi; } vorbis_look_floor1; -/***********************************************/ +typedef struct lsfit_acc{ + long x0; + long x1; + long xa; + long ya; + long x2a; + long y2a; + long xya; + long n; +} lsfit_acc; + +/***********************************************/ + static vorbis_info_floor *floor1_copy_info (vorbis_info_floor *i){ vorbis_info_floor1 *info=(vorbis_info_floor1 *)i; vorbis_info_floor1 *ret=_ogg_malloc(sizeof(vorbis_info_floor1)); @@ -61,6 +83,15 @@ static void floor1_free_look(vorbis_look_floor *i){ } } +static int ilog(unsigned int v){ + int ret=0; + while(v){ + ret++; + v>>=1; + } + return(ret); +} + static int ilog2(unsigned int v){ int ret=0; while(v>1){ @@ -75,73 +106,103 @@ static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){ int j,k; int count=0; int rangebits; + int maxposit=info->postlist[1]; + int maxclass=-1; - oggpack_write(opb,info->maxrange,24); - rangebits=ilog2(info->maxrange); - - oggpack_write(opb,info->ranges,4); /* only 0 to 8 legal now */ - oggpack_write(opb,info->rangeres-1,3); /* only 6,7 or 8 legal now */ - for(j=0,k=0;j<info->ranges;j++){ - oggpack_write(opb,info->rangesizes[j]-1,3); - oggpack_write(opb,info->positionbook[j],8); - oggpack_write(opb,info->ampbook[j],8); - count+=info->rangesizes[j]; + /* save out partitions */ + oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */ + for(j=0;j<info->partitions;j++){ + oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */ + if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j]; + } + + /* save out partition classes */ + for(j=0;j<maxclass+1;j++){ + oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */ + oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */ + if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8); + for(k=0;k<info->class_subs[j];k++) + oggpack_write(opb,info->class_subbook[j][k]+1,8); + } + + /* save out the post list */ + oggpack_write(opb,info->mult-1,2); /* only 1,2,3,4 legal now */ + oggpack_write(opb,ilog2(maxposit),4); + rangebits=ilog2(maxposit); + + for(j=0,k=0;j<info->partitions;j++){ + count+=info->class_dim[info->partitionclass[j]]; for(;k<count;k++) - oggpack_write(opb,info->rangelist[k],rangebits); + oggpack_write(opb,info->postlist[k+2],rangebits); } } + static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){ codec_setup_info *ci=vi->codec_setup; - int j,k,count; - - vorbis_info_floor1 *info=_ogg_malloc(sizeof(vorbis_info_floor1)); - int rangebits; - - info->maxrange=oggpack_read(opb,24); - rangebits=ilog2(info->maxrange); - - info->ranges=oggpack_read(opb,4); - if(info->ranges>8)goto err_out; + int j,k,count,maxclass=-1,rangebits; - info->rangeres=oggpack_read(opb,3)+1; - if(info->rangeres<6)goto err_out; - if(info->rangeres>8)goto err_out; + vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(vorbis_info_floor1)); + /* read partitions */ + info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */ + for(j=0;j<info->partitions;j++){ + info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */ + if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j]; + } - for(j=0,k=0;j<info->ranges;j++){ - info->rangesizes[j]=oggpack_read(opb,3)+1; - info->positionbook[j]=oggpack_read(opb,8); - if(info->positionbook[j]<0 || info->positionbook[j]>=ci->books) + /* read partition classes */ + for(j=0;j<maxclass+1;j++){ + info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */ + info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */ + if(info->class_subs[j]<0) goto err_out; - - info->ampbook[j]=oggpack_read(opb,8); - if(info->ampbook[j]<0 || info->ampbook[j]>=ci->books) + if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8); + if(info->class_book[j]<0 || info->class_book[j]>=ci->books) goto err_out; - count+=info->rangesizes[j]; + for(k=0;k<info->class_subs[j];k++){ + info->class_subbook[j][k]=oggpack_read(opb,8)-1; + if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books) + goto err_out; + } + } + + /* read the post list */ + info->mult=oggpack_read(opb,2)+1; /* only 1,2,3,4 legal now */ + rangebits=oggpack_read(opb,4); + + for(j=0,k=0;j<info->partitions;j++){ + count+=info->class_dim[info->partitionclass[j]]; for(;k<count;k++){ - info->rangelist[k]=oggpack_read(opb,rangebits); - if(info->rangelist[k]<0 || info->rangelist[k]>=maxrange) + int t=info->postlist[k+2]=oggpack_read(opb,rangebits); + if(t<0 || t>=(1<<rangebits)) goto err_out; } } + info->postlist[0]=0; + info->postlist[1]=1<<rangebits; + + return(info); err_out: floor1_free_info(info); return(NULL); } +static int icomp(const void *a,const void *b){ + return(**(int **)a-**(int **)b); +} + static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,vorbis_info_mode *mi, - vorbis_info_floor *i){ + vorbis_info_floor *in){ - vorbis_info *vi=vd->vi; - codec_setup_info *ci=vi->codec_setup; - vorbis_info_floor1 *info=(vorbis_info_floor1 *)i; + int *sortpointer[VIF_POSIT+2]; + vorbis_info_floor1 *info=(vorbis_info_floor1 *)in; vorbis_look_floor1 *look=_ogg_calloc(1,sizeof(vorbis_look_floor1)); int i,j,n=0; - look->vi=vi; - look->n=ci->blocksizes[mi->blockflag]/2; - + look->vi=info; + look->n=info->postlist[1]; + /* we drop each position value in-between already decoded values, and use linear interpolation to predict each new value past the edges. The positions are read in the order of the position @@ -149,37 +210,232 @@ static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,vorbis_info_mode *mi, course, the neighbors can change (if a position is declined), but this is an initial mapping */ - for(i=0;i<info->ranges;i++)n+=info->rangesizes[j]; - for(i=0;i<info->ranges;i++){ - int hi=-1; - int lo=-1; - for(j=0;j<i;j++){ - if(info->rangelist[j]<info->rangelist[i] && - (lo==-1 || info->rangelist[j]>info->rangelist[lo])) + for(i=0;i<info->partitions;i++)n+=info->class_dim[info->partitionclass[i]]; + n+=2; + look->posts=n; + + /* also store a sorted position index */ + for(i=0;i<n;i++)sortpointer[i]=info->postlist+i; + qsort(sortpointer,n,sizeof(int),icomp); + + /* points from sort order back to range number */ + for(i=0;i<n;i++)look->forward_index[i]=sortpointer[i]-info->postlist; + /* points from range order to sorted position */ + for(i=0;i<n;i++)look->reverse_index[look->forward_index[i]]=i; + /* we actually need the post values too */ + for(i=0;i<n;i++)look->sorted_index[i]=info->postlist[look->forward_index[i]]; + + /* quantize values to multiplier spec */ + switch(info->mult){ + case 1: /* 1024 -> 256 */ + look->quant_q=256; + break; + case 2: /* 1024 -> 128 */ + look->quant_q=128; + break; + case 3: /* 1024 -> 86 */ + look->quant_q=86; + break; + case 4: /* 1024 -> 64 */ + look->quant_q=64; + break; + } + + /* discover our neighbors for decode where we don't use fit flags + (that would push the neighbors outward) */ + for(i=0;i<n-2;i++){ + int lo=0; + int hi=1; + int lx=0; + int hx=look->n; + int currentx=info->postlist[i+2]; + for(j=0;j<i+2;j++){ + int x=info->postlist[j]; + if(x>lx && x<currentx){ lo=j; - if(info->rangelist[j]>info->rangelist[i] && - (hi==-1 || info->rangelist[j]<info->rangelist[hi])) + lx=x; + } + if(x<hx && x>currentx){ hi=j; + hx=x; + } } - look->lo_neighbor[i]=lo; - look->hi_neighbor[i]=hi; + look->loneighbor[i]=lo; + look->hineighbor[i]=hi; } return(look); } -static int render_point_dB(int x0,int x1,int y0,int y1,int x){ +static int render_point(int x0,int x1,int y0,int y1,int x){ int dy=y1-y0; int adx=x1-x0; int ady=abs(dy); - int err=(ady>>1)+ady*(x-x0); + int err=ady*(x-x0); - int off=err%adx; + int off=err/adx; if(dy<0)return(y0-off); return(y0+off); } -static int render_line_dB(int x0,int x1,int y0,int y1,int *d){ +static int FLOOR_todB_LOOKUP[560]={ + 1023, 1021, 1019, 1018, 1016, 1014, 1012, 1010, + 1008, 1007, 1005, 1003, 1001, 999, 997, 996, + 994, 992, 990, 988, 986, 985, 983, 981, + 979, 977, 975, 974, 972, 970, 968, 966, + 964, 963, 961, 959, 957, 955, 954, 952, + 950, 948, 946, 944, 943, 941, 939, 937, + 935, 933, 932, 930, 928, 926, 924, 922, + 921, 919, 917, 915, 913, 911, 910, 908, + 906, 904, 902, 900, 899, 897, 895, 893, + 891, 890, 888, 886, 884, 882, 880, 879, + 877, 875, 873, 871, 869, 868, 866, 864, + 862, 860, 858, 857, 855, 853, 851, 849, + 847, 846, 844, 842, 840, 838, 836, 835, + 833, 831, 829, 827, 826, 824, 822, 820, + 818, 816, 815, 813, 811, 809, 807, 805, + 804, 802, 800, 798, 796, 794, 793, 791, + 789, 787, 785, 783, 782, 780, 778, 776, + 774, 772, 771, 769, 767, 765, 763, 762, + 760, 758, 756, 754, 752, 751, 749, 747, + 745, 743, 741, 740, 738, 736, 734, 732, + 730, 729, 727, 725, 723, 721, 719, 718, + 716, 714, 712, 710, 708, 707, 705, 703, + 701, 699, 698, 696, 694, 692, 690, 688, + 687, 685, 683, 681, 679, 677, 676, 674, + 672, 670, 668, 666, 665, 663, 661, 659, + 657, 655, 654, 652, 650, 648, 646, 644, + 643, 641, 639, 637, 635, 634, 632, 630, + 628, 626, 624, 623, 621, 619, 617, 615, + 613, 612, 610, 608, 606, 604, 602, 601, + 599, 597, 595, 593, 591, 590, 588, 586, + 584, 582, 580, 579, 577, 575, 573, 571, + 570, 568, 566, 564, 562, 560, 559, 557, + 555, 553, 551, 549, 548, 546, 544, 542, + 540, 538, 537, 535, 533, 531, 529, 527, + 526, 524, 522, 520, 518, 516, 515, 513, + 511, 509, 507, 506, 504, 502, 500, 498, + 496, 495, 493, 491, 489, 487, 485, 484, + 482, 480, 478, 476, 474, 473, 471, 469, + 467, 465, 463, 462, 460, 458, 456, 454, + 452, 451, 449, 447, 445, 443, 442, 440, + 438, 436, 434, 432, 431, 429, 427, 425, + 423, 421, 420, 418, 416, 414, 412, 410, + 409, 407, 405, 403, 401, 399, 398, 396, + 394, 392, 390, 388, 387, 385, 383, 381, + 379, 378, 376, 374, 372, 370, 368, 367, + 365, 363, 361, 359, 357, 356, 354, 352, + 350, 348, 346, 345, 343, 341, 339, 337, + 335, 334, 332, 330, 328, 326, 324, 323, + 321, 319, 317, 315, 314, 312, 310, 308, + 306, 304, 303, 301, 299, 297, 295, 293, + 292, 290, 288, 286, 284, 282, 281, 279, + 277, 275, 273, 271, 270, 268, 266, 264, + 262, 260, 259, 257, 255, 253, 251, 250, + 248, 246, 244, 242, 240, 239, 237, 235, + 233, 231, 229, 228, 226, 224, 222, 220, + 218, 217, 215, 213, 211, 209, 207, 206, + 204, 202, 200, 198, 196, 195, 193, 191, + 189, 187, 186, 184, 182, 180, 178, 176, + 175, 173, 171, 169, 167, 165, 164, 162, + 160, 158, 156, 154, 153, 151, 149, 147, + 145, 143, 142, 140, 138, 136, 134, 132, + 131, 129, 127, 125, 123, 122, 120, 118, + 116, 114, 112, 111, 109, 107, 105, 103, + 101, 100, 98, 96, 94, 92, 90, 89, + 87, 85, 83, 81, 79, 78, 76, 74, + 72, 70, 68, 67, 65, 63, 61, 59, + 58, 56, 54, 52, 50, 48, 47, 45, + 43, 41, 39, 37, 36, 34, 32, 30, + 28, 26, 25, 23, 21, 19, 17, 15, + 14, 12, 10, 8, 6, 4, 3, 1, +}; + +static float FLOOR_fromdB_LOOKUP[256]={ + 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F, + 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F, + 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F, + 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F, + 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F, + 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F, + 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F, + 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F, + 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F, + 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F, + 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F, + 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F, + 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F, + 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F, + 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F, + 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F, + 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F, + 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F, + 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F, + 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F, + 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F, + 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F, + 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F, + 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F, + 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F, + 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F, + 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F, + 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F, + 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F, + 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F, + 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F, + 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F, + 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F, + 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F, + 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F, + 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F, + 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F, + 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F, + 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F, + 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F, + 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F, + 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F, + 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F, + 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F, + 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F, + 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F, + 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F, + 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F, + 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F, + 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F, + 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F, + 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F, + 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F, + 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F, + 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F, + 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F, + 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F, + 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F, + 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F, + 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F, + 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F, + 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F, + 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F, + 0.82788260F, 0.88168307F, 0.9389798F, 1.F, +}; + +#ifdef VORBIS_IEEE_FLOAT32 +static int vorbis_floor1_dBquant(float *x){ + float temp=*x-256.f; + ogg_uint32_t *i=(ogg_uint32_t *)(&temp); + if(*i<(ogg_uint32_t)0xc3800000)return(1023); + if(*i>(ogg_uint32_t)0xc3c5e000)return(0); + return FLOOR_quantdB_LOOKUP[(*i-0xc3800000)>>13]; +} +#else +static int vorbis_floor1_dBquant(float *x){ + int i= ((*x)+140.)/140.*1024.+.5; + if(i>1023)return(1023); + if(i<0)return(0); + return i; +} +#endif + +static void render_line(int x0,int x1,int y0,int y1,float *d){ int dy=y1-y0; int adx=x1-x0; int ady=abs(dy); @@ -187,95 +443,579 @@ static int render_line_dB(int x0,int x1,int y0,int y1,int *d){ int sy=(dy<0?base-1:base+1); int x=x0; int y=y0; - int err; + int err=0; - ady-=base*adx; - err=ady>>1; + ady-=abs(base*adx); - d[x]=y; + d[x]=FLOOR_fromdB_LOOKUP[y]; while(++x<x1){ err=err+ady; - if(err>adx){ + if(err>=adx){ err-=adx; y+=sy; }else{ y+=base; } - d[x]=y; + d[x]=FLOOR_fromdB_LOOKUP[y]; } } /* the floor has already been filtered to only include relevant sections */ -static int fit_line(int *floor,int x0, int x1,int *y0, int *y1){ - long i,n=0; - long x=0,y=0,x2=0,y2=0,xy=0; - - if(*y0>-900){ /* hint used to break degenerate cases */ - x+=x0; - y+=*y0; - x2+=(x0*x0); +static int accumulate_fit(float *floor,int x0, int x1,lsfit_acc *a){ + long i; + + memset(a,0,sizeof(lsfit_acc)); + a->x0=x0; + a->x1=x1; + + for(i=x0;i<x1;i++){ + int quantized=vorbis_floor1_dBquant(floor+i); + if(quantized){ + a->xa += i; + a->ya += quantized; + a->x2a += (i*i); + a->y2a += quantized*quantized; + a->xya += i*quantized; + a->n++; + } + } + return(a->n); +} + +/* returns < 0 on too few points to fit, >=0 (meansq error) on success */ +static int fit_line(lsfit_acc *a,int fits,int *y0,int *y1){ + long x=0,y=0,x2=0,y2=0,xy=0,n=0,i; + long x0=a[0].x0; + long x1=a[fits-1].x1; + + if(*y0>=0){ /* hint used to break degenerate cases */ + x+= x0; + y+= *y0; + x2+= x0 * x0; y2+= *y0 * *y0; - xy= *y0*x0; + xy+= *y0 * x0; n++; } - if(*y1>-900){ /* hint used to break degenerate cases */ - x+=x1; - y+=*y1; - x2+=(x1*x1); + if(*y1>=0){ /* hint used to break degenerate cases */ + x+= x1; + y+= *y1; + x2+= x1 * x1; y2+= *y1 * *y1; - xy= *y1*x1; + xy+= *y1 * x1; n++; } - for(i=x0;i<x1;i++) - if(floor[i]>-900){ - x+=i; - y+=floor[i]; - x2+=(i*i); - y2+=floor[i]*floor[i]; - xy=i*floor[i]; - n++; - } + for(i=0;i<fits;i++){ + x+=a[i].xa; + y+=a[i].ya; + x2+=a[i].x2a; + y2+=a[i].y2a; + xy+=a[i].xya; + n+=a[i].n; + } if(n<2)return(-1); { - float a=(float)(y*x2-xy*x)/(n*x2-x*x); - float b=(float)(n*xy-x*y)/(n*x2-x*x); + /* need 64 bit multiplies, which C doesn't give portably as int */ + double fx=x; + double fy=y; + double fx2=x2; + double fy2=y2; + double fxy=xy; + double a=(fy*fx2-fxy*fx)/(n*fx2-fx*fx); + double b=(n*fxy-fx*fy)/(n*fx2-fx*fx); + int s=rint((a*a*n + a*b*(2*fx) - a*(2*fy) + b*b*fx2 - b*(2*fxy) + fy2)/(n*n)); + if(s<0)s=0; *y0=rint(a+b*x0); *y1=rint(a+b*x1); + return(s); } +} + +static int inspect_error(int x0,int x1,int y0,int y1,float *flr, + vorbis_info_floor1 *info){ + int dy=y1-y0; + int adx=x1-x0; + int ady=abs(dy); + int base=dy/adx; + int sy=(dy<0?base-1:base+1); + int x=x0; + int y=y0; + int err=0; + int mse=0; + int val=vorbis_floor1_dBquant(flr+x); + int n=0; + + ady-=abs(base*adx); + + if(val){ + if(y-info->maxover>val)return(1); + if(y+info->maxunder<val)return(1); + mse=(y-val); + mse*=mse; + n++; + } + + while(++x<x1){ + err=err+ady; + if(err>=adx){ + err-=adx; + y+=sy; + }else{ + y+=base; + } + + val=vorbis_floor1_dBquant(flr+x); + if(val){ + if(y-info->maxover>val)return(1); + if(y+info->maxunder<val)return(1); + mse+=((y-val)*(y-val)); + n++; + } + } + + if(mse/n>info->maxerr)return(1); return(0); } +int post_Y(int *A,int *B,int pos){ + if(A[pos]==-1) + return B[pos]; + if(B[pos]==-1) + return A[pos]; + return (A[pos]+B[pos])>>1; +} /* didn't need in->out seperation, modifies the flr[] vector; takes in a dB scale floor, puts out linear */ -static int floor1_forward(vorbis_block *vb,vorbis_look_floor *i, - float *flr,float *mdct){ - long i,j; - vorbis_look_floor1 *look=(vorbis_look_floor1 *)i; +static int floor1_forward(vorbis_block *vb,vorbis_look_floor *in, + float *flr){ + long i,j,k,l; + vorbis_look_floor1 *look=(vorbis_look_floor1 *)in; vorbis_info_floor1 *info=look->vi; long n=look->n; - /* first, quantize the floor into the units we're using, - silmultaneously weeding out values for x positions where the - residue will likely quantize to zero */ - int workfloor=alloca(sizeof(int)*n); - memset(workfloor,0,sizeof(int)*n); - for(i=0;i<n;i++){ - if + long posts=look->posts; + long nonzero=0; + lsfit_acc fits[VIF_POSIT+1]; + int fit_valueA[VIF_POSIT+2]; /* index by range list position */ + int fit_valueB[VIF_POSIT+2]; /* index by range list position */ + int fit_flag[VIF_POSIT+2]; + + int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */ + int hineighbor[VIF_POSIT+2]; + int memo[VIF_POSIT+2]; + static_codebook *sbooks=vb->vd->vi->codec_setup->book_param; + codebook *books=vb->vd->backend_state->fullbooks; + + memset(fit_flag,0,sizeof(fit_flag)); + for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */ + for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */ + for(i=0;i<posts;i++)memo[i]=-1; /* no neighbor yet */ + + /* quantize the relevant floor points and collect them into line fit + structures (one per minimal division) at the same time */ + if(posts==0){ + nonzero+=accumulate_fit(flr,0,n,fits); + }else{ + for(i=0;i<posts-1;i++) + nonzero+=accumulate_fit(flr,look->sorted_index[i], + look->sorted_index[i+1],fits+i); + } + + if(nonzero){ + /* start by fitting the implicit base case.... */ + int y0=-999; + int y1=-999; + int mse=fit_line(fits,posts-1,&y0,&y1); + if(mse<0){ + /* Only a single nonzero point */ + y0=-999; + y1=0; + fit_line(fits,posts-1,&y0,&y1); + } + + fit_flag[0]=1; + fit_flag[1]=1; + fit_valueA[0]=y0; + fit_valueB[0]=y0; + fit_valueB[1]=y1; + fit_valueA[1]=y1; + if(mse>=0){ + /* Non degenerate case */ + /* start progressive splitting. This is a greedy, non-optimal + algorithm, but simple and close enough to the best + answer. */ + for(i=2;i<posts;i++){ + int sortpos=look->reverse_index[i]; + int ln=loneighbor[sortpos]; + int hn=hineighbor[sortpos]; - return(val); + /* eliminate repeat searches of a particular range with a memo */ + if(memo[ln]!=hn){ + /* haven't performed this error search yet */ + int lsortpos=look->reverse_index[ln]; + int hsortpos=look->reverse_index[hn]; + + /* if this is an empty segment, its endpoints don't matter. + Mark as such */ + for(j=lsortpos;j<hsortpos;j++) + if(fits[j].n)break; + if(j==hsortpos){ + /* empty segment; important to note that this does not + break 0/n post case */ + fit_valueB[ln]=-1; + fit_valueA[hn]=-1; + }else{ + /* A note: we want to bound/minimize *local*, not global, error */ + int lx=info->postlist[ln]; + int hx=info->postlist[hn]; + int ly=post_Y(fit_valueA,fit_valueB,ln); + int hy=post_Y(fit_valueA,fit_valueB,hn); + memo[ln]=hn; + + if(i<info->searchstart || + inspect_error(lx,hx,ly,hy,flr,info)){ + /* outside error bounds/begin search area. Split it. */ + int ly0=-999; + int ly1=-999; + int hy0=-999; + int hy1=-999; + int lmse=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1); + int hmse=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1); + + /* Handle degeneracy */ + if(lmse<0 && hmse<0){ + ly0=fit_valueA[ln]; + hy1=fit_valueB[hn]; + lmse=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1); + hmse=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1); + } + + if(lmse<0 && hmse<0) continue; + + if(lmse<0){ + ly0=fit_valueA[ln]; + ly1=hy0; + lmse=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1); + } + if(hmse<0){ + hy1=fit_valueB[hn]; + hy0=ly1; + hmse=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1); + } + + /* store new edge values */ + fit_valueB[ln]=ly0; + if(ln==0)fit_valueA[ln]=ly0; + fit_valueA[i]=ly1; + fit_valueB[i]=hy0; + fit_valueA[hn]=hy1; + if(hn==1)fit_valueB[hn]=hy1; + + /* store new neighbor values */ + for(j=sortpos-1;j>=0;j--) + if(hineighbor[j]==hn) + hineighbor[j]=i; + else + break; + for(j=sortpos+1;j<posts;j++) + if(loneighbor[j]==ln) + loneighbor[j]=i; + else + break; + + /* store flag (set) */ + fit_flag[i]=1; + + } + } + } + } + } + + /* quantize values to multiplier spec */ + switch(info->mult){ + case 1: /* 1024 -> 256 */ + for(i=0;i<posts;i++) + if(fit_flag[i]) + fit_valueA[i]=(post_Y(fit_valueA,fit_valueB,i)+2)>>2; + break; + case 2: /* 1024 -> 128 */ + for(i=0;i<posts;i++) + if(fit_flag[i]) + fit_valueA[i]=(post_Y(fit_valueA,fit_valueB,i)+4)>>3; + break; + case 3: /* 1024 -> 86 */ + for(i=0;i<posts;i++) + if(fit_flag[i]) + fit_valueA[i]=(post_Y(fit_valueA,fit_valueB,i)+6)/12; + break; + case 4: /* 1024 -> 64 */ + for(i=0;i<posts;i++) + if(fit_flag[i]) + fit_valueA[i]=(post_Y(fit_valueA,fit_valueB,i)+8)>>4; + break; + } + + /* find prediction values for each post and subtract them */ + /* work backward to avoid a copy; unwind the neighbor arrays */ + for(i=posts-1;i>1;i--){ + int sp=look->reverse_index[i]; + int ln=loneighbor[i]; + int hn=hineighbor[i]; + + int x0=info->postlist[ln]; + int x1=info->postlist[hn]; + int y0=fit_valueA[ln]; + int y1=fit_valueA[hn]; + + int predicted=render_point(x0,x1,y0,y1,info->postlist[i]); + + if(fit_flag[i]){ + int headroom=(look->quant_q-predicted<predicted? + look->quant_q-predicted:predicted); + + int val=fit_valueA[i]-predicted; + + /* at this point the 'deviation' value is in the range +/- max + range, but the real, unique range can be mapped to only + [0-maxrange). So we want to wrap the deviation into this + limited range, but do it in the way that least screws an + essentially gaussian probability distribution. */ + + if(val<0) + if(val<-headroom) + val=headroom-val-1; + else + val=1-(val<<1); + else + if(val>=headroom) + val= val+headroom; + else + val>>=1; + + fit_valueB[i]=val; + + /* unroll the neighbor arrays */ + for(j=sp+1;j<posts;j++) + if(loneighbor[j]==i) + loneighbor[j]=loneighbor[sp]; + else + break; + for(j=sp-1;j>=0;j--) + if(hineighbor[j]==i) + hineighbor[j]=hineighbor[sp]; + else + break; + + }else{ + fit_valueA[i]=predicted; + fit_valueB[i]=0; + } + } + + /* we have everything we need. pack it out */ + /* mark nontrivial floor */ + oggpack_write(&vb->opb,1,1); + + /* beginning/end post */ + oggpack_write(&vb->opb,fit_valueA[0],ilog(look->quant_q-1)); + oggpack_write(&vb->opb,fit_valueA[1],ilog(look->quant_q-1)); + + /* partition by partition */ + for(i=0,j=2;i<info->partitions;i++){ + int class=info->partitionclass[i]; + int cdim=info->class_dim[class]; + int csubbits=info->class_subs[class]; + int csub=1<<csubbits; + int bookas[8]={0,0,0,0,0,0,0,0}; + int cval=0; + int cshift=0; + + /* generate the partition's first stage cascade value */ + if(csubbits){ + int maxval[8]; + for(k=0;k<csub;k++){ + int booknum=info->class_subbook[class][k]; + if(booknum<0){ + maxval[k]=0; + }else{ + maxval[k]=sbooks[info->class_subbook[class][k]].entries; + } + } + for(k=0;k<cdim;k++){ + for(l=0;l<csub;l++){ + val=fit_valueB[j+k]; + if(val<maxval[k]){ + bookas[k]=l; + break; + } + } + cval|= bookas[k]<<cshift; + cshift+=csubbits; + } + /* write it */ + vorbis_book_encode(books+info->class_book[class],cval,&vb->opb); + +#ifdef TRAIN_FLOOR1 + { + FILE *of; + char buffer[80]; + sprintf(buffer,"floor1_class%d.vqd",class); + of=fopen(buffer,"a"); + fprintf(of,"%d\n",cval); + fclose(of); + } +#endif + } + + /* write post values */ + for(k=0;k<cdim;k++){ + int book=info->class_subbook[class][bookas[k]]; + if(book>=0){ + vorbis_book_encode(books+book, + fit_valueB[j+k],&vb->opb); +#ifdef TRAIN_FLOOR1 + { + FILE *of; + char buffer[80]; + sprintf(buffer,"floor1_class%dsub%d.vqd",class,bookas[k]); + of=fopen(buffer,"a"); + fprintf(of,"%d\n",fit_valueB[j+k]); + fclose(of); + } +#endif + } + } + j+=cdim; + } + + /* generate quantized floor equivalent to what we'd unpack in decode */ + { + int hx; + int lx=0; + int ly=fit_valueA[0]*info->mult; + for(j=1;j<posts;j++){ + int current=look->forward_index[j]; + int hy=fit_valueA[current]*info->mult; + hx=info->postlist[current]; + + render_line(lx,hx,ly,hy,out); + + lx=hx; + ly=hy; + } + for(j=hx;j<n;j++)out[j]=0.f; /* be certain */ + } + + }else{ + oggpack_write(&vb->opb,0,1); + memset(flr,0,n*sizeof(float)); + } + return(nonzero); } static int floor1_inverse(vorbis_block *vb,vorbis_look_floor *i,float *out){ vorbis_look_floor1 *look=(vorbis_look_floor1 *)i; vorbis_info_floor1 *info=look->vi; + codec_setup_info *ci=vb->vd->vi->codec_setup; + int n=ci->blocksizes[vb->mode]/2; + int i,j; + codebook *books=vb->vd->backend_state->fullbooks; + + /* unpack wrapped/predicted values from stream */ + if(oggpack_read(&vb->opb,1)==1){ + int fit_value[VIF_POSIT+2]; + + fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1)); + fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1)); + + /* partition by partition */ + /* partition by partition */ + for(i=0,j=2;i<info->partitions;i++){ + int class=info->partitionclass[i]; + int cdim=info->class_dim[class]; + int csubbits=info->class_subs[class]; + int csub=1<<csubbits; + int bookas[8]={0,0,0,0,0,0,0,0}; + + /* decode the partition's first stage cascade value */ + if(csubbits){ + int cval=vorbis_book_decode(books+info->class_book[class],&vb->opb); + if(cval==-1)goto eop; + + for(k=0;k<cdim;k++){ + int book=info->class_subbook[class][val&(csub-1)]; + cval>>=csubbits; + if(book>=0){ + if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1) + goto eop; + }else{ + fit_value[j+k]=0; + } + } + }else{ + for(k=0;k<cdim;k++) + fit_value[j+k]=0; + } + j+=cdim; + } + /* unwrap positive values and reconsitute via linear interpolation */ + for(i=2;i<posts;i++){ + int predicted=render_point(info->postlist[look->loneighbor[i-2]], + info->postlist[look->hineighbor[i-2]], + fit_value[look->loneighbor[i-2]], + fit_value[look->hineighbor[i-2]], + info->postlist[i]); + int hiroom=look->quant_q-predicted; + int loroom=predicted; + int room=(hiroom<loroom?hiroom:loroom)<<1; + int val=fit_value[i]; + + if(val>=room){ + if(hiroom>loroom){ + val = val-loroom; + }else{ + val = -1-(val-hiroom); + } + }else{ + if(val&1){ + val= -1 -(val>>1); + }else{ + val>>=1; + } + } + + fit_value[i]=val+predicted; + } + + /* render the lines */ + { + int hx; + int lx=0; + int ly=fit_value[0]*info->mult; + for(j=1;j<posts;j++){ + int current=look->forward_index[j]; + int hy=fit_value[current]*info->mult; + hx=info->postlist[current]; + + render_line(lx,hx,ly,hy,out); + + lx=hx; + ly=hy; + } + for(j=hx;j<n;j++)out[j]=0.f; /* be certain */ + } + return(1); + } + + /* fall through */ eop: - memset(out,0,sizeof(float)*look->n); + memset(out,0,sizeof(float)*n); return(0); } @@ -285,4 +1025,3 @@ vorbis_func_floor floor1_exportbundle={ &floor1_free_look,&floor1_forward,&floor1_inverse }; - diff --git a/lib/mapping0.c b/lib/mapping0.c new file mode 100644 index 00000000..0e02d9f1 --- /dev/null +++ b/lib/mapping0.c @@ -0,0 +1,401 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * + * by the XIPHOPHORUS Company http://www.xiph.org/ * + + ******************************************************************** + + function: channel mapping 0 implementation + last mod: $Id: mapping0.c,v 1.27.4.1 2001/04/29 22:21:04 xiphmont Exp $ + + ********************************************************************/ + +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <math.h> +#include <ogg/ogg.h> +#include "vorbis/codec.h" +#include "codec_internal.h" +#include "codebook.h" +#include "bitbuffer.h" +#include "registry.h" +#include "psy.h" +#include "misc.h" + +/* simplistic, wasteful way of doing this (unique lookup for each + mode/submapping); there should be a central repository for + identical lookups. That will require minor work, so I'm putting it + off as low priority. + + Why a lookup for each backend in a given mode? Because the + blocksize is set by the mode, and low backend lookups may require + parameters from other areas of the mode/mapping */ + +typedef struct { + drft_lookup fft_look; + vorbis_info_mode *mode; + vorbis_info_mapping0 *map; + + vorbis_look_time **time_look; + vorbis_look_floor **floor_look; + + vorbis_look_residue **residue_look; + vorbis_look_psy *psy_look; + + vorbis_func_time **time_func; + vorbis_func_floor **floor_func; + vorbis_func_residue **residue_func; + + int ch; + long lastframe; /* if a different mode is called, we need to + invalidate decay */ +} vorbis_look_mapping0; + +static vorbis_info_mapping *mapping0_copy_info(vorbis_info_mapping *vm){ + vorbis_info_mapping0 *info=(vorbis_info_mapping0 *)vm; + vorbis_info_mapping0 *ret=_ogg_malloc(sizeof(vorbis_info_mapping0)); + memcpy(ret,info,sizeof(vorbis_info_mapping0)); + return(ret); +} + +static void mapping0_free_info(vorbis_info_mapping *i){ + if(i){ + memset(i,0,sizeof(vorbis_info_mapping0)); + _ogg_free(i); + } +} + +static void mapping0_free_look(vorbis_look_mapping *look){ + int i; + vorbis_look_mapping0 *l=(vorbis_look_mapping0 *)look; + if(l){ + drft_clear(&l->fft_look); + + for(i=0;i<l->map->submaps;i++){ + l->time_func[i]->free_look(l->time_look[i]); + l->floor_func[i]->free_look(l->floor_look[i]); + l->residue_func[i]->free_look(l->residue_look[i]); + if(l->psy_look)_vp_psy_clear(l->psy_look+i); + } + + _ogg_free(l->time_func); + _ogg_free(l->floor_func); + _ogg_free(l->residue_func); + _ogg_free(l->time_look); + _ogg_free(l->floor_look); + _ogg_free(l->residue_look); + if(l->psy_look)_ogg_free(l->psy_look); + memset(l,0,sizeof(vorbis_look_mapping0)); + _ogg_free(l); + } +} + +static vorbis_look_mapping *mapping0_look(vorbis_dsp_state *vd,vorbis_info_mode *vm, + vorbis_info_mapping *m){ + int i; + vorbis_info *vi=vd->vi; + codec_setup_info *ci=vi->codec_setup; + vorbis_look_mapping0 *look=_ogg_calloc(1,sizeof(vorbis_look_mapping0)); + vorbis_info_mapping0 *info=look->map=(vorbis_info_mapping0 *)m; + look->mode=vm; + + look->time_look=_ogg_calloc(info->submaps,sizeof(vorbis_look_time *)); + look->floor_look=_ogg_calloc(info->submaps,sizeof(vorbis_look_floor *)); + + look->residue_look=_ogg_calloc(info->submaps,sizeof(vorbis_look_residue *)); + if(ci->psys)look->psy_look=_ogg_calloc(info->submaps,sizeof(vorbis_look_psy)); + + look->time_func=_ogg_calloc(info->submaps,sizeof(vorbis_func_time *)); + look->floor_func=_ogg_calloc(info->submaps,sizeof(vorbis_func_floor *)); + look->residue_func=_ogg_calloc(info->submaps,sizeof(vorbis_func_residue *)); + + for(i=0;i<info->submaps;i++){ + int timenum=info->timesubmap[i]; + int floornum=info->floorsubmap[i]; + int resnum=info->residuesubmap[i]; + + look->time_func[i]=_time_P[ci->time_type[timenum]]; + look->time_look[i]=look->time_func[i]-> + look(vd,vm,ci->time_param[timenum]); + look->floor_func[i]=_floor_P[ci->floor_type[floornum]]; + look->floor_look[i]=look->floor_func[i]-> + look(vd,vm,ci->floor_param[floornum]); + look->residue_func[i]=_residue_P[ci->residue_type[resnum]]; + look->residue_look[i]=look->residue_func[i]-> + look(vd,vm,ci->residue_param[resnum]); + + if(ci->psys && vd->analysisp){ + int psynum=info->psysubmap[i]; + _vp_psy_init(look->psy_look+i,ci->psy_param[psynum], + ci->blocksizes[vm->blockflag]/2,vi->rate); + } + } + + look->ch=vi->channels; + + if(vd->analysisp)drft_init(&look->fft_look,ci->blocksizes[vm->blockflag]); + return(look); +} + +static void mapping0_pack(vorbis_info *vi,vorbis_info_mapping *vm,oggpack_buffer *opb){ + int i; + vorbis_info_mapping0 *info=(vorbis_info_mapping0 *)vm; + + oggpack_write(opb,info->submaps-1,4); + /* we don't write the channel submappings if we only have one... */ + if(info->submaps>1){ + for(i=0;i<vi->channels;i++) + oggpack_write(opb,info->chmuxlist[i],4); + } + for(i=0;i<info->submaps;i++){ + oggpack_write(opb,info->timesubmap[i],8); + oggpack_write(opb,info->floorsubmap[i],8); + oggpack_write(opb,info->residuesubmap[i],8); + } +} + +/* also responsible for range checking */ +static vorbis_info_mapping *mapping0_unpack(vorbis_info *vi,oggpack_buffer *opb){ + int i; + vorbis_info_mapping0 *info=_ogg_calloc(1,sizeof(vorbis_info_mapping0)); + codec_setup_info *ci=vi->codec_setup; + memset(info,0,sizeof(vorbis_info_mapping0)); + + info->submaps=oggpack_read(opb,4)+1; + + if(info->submaps>1){ + for(i=0;i<vi->channels;i++){ + info->chmuxlist[i]=oggpack_read(opb,4); + if(info->chmuxlist[i]>=info->submaps)goto err_out; + } + } + for(i=0;i<info->submaps;i++){ + info->timesubmap[i]=oggpack_read(opb,8); + if(info->timesubmap[i]>=ci->times)goto err_out; + info->floorsubmap[i]=oggpack_read(opb,8); + if(info->floorsubmap[i]>=ci->floors)goto err_out; + info->residuesubmap[i]=oggpack_read(opb,8); + if(info->residuesubmap[i]>=ci->residues)goto err_out; + } + + return info; + + err_out: + mapping0_free_info(info); + return(NULL); +} + +#include "os.h" +#include "lpc.h" +#include "lsp.h" +#include "envelope.h" +#include "mdct.h" +#include "psy.h" +#include "scales.h" + +/* no time mapping implementation for now */ +static long seq=0; +static int mapping0_forward(vorbis_block *vb,vorbis_look_mapping *l){ + vorbis_dsp_state *vd=vb->vd; + vorbis_info *vi=vd->vi; + backend_lookup_state *b=vb->vd->backend_state; + vorbis_look_mapping0 *look=(vorbis_look_mapping0 *)l; + vorbis_info_mapping0 *info=look->map; + vorbis_info_mode *mode=look->mode; + vorbis_block_internal *vbi=(vorbis_block_internal *)vb->internal; + int n=vb->pcmend; + int i,j; + float *window=b->window[vb->W][vb->lW][vb->nW][mode->windowtype]; + + float **pcmbundle=alloca(sizeof(float *)*vi->channels); + + int *nonzero=alloca(sizeof(int)*vi->channels); + + float **floor=_vorbis_block_alloc(vb,vi->channels*sizeof(float *)); + float *additional=_vorbis_block_alloc(vb,n*sizeof(float)); + float newmax=vbi->ampmax; + + for(i=0;i<vi->channels;i++){ + float *pcm=vb->pcm[i]; + float scale=4.f/n; + int submap=info->chmuxlist[i]; + float ret; + + _analysis_output("pcm",seq,pcm,n,0,0); + + /* window the PCM data */ + for(j=0;j<n;j++) + additional[j]=pcm[j]*=window[j]; + + _analysis_output("windowed",seq,pcm,n,0,0); + + /* transform the PCM data */ + /* only MDCT right now.... */ + mdct_forward(b->transform[vb->W][0],pcm,pcm); + + /* FFT yields more accurate tonal estimation (not phase sensitive) */ + drft_forward(&look->fft_look,additional); + + additional[0]=fabs(additional[0]*scale); + for(j=1;j<n-1;j+=2) + additional[(j+1)>>1]=scale*FAST_HYPOT(additional[j],additional[j+1]); + + /* set up our masking data working vector, and stash unquantized + data for later */ + /*memcpy(pcm+n/2,pcm,n*sizeof(float)/2);*/ + memcpy(additional+n/2,pcm,n*sizeof(float)/2); + + /* begin masking work */ + floor[i]=_vorbis_block_alloc(vb,n*sizeof(float)/2); + + _analysis_output("fft",seq,additional,n/2,0,1); + _analysis_output("mdct",seq,additional+n/2,n/2,0,1); + _analysis_output("lfft",seq,additional,n/2,0,0); + _analysis_output("lmdct",seq,additional+n/2,n/2,0,0); + + /* perform psychoacoustics; do masking */ + ret=_vp_compute_mask(look->psy_look+submap,additional,additional+n/2, + floor[i],NULL,vbi->ampmax); + if(ret>newmax)newmax=ret; + + _analysis_output_always("mask",seq,floor[i],n/2,0,0); + + /* perform floor encoding */ + nonzero[i]=look->floor_func[submap]-> + forward(vb,look->floor_look[submap],floor[i]); + + _analysis_output("floor",seq,floor[i],n/2,0,1); + + /* apply the floor, do optional noise levelling */ + _vp_apply_floor(look->psy_look+submap,pcm,floor[i]); + + _analysis_output("res",seq++,pcm,n/2,0,0); + +#ifdef TRAIN_RES + if(nonzero[i]){ + FILE *of; + char buffer[80]; + int i; + + sprintf(buffer,"residue_%d.vqd",vb->mode); + of=fopen(buffer,"a"); + for(i=0;i<n/2;i++){ + fprintf(of,"%.2f, ",pcm[i]); + if(fabs(pcm[i])>1000){ + fprintf(stderr," %d ",seq-1); + } + } + fprintf(of,"\n"); + fclose(of); + } +#endif + } + + seq++; + vbi->ampmax=newmax; + + /* perform residue encoding with residue mapping; this is + multiplexed. All the channels belonging to one submap are + encoded (values interleaved), then the next submap, etc */ + + for(i=0;i<info->submaps;i++){ + int ch_in_bundle=0; + for(j=0;j<vi->channels;j++){ + if(info->chmuxlist[j]==i && nonzero[j]>0){ + pcmbundle[ch_in_bundle++]=vb->pcm[j]; + } + } + + look->residue_func[i]->forward(vb,look->residue_look[i], + pcmbundle,ch_in_bundle); + } + + look->lastframe=vb->sequence; + return(0); +} + +static int mapping0_inverse(vorbis_block *vb,vorbis_look_mapping *l){ + vorbis_dsp_state *vd=vb->vd; + vorbis_info *vi=vd->vi; + codec_setup_info *ci=vi->codec_setup; + backend_lookup_state *b=vd->backend_state; + vorbis_look_mapping0 *look=(vorbis_look_mapping0 *)l; + vorbis_info_mapping0 *info=look->map; + vorbis_info_mode *mode=look->mode; + int i,j; + long n=vb->pcmend=ci->blocksizes[vb->W]; + + float *window=b->window[vb->W][vb->lW][vb->nW][mode->windowtype]; + float **pcmbundle=alloca(sizeof(float *)*vi->channels); + int *nonzero=alloca(sizeof(int)*vi->channels); + + /* time domain information decode (note that applying the + information would have to happen later; we'll probably add a + function entry to the harness for that later */ + /* NOT IMPLEMENTED */ + + /* recover the spectral envelope; store it in the PCM vector for now */ + for(i=0;i<vi->channels;i++){ + float *pcm=vb->pcm[i]; + int submap=info->chmuxlist[i]; + nonzero[i]=look->floor_func[submap]-> + inverse(vb,look->floor_look[submap],pcm); + _analysis_output("ifloor",seq+i,pcm,n/2,0,1); + } + + /* recover the residue, apply directly to the spectral envelope */ + + for(i=0;i<info->submaps;i++){ + int ch_in_bundle=0; + for(j=0;j<vi->channels;j++){ + if(info->chmuxlist[j]==i && nonzero[j]) + pcmbundle[ch_in_bundle++]=vb->pcm[j]; + } + + look->residue_func[i]->inverse(vb,look->residue_look[i],pcmbundle,ch_in_bundle); + } + + /* transform the PCM data; takes PCM vector, vb; modifies PCM vector */ + /* only MDCT right now.... */ + for(i=0;i<vi->channels;i++){ + float *pcm=vb->pcm[i]; + _analysis_output("out",seq+i,pcm,n/2,0,1); + mdct_backward(b->transform[vb->W][0],pcm,pcm); + } + + /* now apply the decoded pre-window time information */ + /* NOT IMPLEMENTED */ + + /* window the data */ + for(i=0;i<vi->channels;i++){ + float *pcm=vb->pcm[i]; + if(nonzero[i]) + for(j=0;j<n;j++) + pcm[j]*=window[j]; + else + for(j=0;j<n;j++) + pcm[j]=0.f; + _analysis_output("final",seq++,pcm,n,0,0); + } + + /* now apply the decoded post-window time information */ + /* NOT IMPLEMENTED */ + + /* all done! */ + return(0); +} + +/* export hooks */ +vorbis_func_mapping mapping0_exportbundle={ + &mapping0_pack,&mapping0_unpack,&mapping0_look,&mapping0_copy_info, + &mapping0_free_info,&mapping0_free_look,&mapping0_forward,&mapping0_inverse +}; + + + diff --git a/lib/modes/mode_A.h b/lib/modes/mode_A.h new file mode 100644 index 00000000..09694c04 --- /dev/null +++ b/lib/modes/mode_A.h @@ -0,0 +1,316 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * + * by the XIPHOPHORUS Company http://www.xiph.org/ * + + ******************************************************************** + + function: predefined encoding modes + last mod: $Id: mode_A.h,v 1.14.4.1 2001/04/29 22:21:06 xiphmont Exp $ + + ********************************************************************/ + +#ifndef _V_MODES_A_H_ +#define _V_MODES_A_H_ + +#include <stdio.h> +#include "vorbis/codec.h" +#include "backends.h" + +#include "books/lsp12_0.vqh" +#include "books/lsp30_0.vqh" +#include "books/lsp12_1.vqh" +#include "books/lsp30_1.vqh" + +#include "books/res0_128_128aux.vqh" +#include "books/res0_128_1024aux.vqh" + +#include "books/res0_128_128_1.vqh" +#include "books/res0_128_128_2.vqh" +#include "books/res0_128_128_3.vqh" +#include "books/res0_128_128_4.vqh" +#include "books/res0_128_128_5.vqh" + +#include "books/res0_128_1024_1.vqh" +#include "books/res0_128_1024_2.vqh" +#include "books/res0_128_1024_3.vqh" +#include "books/res0_128_1024_4.vqh" +#include "books/res0_128_1024_5.vqh" +#include "books/res0_128_1024_6.vqh" +#include "books/res0_128_1024_7.vqh" +#include "books/res0_128_1024_8.vqh" +#include "books/res0_128_1024_9.vqh" + + +static vorbis_info_psy _psy_set_A0={ + 1,/*athp*/ + 1,/*decayp*/ + + -100., + -140., + + 8, + + /* 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 */ + /* x: 63 88 125 175 250 350 500 700 1k 1.4k 2k 2.8k 4k 5.6k 8k 11.5k 16k Hz */ + /* y: 0 10 20 30 40 50 60 70 80 90 100 dB */ + + 1,/* tonemaskp */ + /* 0 10 20 30 40 50 60 70 80 90 100 */ + { + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*63*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*88*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*125*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*175*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*250*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*350*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*500*/ + {-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}, /*700*/ + + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1000*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1400*/ + {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2000*/ + {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2800*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*4000*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*5600*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*8000*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*11500*/ + {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*16000*/ + }, + + 1,/* peakattp */ + {{-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-22.f,-22.f,-22.f},/*63*/ + {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-22.f,-22.f,-22.f},/*88*/ + {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-22.f,-22.f,-22.f},/*125*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-22.f,-22.f},/*175*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-22.f,-22.f},/*250*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-22.f,-22.f},/*350*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-22.f,-22.f},/*500*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-14.f,-20.f,-22.f,-22.f,-22.f},/*700*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-14.f,-20.f,-22.f,-22.f,-22.f},/*1000*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-22.f,-22.f},/*1400*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-22.f,-22.f},/*2000*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-22.f,-22.f},/*2400*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-22.f,-22.f},/*4000*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-11.f,-12.f,-13.f,-22.f,-22.f,-22.f},/*5600*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-11.f,-12.f,-13.f,-22.f,-22.f,-22.f},/*8000*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-11.f,-22.f,-22.f,-22.f},/*11500*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-20.f,-21.f,-22.f},/*16000*/ + }, + + 1,/*noisemaskp */ + -10.f, /* suppress any noise curve over maxspec+n */ + .5f, /* low window */ + .5f, /* high window */ + 5, + 5, + {.000f, 0.f,/*63*/ + .000f, 0.f,/*88*/ + .000f, 0.f,/*125*/ + .000f, 0.f,/*175*/ + .000f, 0.f,/*250*/ + .000f, 0.f,/*350*/ + .000f, 0.f,/*500*/ + .300f, 0.f,/*700*/ + .500f, 0.f,/*1000*/ + .500f, 0.f,/*1400*/ + .500f, 0.f,/*2000*/ + .500f, 0.f,/*2800*/ + .600f, 0.f,/*4000*/ + .700f, 0.f,/*5600*/ + .850f, 0.f,/*8000*/ + .900f, 0.f,/*11500*/ + .900f, 0.f,/*16000*/ + }, + + 95.f, /* even decade + 5 is important; saves an rint() later in a + tight loop) */ + -22., + + 1, +}; + +static vorbis_info_psy _psy_set_A={ + 1,/*athp*/ + 1,/*decayp*/ + + -100.f, + -140.f, + + 8, + + /* 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 */ + /* x: 63 88 125 175 250 350 500 700 1k 1.4k 2k 2.8k 4k 5.6k 8k 11.5k 16k Hz */ + /* y: 0 10 20 30 40 50 60 70 80 90 100 dB */ + 1,/* tonemaskp */ + /* 0 10 20 30 40 50 60 70 80 90 100 */ + { + {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*63*/ + {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*88*/ + {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*125*/ + + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*175*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*250*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*350*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*500*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*700*/ + + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1000*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1400*/ + {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2000*/ + {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2800*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*4000*/ + {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*5600*/ + + {-30.f,-30.f,-33.f,-35.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*8000*/ + {-30.f,-30.f,-33.f,-35.f,-35.f,-45.f,-50.f,-60.f,-70.f,-90.f,-100.f}, /*11500*/ + {-24.f,-24.f,-26.f,-32.f,-32.f,-42.f,-50.f,-60.f,-70.f,-90.f,-100.f}, /*16000*/ + + }, + + 1,/* peakattp */ + {{-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*63*/ + {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*88*/ + {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*125*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*175*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*250*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*350*/ + {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*500*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*700*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1000*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1400*/ + {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2000*/ + {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2400*/ + {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*4000*/ + {-10.f,-10.f,-10.f,-12.f,-12.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*5600*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*8000*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*11500*/ + {-10.f,-10.f,-10.f,-10.f,-10.f,-12.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*16000*/ + }, + + 1,/*noisemaskp */ + -24.f, /* suppress any noise curve over maxspec+n */ + .5f, /* low window */ + .5f, /* high window */ + 25, + 25, + {.000f, 0.f, /*63*/ + .000f, 0.f, /*88*/ + .000f, 0.f, /*125*/ + .000f, 0.f, /*175*/ + .000f, 0.f, /*250*/ + .000f, 0.f, /*350*/ + .000f, 0.f, /*500*/ + .200f, 0.f, /*700*/ + .300f, 0.f, /*1000*/ + .400f, 0.f, /*1400*/ + .400f, 0.f, /*2000*/ + .400f, 0.f, /*2800*/ + .700f, 0.f, /*4000*/ + .850f, 0.f, /*5600*/ + .900f, 0.f, /*8000*/ + .900f, 0.f, /*11500*/ + .900f, 1.f, /*16000*/ + }, + + 95.f, /* even decade + 5 is important; saves an rint() later in a + tight loop) */ + -28., + + 1, + +}; + +/* with GNUisms, this could be short and readable. Oh well */ +static vorbis_info_time0 _time_set0A={0}; +static vorbis_info_floor0 _floor_set0A={12, 44100, 64, 10,130, 2, {0,1}, + 0.199f, .285f}; +static vorbis_info_floor0 _floor_set1A={30, 44100, 256, 12,150, 2, {2,3}, + .082f, .126f}; +static vorbis_info_residue0 _residue_set0A={0,96,16,6,4, + {0,1,1,1,1,1}, + {6,7,8,9,10}, + + {0,99999,9999,9999,9999}, + {999.f,1.5f,3.5f,15.5f,26.5f}, + {4,4,4,4,4}, + {99,99,99,99,99}}; + +static vorbis_info_residue0 _residue_set1A={0, 960, 32,10,5, + {0,1,1,1,1,1,1,1,1,1}, + {11,12,13,14,15,16,17,18,19}, + + {0,8,9999,16,9999, + 24,9999,9999,9999}, + {999.f,1.5f,1.5f,2.5f,2.5f, + 6.5f,6.5f,12.5f,22.5f}, + {5,5,5,5,5,5,5,5,5}, + {99,99,99,99,99,99,99,99,99}}; + +static vorbis_info_mapping0 _mapping_set0A={1, {0,0}, {0}, {0}, {0}, {0}}; +static vorbis_info_mapping0 _mapping_set1A={1, {0,0}, {0}, {1}, {1}, {1}}; +static vorbis_info_mode _mode_set0A={0,0,0,0}; +static vorbis_info_mode _mode_set1A={1,0,0,1}; + +/* CD quality stereo, no channel coupling */ +codec_setup_info info_A={ + + /* smallblock, largeblock */ + {256, 2048}, + /* modes,maps,times,floors,residues,books,psys */ + 2, 2, 1, 2, 2, 20, 2, + /* modes */ + {&_mode_set0A,&_mode_set1A}, + /* maps */ + {0,0},{&_mapping_set0A,&_mapping_set1A}, + /* times */ + {0,0},{&_time_set0A}, + /* floors */ + {0,0},{&_floor_set0A,&_floor_set1A}, + /* residue */ + {0,0},{&_residue_set0A,&_residue_set1A}, + /* books */ + {&_vq_book_lsp12_0, /* 0 */ + &_vq_book_lsp12_1, /* 1 */ + &_vq_book_lsp30_0, /* 2 */ + &_vq_book_lsp30_1, /* 3 */ + + &_huff_book_res0_128_128aux, + &_huff_book_res0_128_1024aux, + + &_vq_book_res0_128_128_1, + &_vq_book_res0_128_128_2, + &_vq_book_res0_128_128_3, + &_vq_book_res0_128_128_4, + &_vq_book_res0_128_128_5, + + &_vq_book_res0_128_1024_1, + &_vq_book_res0_128_1024_2, + &_vq_book_res0_128_1024_3, + &_vq_book_res0_128_1024_4, + &_vq_book_res0_128_1024_5, + &_vq_book_res0_128_1024_6, + &_vq_book_res0_128_1024_7, + &_vq_book_res0_128_1024_8, + &_vq_book_res0_128_1024_9, + + }, + /* psy */ + {&_psy_set_A0,&_psy_set_A}, + + /* thresh sample period, preecho clamp trigger threshhold, range, minenergy */ + 256, {26.f,26.f,26.f,30.f}, {-90.f,-90.f,-90.f,-90.f}, -90.f, + + -10., + + 0, +}; + +#define PREDEF_INFO_MAX 0 + +#endif diff --git a/lib/psy.c b/lib/psy.c new file mode 100644 index 00000000..8cd4272e --- /dev/null +++ b/lib/psy.c @@ -0,0 +1,727 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * + * by the XIPHOPHORUS Company http://www.xiph.org/ * + + ******************************************************************** + + function: psychoacoustics not including preecho + last mod: $Id: psy.c,v 1.44.4.1 2001/04/29 22:21:04 xiphmont Exp $ + + ********************************************************************/ + +#include <stdlib.h> +#include <math.h> +#include <string.h> +#include "vorbis/codec.h" +#include "codec_internal.h" + +#include "masking.h" +#include "psy.h" +#include "os.h" +#include "lpc.h" +#include "smallft.h" +#include "scales.h" +#include "misc.h" + +#define NEGINF -9999.f + +/* Why Bark scale for encoding but not masking computation? Because + masking has a strong harmonic dependancy */ + +/* the beginnings of real psychoacoustic infrastructure. This is + still not tightly tuned */ +void _vi_psy_free(vorbis_info_psy *i){ + if(i){ + memset(i,0,sizeof(vorbis_info_psy)); + _ogg_free(i); + } +} + +vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i){ + vorbis_info_psy *ret=_ogg_malloc(sizeof(vorbis_info_psy)); + memcpy(ret,i,sizeof(vorbis_info_psy)); + return(ret); +} + +/* Set up decibel threshold slopes on a Bark frequency scale */ +/* ATH is the only bit left on a Bark scale. No reason to change it + right now */ +static void set_curve(float *ref,float *c,int n, float crate){ + int i,j=0; + + for(i=0;i<MAX_BARK-1;i++){ + int endpos=rint(fromBARK(i+1)*2*n/crate); + float base=ref[i]; + if(j<endpos){ + float delta=(ref[i+1]-base)/(endpos-j); + for(;j<endpos && j<n;j++){ + c[j]=base; + base+=delta; + } + } + } +} + +static void min_curve(float *c, + float *c2){ + int i; + for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i]; +} +static void max_curve(float *c, + float *c2){ + int i; + for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i]; +} + +static void attenuate_curve(float *c,float att){ + int i; + for(i=0;i<EHMER_MAX;i++) + c[i]+=att; +} + +static void interp_curve(float *c,float *c1,float *c2,float del){ + int i; + for(i=0;i<EHMER_MAX;i++) + c[i]=c2[i]*del+c1[i]*(1.f-del); +} + +static void setup_curve(float **c, + int band, + float *curveatt_dB){ + int i,j; + float ath[EHMER_MAX]; + float tempc[P_LEVELS][EHMER_MAX]; + + memcpy(c[0]+2,c[4]+2,sizeof(float)*EHMER_MAX); + memcpy(c[2]+2,c[4]+2,sizeof(float)*EHMER_MAX); + + /* we add back in the ATH to avoid low level curves falling off to + -infinity and unneccessarily cutting off high level curves in the + curve limiting (last step). But again, remember... a half-band's + settings must be valid over the whole band, and it's better to + mask too little than too much, so be pessimal. */ + + for(i=0;i<EHMER_MAX;i++){ + float oc_min=band*.5+(i-EHMER_OFFSET)*.125; + float oc_max=band*.5+(i-EHMER_OFFSET+1)*.125; + float bark=toBARK(fromOC(oc_min)); + int ibark=floor(bark); + float del=bark-ibark; + float ath_min,ath_max; + + if(ibark<26) + ath_min=ATH_Bark_dB[ibark]*(1.f-del)+ATH_Bark_dB[ibark+1]*del; + else + ath_min=ATH_Bark_dB[25]; + + bark=toBARK(fromOC(oc_max)); + ibark=floor(bark); + del=bark-ibark; + + if(ibark<26) + ath_max=ATH_Bark_dB[ibark]*(1.f-del)+ATH_Bark_dB[ibark+1]*del; + else + ath_max=ATH_Bark_dB[25]; + + ath[i]=min(ath_min,ath_max); + } + + /* The c array is comes in as dB curves at 20 40 60 80 100 dB. + interpolate intermediate dB curves */ + for(i=1;i<P_LEVELS;i+=2){ + interp_curve(c[i]+2,c[i-1]+2,c[i+1]+2,.5); + } + + /* normalize curves so the driving amplitude is 0dB */ + /* make temp curves with the ATH overlayed */ + for(i=0;i<P_LEVELS;i++){ + attenuate_curve(c[i]+2,curveatt_dB[i]); + memcpy(tempc[i],ath,EHMER_MAX*sizeof(float)); + attenuate_curve(tempc[i],-i*10.f); + max_curve(tempc[i],c[i]+2); + } + + /* Now limit the louder curves. + + the idea is this: We don't know what the playback attenuation + will be; 0dB SL moves every time the user twiddles the volume + knob. So that means we have to use a single 'most pessimal' curve + for all masking amplitudes, right? Wrong. The *loudest* sound + can be in (we assume) a range of ...+100dB] SL. However, sounds + 20dB down will be in a range ...+80], 40dB down is from ...+60], + etc... */ + + for(j=1;j<P_LEVELS;j++){ + min_curve(tempc[j],tempc[j-1]); + min_curve(c[j]+2,tempc[j]); + } + + /* add fenceposts */ + for(j=0;j<P_LEVELS;j++){ + + for(i=0;i<EHMER_MAX;i++) + if(c[j][i+2]>-200.f)break; + c[j][0]=i; + + for(i=EHMER_MAX-1;i>=0;i--) + if(c[j][i+2]>-200.f) + break; + c[j][1]=i; + + } +} + +void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){ + long i,j; + long maxoc; + memset(p,0,sizeof(vorbis_look_psy)); + + + p->eighth_octave_lines=vi->eighth_octave_lines; + p->shiftoc=rint(log(vi->eighth_octave_lines*8)/log(2))-1; + + p->firstoc=toOC(.25f*rate/n)*(1<<(p->shiftoc+1))-vi->eighth_octave_lines; + maxoc=toOC((n*.5f-.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f; + p->total_octave_lines=maxoc-p->firstoc+1; + + p->ath=_ogg_malloc(n*sizeof(float)); + p->octave=_ogg_malloc(n*sizeof(long)); + p->bark=_ogg_malloc(n*sizeof(float)); + p->vi=vi; + p->n=n; + + /* set up the lookups for a given blocksize and sample rate */ + /* Vorbis max sample rate is currently limited by 26 Bark (54kHz) */ + set_curve(ATH_Bark_dB, p->ath,n,rate); + for(i=0;i<n;i++) + p->bark[i]=toBARK(rate/(2*n)*i); + + for(i=0;i<n;i++) + p->octave[i]=toOC((i*.5f+.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f; + + p->tonecurves=_ogg_malloc(P_BANDS*sizeof(float **)); + p->noisemedian=_ogg_malloc(n*sizeof(float)); + p->noiseoffset=_ogg_malloc(n*sizeof(float)); + p->peakatt=_ogg_malloc(P_BANDS*sizeof(float *)); + for(i=0;i<P_BANDS;i++){ + p->tonecurves[i]=_ogg_malloc(P_LEVELS*sizeof(float *)); + p->peakatt[i]=_ogg_malloc(P_LEVELS*sizeof(float)); + } + + for(i=0;i<P_BANDS;i++) + for(j=0;j<P_LEVELS;j++){ + p->tonecurves[i][j]=_ogg_malloc((EHMER_MAX+2)*sizeof(float)); + } + + /* OK, yeah, this was a silly way to do it */ + memcpy(p->tonecurves[0][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[0][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[0][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[0][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[2][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[2][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[2][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[2][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[4][4]+2,tone_250_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[4][6]+2,tone_250_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[4][8]+2,tone_250_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[4][10]+2,tone_250_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[6][4]+2,tone_500_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[6][6]+2,tone_500_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[6][8]+2,tone_500_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[6][10]+2,tone_500_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[8][4]+2,tone_1000_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[8][6]+2,tone_1000_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[8][8]+2,tone_1000_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[8][10]+2,tone_1000_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[10][4]+2,tone_2000_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[10][6]+2,tone_2000_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[10][8]+2,tone_2000_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[10][10]+2,tone_2000_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[12][4]+2,tone_4000_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[12][6]+2,tone_4000_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[12][8]+2,tone_4000_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[12][10]+2,tone_4000_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[14][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[14][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[14][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[14][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX); + + memcpy(p->tonecurves[16][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[16][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[16][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX); + memcpy(p->tonecurves[16][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX); + + /* interpolate curves between */ + for(i=1;i<P_BANDS;i+=2) + for(j=4;j<P_LEVELS;j+=2){ + memcpy(p->tonecurves[i][j]+2,p->tonecurves[i-1][j]+2,EHMER_MAX*sizeof(float)); + /*interp_curve(p->tonecurves[i][j], + p->tonecurves[i-1][j], + p->tonecurves[i+1][j],.5);*/ + min_curve(p->tonecurves[i][j]+2,p->tonecurves[i+1][j]+2); + } + + /* set up the final curves */ + for(i=0;i<P_BANDS;i++) + setup_curve(p->tonecurves[i],i,vi->toneatt[i]); + + /* set up attenuation levels */ + for(i=0;i<P_BANDS;i++) + for(j=0;j<P_LEVELS;j++){ + p->peakatt[i][j]=p->vi->peakatt[i][j]; + } + + /* set up rolling noise median */ + for(i=0;i<n;i++){ + float halfoc=toOC((i+.5)*rate/(2.*n))*2.+2.; + int inthalfoc; + float del; + + if(halfoc<0)halfoc=0; + if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1; + inthalfoc=(int)halfoc; + del=halfoc-inthalfoc; + + p->noisemedian[i]= + p->vi->noisemedian[inthalfoc*2]*(1.-del) + + p->vi->noisemedian[inthalfoc*2+2]*del; + p->noiseoffset[i]= + p->vi->noisemedian[inthalfoc*2+1]*(1.-del) + + p->vi->noisemedian[inthalfoc*2+3]*del; + } + /*_analysis_output("mediancurve",0,p->noisemedian,n,0,0);*/ +} + +void _vp_psy_clear(vorbis_look_psy *p){ + int i,j; + if(p){ + if(p->ath)_ogg_free(p->ath); + if(p->octave)_ogg_free(p->octave); + if(p->bark)_ogg_free(p->bark); + if(p->tonecurves){ + for(i=0;i<P_BANDS;i++){ + for(j=0;j<P_LEVELS;j++){ + _ogg_free(p->tonecurves[i][j]); + } + _ogg_free(p->tonecurves[i]); + _ogg_free(p->peakatt[i]); + } + _ogg_free(p->tonecurves); + _ogg_free(p->noisemedian); + _ogg_free(p->noiseoffset); + _ogg_free(p->peakatt); + } + memset(p,0,sizeof(vorbis_look_psy)); + } +} + +/* octave/(8*eighth_octave_lines) x scale and dB y scale */ +static void seed_curve(float *seed, + float **curves, + float amp, + int oc,int n,int linesper,float dBoffset){ + int i; + long seedptr; + float *posts,*curve; + + int choice=(int)((amp+dBoffset)*.1f); + choice=max(choice,0); + choice=min(choice,P_LEVELS-1); + posts=curves[choice]; + curve=posts+2; + seedptr=oc+(posts[0]-16)*linesper-(linesper>>1); + + for(i=posts[0];i<posts[1];i++){ + if(seedptr>0){ + float lin=amp+curve[i]; + if(seed[seedptr]<lin)seed[seedptr]=lin; + } + seedptr+=linesper; + if(seedptr>=n)break; + } +} + +static void seed_peak(float *seed, + float *att, + float amp, + int oc, + int linesper, + float dBoffset){ + long seedptr; + + int choice=(int)((amp+dBoffset)*.1f); + choice=max(choice,0); + choice=min(choice,P_LEVELS-1); + seedptr=oc-(linesper>>1); + + amp+=att[choice]; + if(seed[seedptr]<amp)seed[seedptr]=amp; + +} + +static void seed_loop(vorbis_look_psy *p, + float ***curves, + float **att, + float *f, + float *flr, + float *minseed, + float *maxseed, + float specmax){ + vorbis_info_psy *vi=p->vi; + long n=p->n,i; + float dBoffset=vi->max_curve_dB-specmax; + + /* prime the working vector with peak values */ + + for(i=0;i<n;i++){ + float max=f[i]; + long oc=p->octave[i]; + while(i+1<n && p->octave[i+1]==oc){ + i++; + if(f[i]>max)max=f[i]; + } + + if(max>flr[i]){ + oc=oc>>p->shiftoc; + if(oc>=P_BANDS)oc=P_BANDS-1; + if(oc<0)oc=0; + if(vi->tonemaskp) + seed_curve(minseed, + curves[oc], + max, + p->octave[i]-p->firstoc, + p->total_octave_lines, + p->eighth_octave_lines, + dBoffset); + if(vi->peakattp) + seed_peak(maxseed, + att[oc], + max, + p->octave[i]-p->firstoc, + p->eighth_octave_lines, + dBoffset); + } + } +} + +static void bound_loop(vorbis_look_psy *p, + float *f, + float *seeds, + float *flr, + float att){ + long n=p->n,i; + + long off=(p->eighth_octave_lines>>1)+p->firstoc; + long *ocp=p->octave; + + for(i=0;i<n;i++){ + long oc=ocp[i]-off; + float v=f[i]+att; + if(seeds[oc]<v)seeds[oc]=v; + } +} + +static void seed_chase(float *seeds, int linesper, long n){ + long *posstack=alloca(n*sizeof(long)); + float *ampstack=alloca(n*sizeof(float)); + long stack=0; + long pos=0; + long i; + + for(i=0;i<n;i++){ + if(stack<2){ + posstack[stack]=i; + ampstack[stack++]=seeds[i]; + }else{ + while(1){ + if(seeds[i]<ampstack[stack-1]){ + posstack[stack]=i; + ampstack[stack++]=seeds[i]; + break; + }else{ + if(i<posstack[stack-1]+linesper){ + if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] && + i<posstack[stack-2]+linesper){ + /* we completely overlap, making stack-1 irrelevant. pop it */ + stack--; + continue; + } + } + posstack[stack]=i; + ampstack[stack++]=seeds[i]; + break; + + } + } + } + } + + /* the stack now contains only the positions that are relevant. Scan + 'em straight through */ + + for(i=0;i<stack;i++){ + long endpos; + if(i<stack-1 && ampstack[i+1]>ampstack[i]){ + endpos=posstack[i+1]; + }else{ + endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is + discarded in short frames */ + } + if(endpos>n)endpos=n; + for(;pos<endpos;pos++) + seeds[pos]=ampstack[i]; + } + + /* there. Linear time. I now remember this was on a problem set I + had in Grad Skool... I didn't solve it at the time ;-) */ + +} + +/* bleaugh, this is more complicated than it needs to be */ +static void max_seeds(vorbis_look_psy *p,float *minseed,float *maxseed, + float *flr){ + long n=p->total_octave_lines; + int linesper=p->eighth_octave_lines; + long linpos=0; + long pos; + + seed_chase(minseed,linesper,n); /* for masking */ + seed_chase(maxseed,linesper,n); /* for peak att */ + + pos=p->octave[0]-p->firstoc-(linesper>>1); + while(linpos+1<p->n){ + float min=minseed[pos]; + float max=maxseed[pos]; + long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc; + while(pos+1<=end){ + pos++; + if((minseed[pos]>NEGINF && minseed[pos]<min) || min==NEGINF) + min=minseed[pos]; + if(maxseed[pos]>max)max=maxseed[pos]; + } + if(max<min)max=min; + + /* seed scale is log. Floor is linear. Map back to it */ + end=pos+p->firstoc; + for(;linpos<p->n && p->octave[linpos]<=end;linpos++) + if(flr[linpos]<max)flr[linpos]=max; + } + + { + float min=minseed[p->total_octave_lines-1]; + float max=maxseed[p->total_octave_lines-1]; + if(max<min)max=min; + for(;linpos<p->n;linpos++) + if(flr[linpos]<max)flr[linpos]=max; + } + +} + +/* quarter-dB bins */ +#define BIN(x) ((int)((x)*negFour)) +#define BINdB(x) ((x)*negQuarter) +#define BINCOUNT (200*4) +#define LASTBIN (BINCOUNT-1) + +static void bark_noise_median(long n,float *b,float *f,float *noise, + float lowidth,float hiwidth, + int lomin,int himin, + float *thresh,float *off){ + long i=0,lo=0,hi=0; + float bi,threshi; + long median=LASTBIN; + float negFour = -4.0f; + float negQuarter = -0.25f; + + /* these are really integral values, but we store them in floats to + avoid excessive float/int conversions, which GCC and MSVC are + farily poor at optimizing. */ + + float radix[BINCOUNT]; + float countabove=0; + float countbelow=0; + + memset(radix,0,sizeof(radix)); + + for(i=0;i<n;i++){ + /* find new lo/hi */ + bi=b[i]+hiwidth; + for(;hi<n && (hi<i+himin || b[hi]<=bi);hi++){ + int bin=BIN(f[hi]); + if(bin>LASTBIN)bin=LASTBIN; + if(bin<0)bin=0; + radix[bin]++; + if(bin<median) + countabove++; + else + countbelow++; + } + bi=b[i]-lowidth; + for(;lo<i && lo+lomin<i && b[lo]<=bi;lo++){ + int bin=BIN(f[lo]); + if(bin>LASTBIN)bin=LASTBIN; + if(bin<0)bin=0; + radix[bin]--; + if(bin<median) + countabove--; + else + countbelow--; + } + + /* move the median if needed */ + if(countabove+countbelow){ + threshi = thresh[i]*(countabove+countbelow); + + while(threshi>countbelow && median>0){ + median--; + countabove-=radix[median]; + countbelow+=radix[median]; + } + + while(threshi<(countbelow-radix[median]) && + median<LASTBIN){ + countabove+=radix[median]; + countbelow-=radix[median]; + median++; + } + } + noise[i]=BINdB(median)+off[i]; + } + +} + +float _vp_compute_mask(vorbis_look_psy *p, + float *fft, + float *mdct, + float *flr, + float *decay, + float specmax){ + int i,n=p->n; + float localmax=NEGINF; + static int seq=0; + + float *minseed=alloca(sizeof(float)*p->total_octave_lines); + float *maxseed=alloca(sizeof(float)*p->total_octave_lines); + for(i=0;i<p->total_octave_lines;i++)minseed[i]=maxseed[i]=NEGINF; + + /* go to dB scale. Also find the highest peak so we know the limits */ + for(i=0;i<n;i++){ + fft[i]=todB_nn(fft[i]); + if(fft[i]>localmax)localmax=fft[i]; + } + if(specmax<localmax)specmax=localmax; + + + for(i=0;i<n;i++){ + mdct[i]=todB(mdct[i]); + } + + _analysis_output("mdct",seq,mdct,n,0,0); + _analysis_output("fft",seq,fft,n,0,0); + + /* noise masking */ + if(p->vi->noisemaskp){ + bark_noise_median(n,p->bark,mdct,flr, + p->vi->noisewindowlo, + p->vi->noisewindowhi, + p->vi->noisewindowlomin, + p->vi->noisewindowhimin, + p->noisemedian, + p->noiseoffset); + /* suppress any noise curve > specmax+p->vi->noisemaxsupp */ + for(i=0;i<n;i++) + if(flr[i]>specmax+p->vi->noisemaxsupp) + flr[i]=specmax+p->vi->noisemaxsupp; + _analysis_output("noise",seq,flr,n,0,0); + }else{ + for(i=0;i<n;i++)flr[i]=NEGINF; + } + + /* set the ATH (floating below localmax, not global max by a + specified att) */ + if(p->vi->athp){ + float att=localmax+p->vi->ath_adjatt; + if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt; + + for(i=0;i<n;i++){ + float av=p->ath[i]+att; + if(av>flr[i])flr[i]=av; + } + } + + _analysis_output("ath",seq,flr,n,0,0); + + /* tone/peak masking */ + + /* XXX apply decay to the fft here */ + + seed_loop(p,p->tonecurves,p->peakatt,fft,flr,minseed,maxseed,specmax); + bound_loop(p,mdct,maxseed,flr,p->vi->bound_att_dB); + _analysis_output("minseed",seq,minseed,p->total_octave_lines,0,0); + _analysis_output("maxseed",seq,maxseed,p->total_octave_lines,0,0); + max_seeds(p,minseed,maxseed,flr); + _analysis_output("final",seq,flr,n,0,0); + + /* doing this here is clean, but we need to find a faster way to do + it than to just tack it on */ + + if(p->vi->floor_cullp){ + for(i=0;i<n;i++) + if(mdct[i]+6.<flr[i]) + flr[i]=NEGINF; + }else{ + for(i=0;i<n;i++)if(mdct[i]>=flr[i])break; + if(i==n)for(i=0;i<n;i++)flr[i]=NEGINF; + } + + seq++; + + return(specmax); +} + + +/* this applies the floor and (optionally) tries to preserve noise + energy in low resolution portions of the spectrum */ +/* f and flr are *linear* scale, not dB */ +void _vp_apply_floor(vorbis_look_psy *p,float *f, float *flr){ + float *work=alloca(p->n*sizeof(float)); + int j; + + /* subtract the floor */ + for(j=0;j<p->n;j++){ + if(flr[j]<=0) + work[j]=0.f; + else + work[j]=f[j]/flr[j]; + } + + memcpy(f,work,p->n*sizeof(float)); +} + +float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){ + vorbis_info *vi=vd->vi; + codec_setup_info *ci=vi->codec_setup; + int n=ci->blocksizes[vd->W]/2; + float secs=(float)n/vi->rate; + + amp+=secs*ci->ampmax_att_per_sec; + if(amp<-9999)amp=-9999; + return(amp); +} + + + diff --git a/lib/psy.h b/lib/psy.h new file mode 100644 index 00000000..1e0afccd --- /dev/null +++ b/lib/psy.h @@ -0,0 +1,101 @@ +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * + * * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * + * by the XIPHOPHORUS Company http://www.xiph.org/ * + + ******************************************************************** + + function: random psychoacoustics (not including preecho) + last mod: $Id: psy.h,v 1.19.4.1 2001/04/29 22:21:05 xiphmont Exp $ + + ********************************************************************/ + +#ifndef _V_PSY_H_ +#define _V_PSY_H_ +#include "smallft.h" + +#ifndef EHMER_MAX +#define EHMER_MAX 56 +#endif + +/* psychoacoustic setup ********************************************/ +#define MAX_BARK 27 +#define P_BANDS 17 +#define P_LEVELS 11 +typedef struct vorbis_info_psy{ + int athp; + int decayp; + + float ath_adjatt; + float ath_maxatt; + + int eighth_octave_lines; + + /* 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 */ + /* x: 63 88 125 175 250 350 500 700 1k 1.4k 2k 2.8k 4k 5.6k 8k 11.5k 16k Hz */ + /* y: 0 10 20 30 40 50 60 70 80 90 100 dB */ + + int tonemaskp; + float toneatt[P_BANDS][P_LEVELS]; + + int peakattp; + float peakatt[P_BANDS][P_LEVELS]; + + int noisemaskp; + float noisemaxsupp; + float noisewindowlo; + float noisewindowhi; + int noisewindowlomin; + int noisewindowhimin; + float noisemedian[P_BANDS*2]; + + float max_curve_dB; + float bound_att_dB; + + int floor_cullp; +} vorbis_info_psy; + +typedef struct { + int n; + struct vorbis_info_psy *vi; + + float ***tonecurves; + float **peakatt; + float *noisemedian; + float *noiseoffset; + + float *ath; + long *octave; /* in n.ocshift format */ + float *bark; + + long firstoc; + long shiftoc; + int eighth_octave_lines; /* power of two, please */ + int total_octave_lines; + +} vorbis_look_psy; + +extern void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate); +extern void _vp_psy_clear(vorbis_look_psy *p); +extern void *_vi_psy_dup(void *source); + +extern void _vi_psy_free(vorbis_info_psy *i); +extern vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i); + +extern float _vp_compute_mask(vorbis_look_psy *p, + float *fft, + float *mdct, + float *floor, + float *decay, + float prev_maxamp); +extern void _vp_apply_floor(vorbis_look_psy *p,float *f,float *flr); +extern float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd); + +#endif + + |