1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: psychoacoustics not including preecho
14 last mod: $Id: psy.c,v 1.58 2001/12/18 01:58:15 segher Exp $
16 ********************************************************************/
21 #include "vorbis/codec.h"
22 #include "codec_internal.h"
32 #define NEGINF -9999.f
34 /* Why Bark scale for encoding but not masking computation? Because
35 masking has a strong harmonic dependency */
37 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
38 codec_setup_info *ci=vi->codec_setup;
39 vorbis_info_psy_global *gi=&ci->psy_g_param;
40 vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
42 look->channels=vi->channels;
49 void _vp_global_free(vorbis_look_psy_global *look){
51 memset(look,0,sizeof(*look));
56 void _vi_gpsy_free(vorbis_info_psy_global *i){
58 memset(i,0,sizeof(*i));
63 void _vi_psy_free(vorbis_info_psy *i){
65 memset(i,0,sizeof(*i));
70 vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i){
71 vorbis_info_psy *ret=_ogg_malloc(sizeof(*ret));
72 memcpy(ret,i,sizeof(*ret));
76 /* Set up decibel threshold slopes on a Bark frequency scale */
77 /* ATH is the only bit left on a Bark scale. No reason to change it
79 static void set_curve(float *ref,float *c,int n, float crate){
82 for(i=0;i<MAX_BARK-1;i++){
83 int endpos=rint(fromBARK(i+1)*2*n/crate);
86 float delta=(ref[i+1]-base)/(endpos-j);
87 for(;j<endpos && j<n;j++){
95 static void min_curve(float *c,
98 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
100 static void max_curve(float *c,
103 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
106 static void attenuate_curve(float *c,float att){
108 for(i=0;i<EHMER_MAX;i++)
112 static void interp_curve(float *c,float *c1,float *c2,float del){
114 for(i=0;i<EHMER_MAX;i++)
115 c[i]=c2[i]*del+c1[i]*(1.f-del);
118 extern int analysis_noisy;
119 static void setup_curve(float **c,
123 float ath[EHMER_MAX];
124 float tempc[P_LEVELS][EHMER_MAX];
125 float *ATH=ATH_Bark_dB_lspconservative; /* just for limiting here */
127 memcpy(c[0]+2,c[4]+2,sizeof(*c[0])*EHMER_MAX);
128 memcpy(c[2]+2,c[4]+2,sizeof(*c[2])*EHMER_MAX);
130 /* we add back in the ATH to avoid low level curves falling off to
131 -infinity and unnecessarily cutting off high level curves in the
132 curve limiting (last step). But again, remember... a half-band's
133 settings must be valid over the whole band, and it's better to
134 mask too little than too much, so be pessimistical. */
136 for(i=0;i<EHMER_MAX;i++){
137 float oc_min=band*.5+(i-EHMER_OFFSET)*.125;
138 float oc_max=band*.5+(i-EHMER_OFFSET+1)*.125;
139 float bark=toBARK(fromOC(oc_min));
140 int ibark=floor(bark);
141 float del=bark-ibark;
142 float ath_min,ath_max;
145 ath_min=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
149 bark=toBARK(fromOC(oc_max));
154 ath_max=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
158 ath[i]=min(ath_min,ath_max);
161 /* The c array comes in as dB curves at 20 40 60 80 100 dB.
162 interpolate intermediate dB curves */
163 for(i=1;i<P_LEVELS;i+=2){
164 interp_curve(c[i]+2,c[i-1]+2,c[i+1]+2,.5);
167 /* normalize curves so the driving amplitude is 0dB */
168 /* make temp curves with the ATH overlayed */
169 for(i=0;i<P_LEVELS;i++){
170 attenuate_curve(c[i]+2,curveatt_dB[i]);
171 memcpy(tempc[i],ath,EHMER_MAX*sizeof(*tempc[i]));
172 attenuate_curve(tempc[i],-i*10.f);
173 max_curve(tempc[i],c[i]+2);
176 /* Now limit the louder curves.
178 the idea is this: We don't know what the playback attenuation
179 will be; 0dB SL moves every time the user twiddles the volume
180 knob. So that means we have to use a single 'most pessimal' curve
181 for all masking amplitudes, right? Wrong. The *loudest* sound
182 can be in (we assume) a range of ...+100dB] SL. However, sounds
183 20dB down will be in a range ...+80], 40dB down is from ...+60],
186 for(j=1;j<P_LEVELS;j++){
187 min_curve(tempc[j],tempc[j-1]);
188 min_curve(c[j]+2,tempc[j]);
192 for(j=0;j<P_LEVELS;j++){
194 for(i=0;i<EHMER_OFFSET;i++)
195 if(c[j][i+2]>-200.f)break;
198 for(i=EHMER_MAX-1;i>EHMER_OFFSET+1;i--)
206 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
207 vorbis_info_psy_global *gi,int n,long rate){
208 long i,j,k,lo=-99,hi=0;
210 memset(p,0,sizeof(*p));
213 p->eighth_octave_lines=gi->eighth_octave_lines;
214 p->shiftoc=rint(log(gi->eighth_octave_lines*8)/log(2))-1;
216 p->firstoc=toOC(.25f*rate/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
217 maxoc=toOC((n*.5f-.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
218 p->total_octave_lines=maxoc-p->firstoc+1;
221 p->ath=_ogg_malloc(n*sizeof(*p->ath));
222 p->octave=_ogg_malloc(n*sizeof(*p->octave));
223 p->bark=_ogg_malloc(n*sizeof(*p->bark));
228 /* set up the lookups for a given blocksize and sample rate */
230 set_curve(vi->ath, p->ath,n,rate);
232 float bark=toBARK(rate/(2*n)*i);
234 for(;lo+vi->noisewindowlomin<i &&
235 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
237 for(;hi<n && (hi<i+vi->noisewindowhimin ||
238 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
240 p->bark[i]=(lo<<16)+hi;
245 p->octave[i]=toOC((i*.5f+.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
247 p->tonecurves=_ogg_malloc(P_BANDS*sizeof(*p->tonecurves));
248 p->noisethresh=_ogg_malloc(n*sizeof(*p->noisethresh));
249 p->noiseoffset=_ogg_malloc(n*sizeof(*p->noiseoffset));
250 for(i=0;i<P_BANDS;i++)
251 p->tonecurves[i]=_ogg_malloc(P_LEVELS*sizeof(*p->tonecurves[i]));
253 for(i=0;i<P_BANDS;i++)
254 for(j=0;j<P_LEVELS;j++)
255 p->tonecurves[i][j]=_ogg_malloc((EHMER_MAX+2)*sizeof(*p->tonecurves[i][j]));
258 /* OK, yeah, this was a silly way to do it */
259 memcpy(p->tonecurves[0][4]+2,tone_125_40dB_SL,sizeof(*p->tonecurves[0][4])*EHMER_MAX);
260 memcpy(p->tonecurves[0][6]+2,tone_125_60dB_SL,sizeof(*p->tonecurves[0][6])*EHMER_MAX);
261 memcpy(p->tonecurves[0][8]+2,tone_125_80dB_SL,sizeof(*p->tonecurves[0][8])*EHMER_MAX);
262 memcpy(p->tonecurves[0][10]+2,tone_125_100dB_SL,sizeof(*p->tonecurves[0][10])*EHMER_MAX);
264 memcpy(p->tonecurves[2][4]+2,tone_125_40dB_SL,sizeof(*p->tonecurves[2][4])*EHMER_MAX);
265 memcpy(p->tonecurves[2][6]+2,tone_125_60dB_SL,sizeof(*p->tonecurves[2][6])*EHMER_MAX);
266 memcpy(p->tonecurves[2][8]+2,tone_125_80dB_SL,sizeof(*p->tonecurves[2][8])*EHMER_MAX);
267 memcpy(p->tonecurves[2][10]+2,tone_125_100dB_SL,sizeof(*p->tonecurves[2][10])*EHMER_MAX);
269 memcpy(p->tonecurves[4][4]+2,tone_250_40dB_SL,sizeof(*p->tonecurves[4][4])*EHMER_MAX);
270 memcpy(p->tonecurves[4][6]+2,tone_250_60dB_SL,sizeof(*p->tonecurves[4][6])*EHMER_MAX);
271 memcpy(p->tonecurves[4][8]+2,tone_250_80dB_SL,sizeof(*p->tonecurves[4][8])*EHMER_MAX);
272 memcpy(p->tonecurves[4][10]+2,tone_250_100dB_SL,sizeof(*p->tonecurves[4][10])*EHMER_MAX);
274 memcpy(p->tonecurves[6][4]+2,tone_500_40dB_SL,sizeof(*p->tonecurves[6][4])*EHMER_MAX);
275 memcpy(p->tonecurves[6][6]+2,tone_500_60dB_SL,sizeof(*p->tonecurves[6][6])*EHMER_MAX);
276 memcpy(p->tonecurves[6][8]+2,tone_500_80dB_SL,sizeof(*p->tonecurves[6][8])*EHMER_MAX);
277 memcpy(p->tonecurves[6][10]+2,tone_500_100dB_SL,sizeof(*p->tonecurves[6][10])*EHMER_MAX);
279 memcpy(p->tonecurves[8][4]+2,tone_1000_40dB_SL,sizeof(*p->tonecurves[8][4])*EHMER_MAX);
280 memcpy(p->tonecurves[8][6]+2,tone_1000_60dB_SL,sizeof(*p->tonecurves[8][6])*EHMER_MAX);
281 memcpy(p->tonecurves[8][8]+2,tone_1000_80dB_SL,sizeof(*p->tonecurves[8][8])*EHMER_MAX);
282 memcpy(p->tonecurves[8][10]+2,tone_1000_100dB_SL,sizeof(*p->tonecurves[8][10])*EHMER_MAX);
284 memcpy(p->tonecurves[10][4]+2,tone_2000_40dB_SL,sizeof(*p->tonecurves[10][4])*EHMER_MAX);
285 memcpy(p->tonecurves[10][6]+2,tone_2000_60dB_SL,sizeof(*p->tonecurves[10][6])*EHMER_MAX);
286 memcpy(p->tonecurves[10][8]+2,tone_2000_80dB_SL,sizeof(*p->tonecurves[10][8])*EHMER_MAX);
287 memcpy(p->tonecurves[10][10]+2,tone_2000_100dB_SL,sizeof(*p->tonecurves[10][10])*EHMER_MAX);
289 memcpy(p->tonecurves[12][4]+2,tone_4000_40dB_SL,sizeof(*p->tonecurves[12][4])*EHMER_MAX);
290 memcpy(p->tonecurves[12][6]+2,tone_4000_60dB_SL,sizeof(*p->tonecurves[12][6])*EHMER_MAX);
291 memcpy(p->tonecurves[12][8]+2,tone_4000_80dB_SL,sizeof(*p->tonecurves[12][8])*EHMER_MAX);
292 memcpy(p->tonecurves[12][10]+2,tone_4000_100dB_SL,sizeof(*p->tonecurves[12][10])*EHMER_MAX);
294 memcpy(p->tonecurves[14][4]+2,tone_8000_40dB_SL,sizeof(*p->tonecurves[14][4])*EHMER_MAX);
295 memcpy(p->tonecurves[14][6]+2,tone_8000_60dB_SL,sizeof(*p->tonecurves[14][6])*EHMER_MAX);
296 memcpy(p->tonecurves[14][8]+2,tone_8000_80dB_SL,sizeof(*p->tonecurves[14][8])*EHMER_MAX);
297 memcpy(p->tonecurves[14][10]+2,tone_8000_100dB_SL,sizeof(*p->tonecurves[14][10])*EHMER_MAX);
299 memcpy(p->tonecurves[16][4]+2,tone_8000_40dB_SL,sizeof(*p->tonecurves[16][4])*EHMER_MAX);
300 memcpy(p->tonecurves[16][6]+2,tone_8000_60dB_SL,sizeof(*p->tonecurves[16][6])*EHMER_MAX);
301 memcpy(p->tonecurves[16][8]+2,tone_8000_80dB_SL,sizeof(*p->tonecurves[16][8])*EHMER_MAX);
302 memcpy(p->tonecurves[16][10]+2,tone_8000_100dB_SL,sizeof(*p->tonecurves[16][10])*EHMER_MAX);
304 for(i=0;i<P_BANDS;i+=2)
305 for(j=4;j<P_LEVELS;j+=2)
306 for(k=2;k<EHMER_MAX+2;k++)
307 p->tonecurves[i][j][k]+=vi->tone_masteratt;
309 /* interpolate curves between */
310 for(i=1;i<P_BANDS;i+=2)
311 for(j=4;j<P_LEVELS;j+=2){
312 memcpy(p->tonecurves[i][j]+2,p->tonecurves[i-1][j]+2,EHMER_MAX*sizeof(*p->tonecurves[i][j]));
313 /*interp_curve(p->tonecurves[i][j],
314 p->tonecurves[i-1][j],
315 p->tonecurves[i+1][j],.5);*/
316 min_curve(p->tonecurves[i][j]+2,p->tonecurves[i+1][j]+2);
319 /* set up the final curves */
320 for(i=0;i<P_BANDS;i++)
321 setup_curve(p->tonecurves[i],i,vi->toneatt.block[i]);
323 for(i=0;i<P_LEVELS;i++)
324 _analysis_output("curve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
325 for(i=0;i<P_LEVELS;i++)
326 _analysis_output("curve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
327 for(i=0;i<P_LEVELS;i++)
328 _analysis_output("curve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
329 for(i=0;i<P_LEVELS;i++)
330 _analysis_output("curve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
331 for(i=0;i<P_LEVELS;i++)
332 _analysis_output("curve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
333 for(i=0;i<P_LEVELS;i++)
334 _analysis_output("curve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
335 for(i=0;i<P_LEVELS;i++)
336 _analysis_output("curve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
337 for(i=0;i<P_LEVELS;i++)
338 _analysis_output("curve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
339 for(i=0;i<P_LEVELS;i++)
340 _analysis_output("curve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
341 for(i=0;i<P_LEVELS;i++)
342 _analysis_output("curve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
343 for(i=0;i<P_LEVELS;i++)
344 _analysis_output("curve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
345 for(i=0;i<P_LEVELS;i++)
346 _analysis_output("curve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
347 for(i=0;i<P_LEVELS;i++)
348 _analysis_output("curve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
349 for(i=0;i<P_LEVELS;i++)
350 _analysis_output("curve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
351 for(i=0;i<P_LEVELS;i++)
352 _analysis_output("curve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
353 for(i=0;i<P_LEVELS;i++)
354 _analysis_output("curve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
355 for(i=0;i<P_LEVELS;i++)
356 _analysis_output("curve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
359 /* value limit the tonal masking curves; the peakatt not only
360 optionally specifies maximum dynamic depth, but also
361 limits the masking curves to a minimum depth */
362 for(i=0;i<P_BANDS;i++)
363 for(j=0;j<P_LEVELS;j++){
364 for(k=2;k<EHMER_OFFSET+2+vi->curvelimitp;k++)
365 if(p->tonecurves[i][j][k]> vi->peakatt.block[i][j])
366 p->tonecurves[i][j][k]= vi->peakatt.block[i][j];
372 for(i=0;i<P_LEVELS;i++)
373 _analysis_output("licurve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
374 for(i=0;i<P_LEVELS;i++)
375 _analysis_output("licurve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
376 for(i=0;i<P_LEVELS;i++)
377 _analysis_output("licurve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
378 for(i=0;i<P_LEVELS;i++)
379 _analysis_output("licurve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
380 for(i=0;i<P_LEVELS;i++)
381 _analysis_output("licurve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
382 for(i=0;i<P_LEVELS;i++)
383 _analysis_output("licurve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
384 for(i=0;i<P_LEVELS;i++)
385 _analysis_output("licurve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
386 for(i=0;i<P_LEVELS;i++)
387 _analysis_output("licurve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
388 for(i=0;i<P_LEVELS;i++)
389 _analysis_output("licurve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
390 for(i=0;i<P_LEVELS;i++)
391 _analysis_output("licurve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
392 for(i=0;i<P_LEVELS;i++)
393 _analysis_output("licurve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
394 for(i=0;i<P_LEVELS;i++)
395 _analysis_output("licurve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
396 for(i=0;i<P_LEVELS;i++)
397 _analysis_output("licurve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
398 for(i=0;i<P_LEVELS;i++)
399 _analysis_output("licurve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
400 for(i=0;i<P_LEVELS;i++)
401 _analysis_output("licurve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
402 for(i=0;i<P_LEVELS;i++)
403 _analysis_output("licurve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
404 for(i=0;i<P_LEVELS;i++)
405 _analysis_output("licurve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
407 if(vi->peakattp) /* we limit maximum depth only optionally */
408 for(i=0;i<P_BANDS;i++)
409 for(j=0;j<P_LEVELS;j++)
410 if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->peakatt.block[i][j])
411 p->tonecurves[i][j][EHMER_OFFSET+2]= vi->peakatt.block[i][j];
413 for(i=0;i<P_LEVELS;i++)
414 _analysis_output("pcurve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
415 for(i=0;i<P_LEVELS;i++)
416 _analysis_output("pcurve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
417 for(i=0;i<P_LEVELS;i++)
418 _analysis_output("pcurve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
419 for(i=0;i<P_LEVELS;i++)
420 _analysis_output("pcurve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
421 for(i=0;i<P_LEVELS;i++)
422 _analysis_output("pcurve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
423 for(i=0;i<P_LEVELS;i++)
424 _analysis_output("pcurve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
425 for(i=0;i<P_LEVELS;i++)
426 _analysis_output("pcurve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
427 for(i=0;i<P_LEVELS;i++)
428 _analysis_output("pcurve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
429 for(i=0;i<P_LEVELS;i++)
430 _analysis_output("pcurve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
431 for(i=0;i<P_LEVELS;i++)
432 _analysis_output("pcurve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
433 for(i=0;i<P_LEVELS;i++)
434 _analysis_output("pcurve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
435 for(i=0;i<P_LEVELS;i++)
436 _analysis_output("pcurve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
437 for(i=0;i<P_LEVELS;i++)
438 _analysis_output("pcurve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
439 for(i=0;i<P_LEVELS;i++)
440 _analysis_output("pcurve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
441 for(i=0;i<P_LEVELS;i++)
442 _analysis_output("pcurve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
443 for(i=0;i<P_LEVELS;i++)
444 _analysis_output("pcurve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
445 for(i=0;i<P_LEVELS;i++)
446 _analysis_output("pcurve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
448 /* but guarding is mandatory */
449 for(i=0;i<P_BANDS;i++)
450 for(j=0;j<P_LEVELS;j++)
451 if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->tone_guard)
452 p->tonecurves[i][j][EHMER_OFFSET+2]= vi->tone_guard;
454 for(i=0;i<P_LEVELS;i++)
455 _analysis_output("fcurve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
456 for(i=0;i<P_LEVELS;i++)
457 _analysis_output("fcurve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
458 for(i=0;i<P_LEVELS;i++)
459 _analysis_output("fcurve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
460 for(i=0;i<P_LEVELS;i++)
461 _analysis_output("fcurve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
462 for(i=0;i<P_LEVELS;i++)
463 _analysis_output("fcurve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
464 for(i=0;i<P_LEVELS;i++)
465 _analysis_output("fcurve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
466 for(i=0;i<P_LEVELS;i++)
467 _analysis_output("fcurve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
468 for(i=0;i<P_LEVELS;i++)
469 _analysis_output("fcurve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
470 for(i=0;i<P_LEVELS;i++)
471 _analysis_output("fcurve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
472 for(i=0;i<P_LEVELS;i++)
473 _analysis_output("fcurve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
474 for(i=0;i<P_LEVELS;i++)
475 _analysis_output("fcurve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
476 for(i=0;i<P_LEVELS;i++)
477 _analysis_output("fcurve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
478 for(i=0;i<P_LEVELS;i++)
479 _analysis_output("fcurve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
480 for(i=0;i<P_LEVELS;i++)
481 _analysis_output("fcurve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
482 for(i=0;i<P_LEVELS;i++)
483 _analysis_output("fcurve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
484 for(i=0;i<P_LEVELS;i++)
485 _analysis_output("fcurve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
486 for(i=0;i<P_LEVELS;i++)
487 _analysis_output("fcurve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
489 /* set up rolling noise median */
491 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
495 if(halfoc<0)halfoc=0;
496 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
497 inthalfoc=(int)halfoc;
498 del=halfoc-inthalfoc;
500 p->vi->noiseoff[inthalfoc]*(1.-del) +
501 p->vi->noiseoff[inthalfoc+1]*del;
505 _analysis_output("noiseoff",0,p->noiseoffset,n,1,0);
506 _analysis_output("noisethresh",0,p->noisethresh,n,1,0);
511 void _vp_psy_clear(vorbis_look_psy *p){
514 if(p->ath)_ogg_free(p->ath);
515 if(p->octave)_ogg_free(p->octave);
516 if(p->bark)_ogg_free(p->bark);
518 for(i=0;i<P_BANDS;i++){
519 for(j=0;j<P_LEVELS;j++){
520 _ogg_free(p->tonecurves[i][j]);
522 _ogg_free(p->tonecurves[i]);
524 _ogg_free(p->tonecurves);
526 _ogg_free(p->noiseoffset);
527 _ogg_free(p->noisethresh);
528 memset(p,0,sizeof(*p));
532 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
533 static void seed_curve(float *seed,
534 const float **curves,
537 int linesper,float dBoffset){
540 const float *posts,*curve;
542 int choice=(int)((amp+dBoffset)*.1f);
543 choice=max(choice,0);
544 choice=min(choice,P_LEVELS-1);
545 posts=curves[choice];
548 seedptr=oc+(posts[0]-16)*linesper-(linesper>>1);
550 for(i=posts[0];i<post1;i++){
552 float lin=amp+curve[i];
553 if(seed[seedptr]<lin)seed[seedptr]=lin;
560 static void seed_loop(vorbis_look_psy *p,
561 const float ***curves,
566 vorbis_info_psy *vi=p->vi;
568 float dBoffset=vi->max_curve_dB-specmax;
570 /* prime the working vector with peak values */
574 long oc=p->octave[i];
575 while(i+1<n && p->octave[i+1]==oc){
577 if(f[i]>max)max=f[i];
582 if(oc>=P_BANDS)oc=P_BANDS-1;
587 p->octave[i]-p->firstoc,
588 p->total_octave_lines,
589 p->eighth_octave_lines,
595 static void seed_chase(float *seeds, int linesper, long n){
596 long *posstack=alloca(n*sizeof(*posstack));
597 float *ampstack=alloca(n*sizeof(*ampstack));
605 ampstack[stack++]=seeds[i];
608 if(seeds[i]<ampstack[stack-1]){
610 ampstack[stack++]=seeds[i];
613 if(i<posstack[stack-1]+linesper){
614 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
615 i<posstack[stack-2]+linesper){
616 /* we completely overlap, making stack-1 irrelevant. pop it */
622 ampstack[stack++]=seeds[i];
630 /* the stack now contains only the positions that are relevant. Scan
631 'em straight through */
633 for(i=0;i<stack;i++){
635 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
636 endpos=posstack[i+1];
638 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
639 discarded in short frames */
641 if(endpos>n)endpos=n;
642 for(;pos<endpos;pos++)
643 seeds[pos]=ampstack[i];
646 /* there. Linear time. I now remember this was on a problem set I
647 had in Grad Skool... I didn't solve it at the time ;-) */
651 /* bleaugh, this is more complicated than it needs to be */
652 static void max_seeds(vorbis_look_psy *p,
653 vorbis_look_psy_global *g,
657 long n=p->total_octave_lines;
658 int linesper=p->eighth_octave_lines;
662 seed_chase(seed,linesper,n); /* for masking */
664 pos=p->octave[0]-p->firstoc-(linesper>>1);
665 while(linpos+1<p->n){
666 float minV=seed[pos];
667 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
668 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
671 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
675 /* seed scale is log. Floor is linear. Map back to it */
677 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
678 if(flr[linpos]<minV)flr[linpos]=minV;
682 float minV=seed[p->total_octave_lines-1];
683 for(;linpos<p->n;linpos++)
684 if(flr[linpos]<minV)flr[linpos]=minV;
689 static void bark_noise_hybridmp(int n,const long *b,
694 long i,hi=b[0]>>16,lo=b[0]>>16,hif=0,lof=0;
707 int ii=(hi<0?-hi:hi);
708 double bin=(f[ii]<-offset?1.:f[ii]+offset);
719 int ii=(lo<0?-lo:lo);
720 double bin=(f[ii]<-offset?1.:f[ii]+offset);
731 if(hif<n && fixed>0){
736 int ii=(hif<0?-hif:hif);
737 double bin=(f[ii]<-offset?1.:f[ii]+offset);
748 int ii=(lof<0?-lof:lof);
749 double bin=(f[ii]<-offset?1.:f[ii]+offset);
764 double denom=1./(na*x2a-xa*xa);
765 double a=(ya*x2a-xya*xa)*denom;
766 double b=(na*xya-xa*ya)*denom;
775 double denomf=1./(nb*x2b-xb*xb);
776 double af=(yb*x2b-xyb*xb)*denomf;
777 double bf=(nb*xyb-xb*yb)*denomf;
781 if(va>vb && vb>0.)va=vb;
791 void _vp_remove_floor(vorbis_look_psy *p,
792 vorbis_look_psy_global *g,
797 float local_specmax){
802 residue[i]=mdct[i]/codedflr[i];
808 void _vp_compute_mask(vorbis_look_psy *p,
809 vorbis_look_psy_global *g,
814 float global_specmax,
817 float bitrate_noise_offset){
821 float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
822 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
825 if(p->vi->noisemaskp){
826 float *work=alloca(n*sizeof(*work));
828 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
831 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
833 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
834 p->vi->noisewindowfixed);
836 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
838 /* work[i] holds the median line (.5), logmask holds the upper
839 envelope line (1.) */
840 _analysis_output("noisemedian",seq,work,n,1,0);
842 for(i=0;i<n;i++)logmask[i]+=work[i];
843 _analysis_output("noiseenvelope",seq,logmask,n,1,0);
844 for(i=0;i<n;i++)logmask[i]-=work[i];
847 int dB=logmask[i]+.5;
848 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
849 logmask[i]= work[i]+p->vi->noisecompand[dB]+p->noiseoffset[i]+bitrate_noise_offset;
850 if(logmask[i]>p->vi->noisemaxsupp)logmask[i]=p->vi->noisemaxsupp;
852 _analysis_output("noise",seq,logmask,n,1,0);
855 for(i=0;i<n;i++)logmask[i]=NEGINF;
858 /* set the ATH (floating below localmax, not global max by a
861 float att=local_specmax+p->vi->ath_adjatt;
862 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
865 float av=p->ath[i]+att;
866 if(av>logmask[i])logmask[i]=av;
871 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
872 max_seeds(p,g,channel,seed,logmask);
874 /* doing this here is clean, but we need to find a faster way to do
875 it than to just tack it on */
877 for(i=0;i<n;i++)if(logmdct[i]>=logmask[i])break;
879 for(i=0;i<n;i++)logmask[i]=NEGINF;
882 logfft[i]=max(logmdct[i],logfft[i]);
888 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
889 vorbis_info *vi=vd->vi;
890 codec_setup_info *ci=vi->codec_setup;
891 vorbis_info_psy_global *gi=&ci->psy_g_param;
893 int n=ci->blocksizes[vd->W]/2;
894 float secs=(float)n/vi->rate;
896 amp+=secs*gi->ampmax_att_per_sec;
897 if(amp<-9999)amp=-9999;
901 static void couple_lossless(float A, float B,
902 float granule,float igranule,
903 float *mag, float *ang,
907 A=rint(A*igranule)*granule; /* must be done *after* the comparison */
908 B=rint(B*igranule)*granule;
910 *mag=A; *ang=(A>0.f?A-B:B-A);
912 A=rint(A*igranule)*granule;
913 B=rint(B*igranule)*granule;
915 *mag=B; *ang=(B>0.f?A-B:B-A);
918 if(flip_p && *ang>fabs(*mag)*1.9999f){
919 *ang= -fabs(*mag)*2.f;
924 static void couple_point(float A, float B, float fA, float fB,
925 float granule,float igranule,
926 float fmag, float *mag, float *ang){
937 corr=sqrt((A*A*fA*fA+B*B*fB*fB)/(fA*fA+fB*fB))/fmag;
938 *mag=rint(*mag*corr*igranule)*granule;
947 void _vp_quantize_couple(vorbis_look_psy *p,
948 vorbis_info_mapping0 *vi,
956 vorbis_info_psy *info=p->vi;
958 /* perform any requested channel coupling */
959 for(i=0;i<vi->coupling_steps;i++){
960 float granulem=info->couple_pass[passno].granulem;
961 float igranulem=info->couple_pass[passno].igranulem;
963 /* make sure coupling a zero and a nonzero channel results in two
965 if(nonzero[vi->coupling_mag[i]] ||
966 nonzero[vi->coupling_ang[i]]){
968 float *pcmM=pcm[vi->coupling_mag[i]];
969 float *pcmA=pcm[vi->coupling_ang[i]];
970 float *floorM=pcm[vi->coupling_mag[i]]+n;
971 float *floorA=pcm[vi->coupling_ang[i]]+n;
972 float *sofarM=sofar[vi->coupling_mag[i]];
973 float *sofarA=sofar[vi->coupling_ang[i]];
974 float *qM=quantized[vi->coupling_mag[i]];
975 float *qA=quantized[vi->coupling_ang[i]];
977 nonzero[vi->coupling_mag[i]]=1;
978 nonzero[vi->coupling_ang[i]]=1;
980 for(j=0,k=0;j<n;k++){
981 vp_couple *part=info->couple_pass[passno].couple_pass+k;
982 float rqlimit=part->outofphase_requant_limit;
983 float flip_p=part->outofphase_redundant_flip_p;
985 for(;j<part->limit && j<p->n;j++){
986 /* partition by partition; k is our by-location partition
988 float ang,mag,fmag=max(fabs(pcmM[j]),fabs(pcmA[j]));
990 if(fmag<part->amppost_point){
991 couple_point(pcmM[j],pcmA[j],floorM[j],floorA[j],
992 granulem,igranulem,fmag,&mag,&ang);
995 couple_lossless(pcmM[j],pcmA[j],
996 granulem,igranulem,&mag,&ang,flip_p);
999 /* executive decision time: when requantizing and recoupling
1000 residue in order to progressively encode at finer
1001 resolution, an out of phase component that originally
1002 quntized to 2*mag can flip flop magnitude/angle if it
1003 requantizes to not-quite out of phase. If that happens,
1004 we opt not to fill in additional resolution (in order to
1005 simplify the iterative codebook design and
1008 if(ang<-rqlimit || ang>rqlimit){
1012 qM[j]=mag-sofarM[j];
1013 qA[j]=ang-sofarA[j];