psytune.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357
  1. /********************************************************************
  2. * *
  3. * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
  4. * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
  5. * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH *
  6. * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
  7. * *
  8. * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
  9. * by Monty <monty@xiph.org> and the XIPHOPHORUS Company *
  10. * http://www.xiph.org/ *
  11. * *
  12. ********************************************************************
  13. function: simple utility that runs audio through the psychoacoustics
  14. without encoding
  15. last mod: $Id: psytune.c,v 1.7.2.2 2000/11/04 06:43:50 xiphmont Exp $
  16. ********************************************************************/
  17. #include <stdio.h>
  18. #include <stdlib.h>
  19. #include <string.h>
  20. #include <math.h>
  21. #include "vorbis/codec.h"
  22. #include "os.h"
  23. #include "psy.h"
  24. #include "mdct.h"
  25. #include "window.h"
  26. #include "scales.h"
  27. #include "lpc.h"
  28. #include "lsp.h"
  29. static vorbis_info_psy _psy_set0={
  30. 1,/*athp*/
  31. 1,/*decayp*/
  32. 1,/*smoothp*/
  33. 0,.2,
  34. -100.,
  35. -140.,
  36. /* 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 */
  37. /* x: 63 88 125 175 250 350 500 700 1k 1.4k 2k 2.8k 4k 5.6k 8k 11.5k 16k Hz */
  38. /* y: 0 10 20 30 40 50 60 70 80 90 100 dB */
  39. 1,/* tonemaskp */
  40. /* 0 10 20 30 40 50 60 70 80 90 100 */
  41. {{-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*63*/
  42. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*88*/
  43. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*125*/
  44. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*175*/
  45. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*250*/
  46. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*350*/
  47. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*500*/
  48. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*700*/
  49. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1000*/
  50. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1400*/
  51. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2000*/
  52. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2800*/
  53. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*4000*/
  54. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*5600*/
  55. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*8000*/
  56. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*11500*/
  57. {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*16000*/
  58. },
  59. 1,/* peakattp */
  60. {{-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*63*/
  61. {-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*88*/
  62. {-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*125*/
  63. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  64. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  65. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  66. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  67. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  68. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  69. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  70. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  71. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  72. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  73. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  74. { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-16.,-17.}, /*8000*/
  75. { -6., -7., -9., -9., -9., -9.,-10.,-11.,-12.,-13.,-14.}, /*11500*/
  76. { -6., -6., -9., -9., -9., -9., -9., -9.,-10.,-11.,-12.}, /*16000*/
  77. },
  78. 1,/*noisemaskp */
  79. /* 0 10 20 30 40 50 60 70 80 90 100 */
  80. {{-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*63*/
  81. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*88*/
  82. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*125*/
  83. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*175*/
  84. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*250*/
  85. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  86. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  87. {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  88. { -6., -6., -6., -6., -6., -6., -6., -6., -6., -6., -6.}, /*2000*/
  89. { -6., -6., -6., -6., -6., -6., -6., -6., -6., -6., -6.}, /*2000*/
  90. { -6., -6., -6., -6., -6., -6., -6., -6., -6., -6., -6.}, /*2000*/
  91. { -6., -6., -6., -6., -6., -6., -6., -6., -6., -6., -6.}, /*2800*/
  92. { -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.}, /*4000*/
  93. { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, /*5600*/
  94. { 4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.}, /*8000*/
  95. { 4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.}, /*11500*/
  96. { 4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.}, /*16000*/
  97. },
  98. 110.,
  99. -0, -.005, /* attack/decay control */
  100. };
  101. static int noisy=0;
  102. void analysis(char *base,int i,float *v,int n,int bark,int dB){
  103. if(noisy){
  104. int j;
  105. FILE *of;
  106. char buffer[80];
  107. sprintf(buffer,"%s_%d.m",base,i);
  108. of=fopen(buffer,"w");
  109. for(j=0;j<n;j++){
  110. if(dB && v[j]==0)
  111. fprintf(of,"\n\n");
  112. else{
  113. if(bark)
  114. fprintf(of,"%g ",toBARK(22050.*j/n));
  115. else
  116. fprintf(of,"%g ",(float)j);
  117. if(dB){
  118. fprintf(of,"%g\n",todB(fabs(v[j])));
  119. }else{
  120. fprintf(of,"%g\n",v[j]);
  121. }
  122. }
  123. }
  124. fclose(of);
  125. }
  126. }
  127. typedef struct {
  128. long n;
  129. int ln;
  130. int m;
  131. int *linearmap;
  132. vorbis_info_floor0 *vi;
  133. lpc_lookup lpclook;
  134. } vorbis_look_floor0;
  135. long granulepos=0;
  136. /* hacked from floor0.c */
  137. static void floorinit(vorbis_look_floor0 *look,int n,int m,int ln){
  138. int j;
  139. float scale;
  140. look->m=m;
  141. look->n=n;
  142. look->ln=ln;
  143. lpc_init(&look->lpclook,look->ln,look->m);
  144. scale=look->ln/toBARK(22050.);
  145. look->linearmap=_ogg_malloc(look->n*sizeof(int));
  146. for(j=0;j<look->n;j++){
  147. int val=floor( toBARK(22050./n*j) *scale);
  148. if(val>look->ln)val=look->ln;
  149. look->linearmap[j]=val;
  150. }
  151. }
  152. int main(int argc,char *argv[]){
  153. int eos=0;
  154. float nonz=0.;
  155. float acc=0.;
  156. float tot=0.;
  157. int framesize=2048;
  158. int order=32;
  159. int map=256;
  160. float *pcm[2],*out[2],*window,*decay[2],*lpc,*floor;
  161. signed char *buffer,*buffer2;
  162. mdct_lookup m_look;
  163. vorbis_look_psy p_look;
  164. long i,j,k;
  165. vorbis_look_floor0 floorlook;
  166. int ath=0;
  167. int decayp=0;
  168. argv++;
  169. while(*argv){
  170. if(*argv[0]=='-'){
  171. /* option */
  172. if(argv[0][1]=='v'){
  173. noisy=0;
  174. }
  175. if(argv[0][1]=='A'){
  176. ath=0;
  177. }
  178. if(argv[0][1]=='D'){
  179. decayp=0;
  180. }
  181. if(argv[0][1]=='X'){
  182. ath=0;
  183. decayp=0;
  184. }
  185. }else
  186. if(*argv[0]=='+'){
  187. /* option */
  188. if(argv[0][1]=='v'){
  189. noisy=1;
  190. }
  191. if(argv[0][1]=='A'){
  192. ath=1;
  193. }
  194. if(argv[0][1]=='D'){
  195. decayp=1;
  196. }
  197. if(argv[0][1]=='X'){
  198. ath=1;
  199. decayp=1;
  200. }
  201. }else
  202. framesize=atoi(argv[0]);
  203. argv++;
  204. }
  205. pcm[0]=_ogg_malloc(framesize*sizeof(float));
  206. pcm[1]=_ogg_malloc(framesize*sizeof(float));
  207. out[0]=_ogg_calloc(framesize/2,sizeof(float));
  208. out[1]=_ogg_calloc(framesize/2,sizeof(float));
  209. decay[0]=_ogg_calloc(framesize/2,sizeof(float));
  210. decay[1]=_ogg_calloc(framesize/2,sizeof(float));
  211. floor=_ogg_malloc(framesize*sizeof(float));
  212. lpc=_ogg_malloc(order*sizeof(float));
  213. buffer=_ogg_malloc(framesize*4);
  214. buffer2=buffer+framesize*2;
  215. window=_vorbis_window(0,framesize,framesize/2,framesize/2);
  216. mdct_init(&m_look,framesize);
  217. _vp_psy_init(&p_look,&_psy_set0,framesize/2,44100);
  218. floorinit(&floorlook,framesize/2,order,map);
  219. for(i=0;i<P_BANDS;i++)
  220. for(j=0;j<P_LEVELS;j++)
  221. analysis("Ptonecurve",i*100+j,p_look.tonecurves[i][j],EHMER_MAX,0,1);
  222. /* we cheat on the WAV header; we just bypass 44 bytes and never
  223. verify that it matches 16bit/stereo/44.1kHz. */
  224. fread(buffer,1,44,stdin);
  225. fwrite(buffer,1,44,stdout);
  226. memset(buffer,0,framesize*2);
  227. analysis("window",0,window,framesize,0,0);
  228. fprintf(stderr,"Processing for frame size %d...\n",framesize);
  229. while(!eos){
  230. long bytes=fread(buffer2,1,framesize*2,stdin);
  231. if(bytes<framesize*2)
  232. memset(buffer2+bytes,0,framesize*2-bytes);
  233. if(bytes!=0){
  234. /* uninterleave samples */
  235. for(i=0;i<framesize;i++){
  236. pcm[0][i]=((buffer[i*4+1]<<8)|
  237. (0x00ff&(int)buffer[i*4]))/32768.;
  238. pcm[1][i]=((buffer[i*4+3]<<8)|
  239. (0x00ff&(int)buffer[i*4+2]))/32768.;
  240. }
  241. for(i=0;i<2;i++){
  242. float amp;
  243. analysis("pre",granulepos,pcm[i],framesize,0,0);
  244. /* do the psychacoustics */
  245. for(j=0;j<framesize;j++)
  246. pcm[i][j]*=window[j];
  247. mdct_forward(&m_look,pcm[i],pcm[i]);
  248. analysis("mdct",granulepos,pcm[i],framesize/2,1,1);
  249. _vp_compute_mask(&p_look,pcm[i],floor,decay[i]);
  250. analysis("floor",frameno,floor,framesize/2,1,1);
  251. analysis("decay",frameno,decay[i],framesize/2,1,1);
  252. _vp_apply_floor(&p_look,pcm[i],floor);
  253. /*r(j=0;j<framesize/2;j++)
  254. if(fabs(pcm[i][j])<1.)pcm[i][j]=0;*/
  255. analysis("quant",granulepos,pcm[i],framesize/2,1,1);
  256. /* re-add floor */
  257. for(j=0;j<framesize/2;j++){
  258. float val=rint(pcm[i][j]);
  259. tot++;
  260. if(val){
  261. nonz++;
  262. acc+=log(fabs(val)*2.+1.)/log(2);
  263. pcm[i][j]=val*floor[j];
  264. }else{
  265. pcm[i][j]=0;
  266. }
  267. }
  268. analysis("final",granulepos,pcm[i],framesize/2,1,1);
  269. /* take it back to time */
  270. mdct_backward(&m_look,pcm[i],pcm[i]);
  271. for(j=0;j<framesize/2;j++)
  272. out[i][j]+=pcm[i][j]*window[j];
  273. granulepos++;
  274. }
  275. /* write data. Use the part of buffer we're about to shift out */
  276. for(i=0;i<2;i++){
  277. char *ptr=buffer+i*2;
  278. float *mono=out[i];
  279. for(j=0;j<framesize/2;j++){
  280. int val=mono[j]*32767.;
  281. /* might as well guard against clipping */
  282. if(val>32767)val=32767;
  283. if(val<-32768)val=-32768;
  284. ptr[0]=val&0xff;
  285. ptr[1]=(val>>8)&0xff;
  286. ptr+=4;
  287. }
  288. }
  289. fprintf(stderr,"*");
  290. fwrite(buffer,1,framesize*2,stdout);
  291. memmove(buffer,buffer2,framesize*2);
  292. for(i=0;i<2;i++){
  293. for(j=0,k=framesize/2;j<framesize/2;j++,k++)
  294. out[i][j]=pcm[i][k]*window[k];
  295. }
  296. }else
  297. eos=1;
  298. }
  299. fprintf(stderr,"average raw bits of entropy: %.03g/sample\n",acc/tot);
  300. fprintf(stderr,"average nonzero samples: %.03g/%d\n",nonz/tot*framesize/2,
  301. framesize/2);
  302. fprintf(stderr,"Done\n\n");
  303. return 0;
  304. }