Line | Count | Source (jump to first uncovered line) |
1 | | /******************************************************************** |
2 | | * * |
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. * |
7 | | * * |
8 | | * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2010 * |
9 | | * by the Xiph.Org Foundation https://xiph.org/ * |
10 | | * * |
11 | | ******************************************************************** |
12 | | |
13 | | function: psychoacoustics not including preecho |
14 | | |
15 | | ********************************************************************/ |
16 | | |
17 | | #include <stdlib.h> |
18 | | #include <math.h> |
19 | | #include <string.h> |
20 | | #include "vorbis/codec.h" |
21 | | #include "codec_internal.h" |
22 | | |
23 | | #include "masking.h" |
24 | | #include "psy.h" |
25 | | #include "os.h" |
26 | | #include "lpc.h" |
27 | | #include "smallft.h" |
28 | | #include "scales.h" |
29 | | #include "misc.h" |
30 | | |
31 | 0 | #define NEGINF -9999.f |
32 | | static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10}; |
33 | | static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10}; |
34 | | |
35 | 0 | vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){ |
36 | 0 | codec_setup_info *ci=vi->codec_setup; |
37 | 0 | vorbis_info_psy_global *gi=&ci->psy_g_param; |
38 | 0 | vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look)); |
39 | |
|
40 | 0 | look->channels=vi->channels; |
41 | |
|
42 | 0 | look->ampmax=-9999.; |
43 | 0 | look->gi=gi; |
44 | 0 | return(look); |
45 | 0 | } |
46 | | |
47 | 0 | void _vp_global_free(vorbis_look_psy_global *look){ |
48 | 0 | if(look){ |
49 | 0 | memset(look,0,sizeof(*look)); |
50 | 0 | _ogg_free(look); |
51 | 0 | } |
52 | 0 | } |
53 | | |
54 | 0 | void _vi_gpsy_free(vorbis_info_psy_global *i){ |
55 | 0 | if(i){ |
56 | 0 | memset(i,0,sizeof(*i)); |
57 | 0 | _ogg_free(i); |
58 | 0 | } |
59 | 0 | } |
60 | | |
61 | 0 | void _vi_psy_free(vorbis_info_psy *i){ |
62 | 0 | if(i){ |
63 | 0 | memset(i,0,sizeof(*i)); |
64 | 0 | _ogg_free(i); |
65 | 0 | } |
66 | 0 | } |
67 | | |
68 | | static void min_curve(float *c, |
69 | 0 | float *c2){ |
70 | 0 | int i; |
71 | 0 | for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i]; |
72 | 0 | } |
73 | | static void max_curve(float *c, |
74 | 0 | float *c2){ |
75 | 0 | int i; |
76 | 0 | for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i]; |
77 | 0 | } |
78 | | |
79 | 0 | static void attenuate_curve(float *c,float att){ |
80 | 0 | int i; |
81 | 0 | for(i=0;i<EHMER_MAX;i++) |
82 | 0 | c[i]+=att; |
83 | 0 | } |
84 | | |
85 | | static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n, |
86 | 0 | float center_boost, float center_decay_rate){ |
87 | 0 | int i,j,k,m; |
88 | 0 | float ath[EHMER_MAX]; |
89 | 0 | float workc[P_BANDS][P_LEVELS][EHMER_MAX]; |
90 | 0 | float athc[P_LEVELS][EHMER_MAX]; |
91 | 0 | float *brute_buffer=alloca(n*sizeof(*brute_buffer)); |
92 | |
|
93 | 0 | float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS); |
94 | |
|
95 | 0 | memset(workc,0,sizeof(workc)); |
96 | |
|
97 | 0 | for(i=0;i<P_BANDS;i++){ |
98 | | /* we add back in the ATH to avoid low level curves falling off to |
99 | | -infinity and unnecessarily cutting off high level curves in the |
100 | | curve limiting (last step). */ |
101 | | |
102 | | /* A half-band's settings must be valid over the whole band, and |
103 | | it's better to mask too little than too much */ |
104 | 0 | int ath_offset=i*4; |
105 | 0 | for(j=0;j<EHMER_MAX;j++){ |
106 | 0 | float min=999.; |
107 | 0 | for(k=0;k<4;k++) |
108 | 0 | if(j+k+ath_offset<MAX_ATH){ |
109 | 0 | if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset]; |
110 | 0 | }else{ |
111 | 0 | if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1]; |
112 | 0 | } |
113 | 0 | ath[j]=min; |
114 | 0 | } |
115 | | |
116 | | /* copy curves into working space, replicate the 50dB curve to 30 |
117 | | and 40, replicate the 100dB curve to 110 */ |
118 | 0 | for(j=0;j<6;j++) |
119 | 0 | memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j])); |
120 | 0 | memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); |
121 | 0 | memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); |
122 | | |
123 | | /* apply centered curve boost/decay */ |
124 | 0 | for(j=0;j<P_LEVELS;j++){ |
125 | 0 | for(k=0;k<EHMER_MAX;k++){ |
126 | 0 | float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate; |
127 | 0 | if(adj<0. && center_boost>0)adj=0.; |
128 | 0 | if(adj>0. && center_boost<0)adj=0.; |
129 | 0 | workc[i][j][k]+=adj; |
130 | 0 | } |
131 | 0 | } |
132 | | |
133 | | /* normalize curves so the driving amplitude is 0dB */ |
134 | | /* make temp curves with the ATH overlayed */ |
135 | 0 | for(j=0;j<P_LEVELS;j++){ |
136 | 0 | attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0); |
137 | 0 | memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc)); |
138 | 0 | attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0); |
139 | 0 | max_curve(athc[j],workc[i][j]); |
140 | 0 | } |
141 | | |
142 | | /* Now limit the louder curves. |
143 | | |
144 | | the idea is this: We don't know what the playback attenuation |
145 | | will be; 0dB SL moves every time the user twiddles the volume |
146 | | knob. So that means we have to use a single 'most pessimal' curve |
147 | | for all masking amplitudes, right? Wrong. The *loudest* sound |
148 | | can be in (we assume) a range of ...+100dB] SL. However, sounds |
149 | | 20dB down will be in a range ...+80], 40dB down is from ...+60], |
150 | | etc... */ |
151 | |
|
152 | 0 | for(j=1;j<P_LEVELS;j++){ |
153 | 0 | min_curve(athc[j],athc[j-1]); |
154 | 0 | min_curve(workc[i][j],athc[j]); |
155 | 0 | } |
156 | 0 | } |
157 | |
|
158 | 0 | for(i=0;i<P_BANDS;i++){ |
159 | 0 | int hi_curve,lo_curve,bin; |
160 | 0 | ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS); |
161 | | |
162 | | /* low frequency curves are measured with greater resolution than |
163 | | the MDCT/FFT will actually give us; we want the curve applied |
164 | | to the tone data to be pessimistic and thus apply the minimum |
165 | | masking possible for a given bin. That means that a single bin |
166 | | could span more than one octave and that the curve will be a |
167 | | composite of multiple octaves. It also may mean that a single |
168 | | bin may span > an eighth of an octave and that the eighth |
169 | | octave values may also be composited. */ |
170 | | |
171 | | /* which octave curves will we be compositing? */ |
172 | 0 | bin=floor(fromOC(i*.5)/binHz); |
173 | 0 | lo_curve= ceil(toOC(bin*binHz+1)*2); |
174 | 0 | hi_curve= floor(toOC((bin+1)*binHz)*2); |
175 | 0 | if(lo_curve>i)lo_curve=i; |
176 | 0 | if(lo_curve<0)lo_curve=0; |
177 | 0 | if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1; |
178 | |
|
179 | 0 | for(m=0;m<P_LEVELS;m++){ |
180 | 0 | ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2)); |
181 | |
|
182 | 0 | for(j=0;j<n;j++)brute_buffer[j]=999.; |
183 | | |
184 | | /* render the curve into bins, then pull values back into curve. |
185 | | The point is that any inherent subsampling aliasing results in |
186 | | a safe minimum */ |
187 | 0 | for(k=lo_curve;k<=hi_curve;k++){ |
188 | 0 | int l=0; |
189 | |
|
190 | 0 | for(j=0;j<EHMER_MAX;j++){ |
191 | 0 | int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz; |
192 | 0 | int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1; |
193 | |
|
194 | 0 | if(lo_bin<0)lo_bin=0; |
195 | 0 | if(lo_bin>n)lo_bin=n; |
196 | 0 | if(lo_bin<l)l=lo_bin; |
197 | 0 | if(hi_bin<0)hi_bin=0; |
198 | 0 | if(hi_bin>n)hi_bin=n; |
199 | |
|
200 | 0 | for(;l<hi_bin && l<n;l++) |
201 | 0 | if(brute_buffer[l]>workc[k][m][j]) |
202 | 0 | brute_buffer[l]=workc[k][m][j]; |
203 | 0 | } |
204 | |
|
205 | 0 | for(;l<n;l++) |
206 | 0 | if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) |
207 | 0 | brute_buffer[l]=workc[k][m][EHMER_MAX-1]; |
208 | |
|
209 | 0 | } |
210 | | |
211 | | /* be equally paranoid about being valid up to next half ocatve */ |
212 | 0 | if(i+1<P_BANDS){ |
213 | 0 | int l=0; |
214 | 0 | k=i+1; |
215 | 0 | for(j=0;j<EHMER_MAX;j++){ |
216 | 0 | int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz; |
217 | 0 | int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1; |
218 | |
|
219 | 0 | if(lo_bin<0)lo_bin=0; |
220 | 0 | if(lo_bin>n)lo_bin=n; |
221 | 0 | if(lo_bin<l)l=lo_bin; |
222 | 0 | if(hi_bin<0)hi_bin=0; |
223 | 0 | if(hi_bin>n)hi_bin=n; |
224 | |
|
225 | 0 | for(;l<hi_bin && l<n;l++) |
226 | 0 | if(brute_buffer[l]>workc[k][m][j]) |
227 | 0 | brute_buffer[l]=workc[k][m][j]; |
228 | 0 | } |
229 | |
|
230 | 0 | for(;l<n;l++) |
231 | 0 | if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) |
232 | 0 | brute_buffer[l]=workc[k][m][EHMER_MAX-1]; |
233 | |
|
234 | 0 | } |
235 | | |
236 | |
|
237 | 0 | for(j=0;j<EHMER_MAX;j++){ |
238 | 0 | int bin=fromOC(j*.125+i*.5-2.)/binHz; |
239 | 0 | if(bin<0){ |
240 | 0 | ret[i][m][j+2]=-999.; |
241 | 0 | }else{ |
242 | 0 | if(bin>=n){ |
243 | 0 | ret[i][m][j+2]=-999.; |
244 | 0 | }else{ |
245 | 0 | ret[i][m][j+2]=brute_buffer[bin]; |
246 | 0 | } |
247 | 0 | } |
248 | 0 | } |
249 | | |
250 | | /* add fenceposts */ |
251 | 0 | for(j=0;j<EHMER_OFFSET;j++) |
252 | 0 | if(ret[i][m][j+2]>-200.f)break; |
253 | 0 | ret[i][m][0]=j; |
254 | |
|
255 | 0 | for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--) |
256 | 0 | if(ret[i][m][j+2]>-200.f) |
257 | 0 | break; |
258 | 0 | ret[i][m][1]=j; |
259 | |
|
260 | 0 | } |
261 | 0 | } |
262 | |
|
263 | 0 | return(ret); |
264 | 0 | } |
265 | | |
266 | | void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi, |
267 | 0 | vorbis_info_psy_global *gi,int n,long rate){ |
268 | 0 | long i,j,lo=-99,hi=1; |
269 | 0 | long maxoc; |
270 | 0 | memset(p,0,sizeof(*p)); |
271 | |
|
272 | 0 | p->eighth_octave_lines=gi->eighth_octave_lines; |
273 | 0 | p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1; |
274 | |
|
275 | 0 | p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines; |
276 | 0 | maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f; |
277 | 0 | p->total_octave_lines=maxoc-p->firstoc+1; |
278 | 0 | p->ath=_ogg_malloc(n*sizeof(*p->ath)); |
279 | |
|
280 | 0 | p->octave=_ogg_malloc(n*sizeof(*p->octave)); |
281 | 0 | p->bark=_ogg_malloc(n*sizeof(*p->bark)); |
282 | 0 | p->vi=vi; |
283 | 0 | p->n=n; |
284 | 0 | p->rate=rate; |
285 | | |
286 | | /* AoTuV HF weighting */ |
287 | 0 | p->m_val = 1.; |
288 | 0 | if(rate < 26000) p->m_val = 0; |
289 | 0 | else if(rate < 38000) p->m_val = .94; /* 32kHz */ |
290 | 0 | else if(rate > 46000) p->m_val = 1.275; /* 48kHz */ |
291 | | |
292 | | /* set up the lookups for a given blocksize and sample rate */ |
293 | |
|
294 | 0 | for(i=0,j=0;i<MAX_ATH-1;i++){ |
295 | 0 | int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate); |
296 | 0 | float base=ATH[i]; |
297 | 0 | if(j<endpos){ |
298 | 0 | float delta=(ATH[i+1]-base)/(endpos-j); |
299 | 0 | for(;j<endpos && j<n;j++){ |
300 | 0 | p->ath[j]=base+100.; |
301 | 0 | base+=delta; |
302 | 0 | } |
303 | 0 | } |
304 | 0 | } |
305 | |
|
306 | 0 | for(;j<n;j++){ |
307 | 0 | p->ath[j]=p->ath[j-1]; |
308 | 0 | } |
309 | |
|
310 | 0 | for(i=0;i<n;i++){ |
311 | 0 | float bark=toBARK(rate/(2*n)*i); |
312 | |
|
313 | 0 | for(;lo+vi->noisewindowlomin<i && |
314 | 0 | toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++); |
315 | |
|
316 | 0 | for(;hi<=n && (hi<i+vi->noisewindowhimin || |
317 | 0 | toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++); |
318 | |
|
319 | 0 | p->bark[i]=((lo-1)<<16)+(hi-1); |
320 | |
|
321 | 0 | } |
322 | |
|
323 | 0 | for(i=0;i<n;i++) |
324 | 0 | p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f; |
325 | |
|
326 | 0 | p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n, |
327 | 0 | vi->tone_centerboost,vi->tone_decay); |
328 | | |
329 | | /* set up rolling noise median */ |
330 | 0 | p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset)); |
331 | 0 | for(i=0;i<P_NOISECURVES;i++) |
332 | 0 | p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset)); |
333 | |
|
334 | 0 | for(i=0;i<n;i++){ |
335 | 0 | float halfoc=toOC((i+.5)*rate/(2.*n))*2.; |
336 | 0 | int inthalfoc; |
337 | 0 | float del; |
338 | |
|
339 | 0 | if(halfoc<0)halfoc=0; |
340 | 0 | if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1; |
341 | 0 | inthalfoc=(int)halfoc; |
342 | 0 | del=halfoc-inthalfoc; |
343 | |
|
344 | 0 | for(j=0;j<P_NOISECURVES;j++) |
345 | 0 | p->noiseoffset[j][i]= |
346 | 0 | p->vi->noiseoff[j][inthalfoc]*(1.-del) + |
347 | 0 | p->vi->noiseoff[j][inthalfoc+1]*del; |
348 | |
|
349 | 0 | } |
350 | | #if 0 |
351 | | { |
352 | | static int ls=0; |
353 | | _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0); |
354 | | _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0); |
355 | | _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0); |
356 | | } |
357 | | #endif |
358 | 0 | } |
359 | | |
360 | 0 | void _vp_psy_clear(vorbis_look_psy *p){ |
361 | 0 | int i,j; |
362 | 0 | if(p){ |
363 | 0 | if(p->ath)_ogg_free(p->ath); |
364 | 0 | if(p->octave)_ogg_free(p->octave); |
365 | 0 | if(p->bark)_ogg_free(p->bark); |
366 | 0 | if(p->tonecurves){ |
367 | 0 | for(i=0;i<P_BANDS;i++){ |
368 | 0 | for(j=0;j<P_LEVELS;j++){ |
369 | 0 | _ogg_free(p->tonecurves[i][j]); |
370 | 0 | } |
371 | 0 | _ogg_free(p->tonecurves[i]); |
372 | 0 | } |
373 | 0 | _ogg_free(p->tonecurves); |
374 | 0 | } |
375 | 0 | if(p->noiseoffset){ |
376 | 0 | for(i=0;i<P_NOISECURVES;i++){ |
377 | 0 | _ogg_free(p->noiseoffset[i]); |
378 | 0 | } |
379 | 0 | _ogg_free(p->noiseoffset); |
380 | 0 | } |
381 | 0 | memset(p,0,sizeof(*p)); |
382 | 0 | } |
383 | 0 | } |
384 | | |
385 | | /* octave/(8*eighth_octave_lines) x scale and dB y scale */ |
386 | | static void seed_curve(float *seed, |
387 | | const float **curves, |
388 | | float amp, |
389 | | int oc, int n, |
390 | 0 | int linesper,float dBoffset){ |
391 | 0 | int i,post1; |
392 | 0 | int seedptr; |
393 | 0 | const float *posts,*curve; |
394 | |
|
395 | 0 | int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f); |
396 | 0 | choice=max(choice,0); |
397 | 0 | choice=min(choice,P_LEVELS-1); |
398 | 0 | posts=curves[choice]; |
399 | 0 | curve=posts+2; |
400 | 0 | post1=(int)posts[1]; |
401 | 0 | seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1); |
402 | |
|
403 | 0 | for(i=posts[0];i<post1;i++){ |
404 | 0 | if(seedptr>0){ |
405 | 0 | float lin=amp+curve[i]; |
406 | 0 | if(seed[seedptr]<lin)seed[seedptr]=lin; |
407 | 0 | } |
408 | 0 | seedptr+=linesper; |
409 | 0 | if(seedptr>=n)break; |
410 | 0 | } |
411 | 0 | } |
412 | | |
413 | | static void seed_loop(vorbis_look_psy *p, |
414 | | const float ***curves, |
415 | | const float *f, |
416 | | const float *flr, |
417 | | float *seed, |
418 | 0 | float specmax){ |
419 | 0 | vorbis_info_psy *vi=p->vi; |
420 | 0 | long n=p->n,i; |
421 | 0 | float dBoffset=vi->max_curve_dB-specmax; |
422 | | |
423 | | /* prime the working vector with peak values */ |
424 | |
|
425 | 0 | for(i=0;i<n;i++){ |
426 | 0 | float max=f[i]; |
427 | 0 | long oc=p->octave[i]; |
428 | 0 | while(i+1<n && p->octave[i+1]==oc){ |
429 | 0 | i++; |
430 | 0 | if(f[i]>max)max=f[i]; |
431 | 0 | } |
432 | |
|
433 | 0 | if(max+6.f>flr[i]){ |
434 | 0 | oc=oc>>p->shiftoc; |
435 | |
|
436 | 0 | if(oc>=P_BANDS)oc=P_BANDS-1; |
437 | 0 | if(oc<0)oc=0; |
438 | |
|
439 | 0 | seed_curve(seed, |
440 | 0 | curves[oc], |
441 | 0 | max, |
442 | 0 | p->octave[i]-p->firstoc, |
443 | 0 | p->total_octave_lines, |
444 | 0 | p->eighth_octave_lines, |
445 | 0 | dBoffset); |
446 | 0 | } |
447 | 0 | } |
448 | 0 | } |
449 | | |
450 | 0 | static void seed_chase(float *seeds, int linesper, long n){ |
451 | 0 | long *posstack=alloca(n*sizeof(*posstack)); |
452 | 0 | float *ampstack=alloca(n*sizeof(*ampstack)); |
453 | 0 | long stack=0; |
454 | 0 | long pos=0; |
455 | 0 | long i; |
456 | |
|
457 | 0 | for(i=0;i<n;i++){ |
458 | 0 | if(stack<2){ |
459 | 0 | posstack[stack]=i; |
460 | 0 | ampstack[stack++]=seeds[i]; |
461 | 0 | }else{ |
462 | 0 | while(1){ |
463 | 0 | if(seeds[i]<ampstack[stack-1]){ |
464 | 0 | posstack[stack]=i; |
465 | 0 | ampstack[stack++]=seeds[i]; |
466 | 0 | break; |
467 | 0 | }else{ |
468 | 0 | if(i<posstack[stack-1]+linesper){ |
469 | 0 | if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] && |
470 | 0 | i<posstack[stack-2]+linesper){ |
471 | | /* we completely overlap, making stack-1 irrelevant. pop it */ |
472 | 0 | stack--; |
473 | 0 | continue; |
474 | 0 | } |
475 | 0 | } |
476 | 0 | posstack[stack]=i; |
477 | 0 | ampstack[stack++]=seeds[i]; |
478 | 0 | break; |
479 | |
|
480 | 0 | } |
481 | 0 | } |
482 | 0 | } |
483 | 0 | } |
484 | | |
485 | | /* the stack now contains only the positions that are relevant. Scan |
486 | | 'em straight through */ |
487 | |
|
488 | 0 | for(i=0;i<stack;i++){ |
489 | 0 | long endpos; |
490 | 0 | if(i<stack-1 && ampstack[i+1]>ampstack[i]){ |
491 | 0 | endpos=posstack[i+1]; |
492 | 0 | }else{ |
493 | 0 | endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is |
494 | | discarded in short frames */ |
495 | 0 | } |
496 | 0 | if(endpos>n)endpos=n; |
497 | 0 | for(;pos<endpos;pos++) |
498 | 0 | seeds[pos]=ampstack[i]; |
499 | 0 | } |
500 | | |
501 | | /* there. Linear time. I now remember this was on a problem set I |
502 | | had in Grad Skool... I didn't solve it at the time ;-) */ |
503 | |
|
504 | 0 | } |
505 | | |
506 | | /* bleaugh, this is more complicated than it needs to be */ |
507 | | #include<stdio.h> |
508 | | static void max_seeds(vorbis_look_psy *p, |
509 | | float *seed, |
510 | 0 | float *flr){ |
511 | 0 | long n=p->total_octave_lines; |
512 | 0 | int linesper=p->eighth_octave_lines; |
513 | 0 | long linpos=0; |
514 | 0 | long pos; |
515 | |
|
516 | 0 | seed_chase(seed,linesper,n); /* for masking */ |
517 | |
|
518 | 0 | pos=p->octave[0]-p->firstoc-(linesper>>1); |
519 | |
|
520 | 0 | while(linpos+1<p->n){ |
521 | 0 | float minV=seed[pos]; |
522 | 0 | long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc; |
523 | 0 | if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit; |
524 | 0 | while(pos+1<=end){ |
525 | 0 | pos++; |
526 | 0 | if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF) |
527 | 0 | minV=seed[pos]; |
528 | 0 | } |
529 | |
|
530 | 0 | end=pos+p->firstoc; |
531 | 0 | for(;linpos<p->n && p->octave[linpos]<=end;linpos++) |
532 | 0 | if(flr[linpos]<minV)flr[linpos]=minV; |
533 | 0 | } |
534 | |
|
535 | 0 | { |
536 | 0 | float minV=seed[p->total_octave_lines-1]; |
537 | 0 | for(;linpos<p->n;linpos++) |
538 | 0 | if(flr[linpos]<minV)flr[linpos]=minV; |
539 | 0 | } |
540 | |
|
541 | 0 | } |
542 | | |
543 | | static void bark_noise_hybridmp(int n,const long *b, |
544 | | const float *f, |
545 | | float *noise, |
546 | | const float offset, |
547 | 0 | const int fixed){ |
548 | |
|
549 | 0 | float *N=alloca(n*sizeof(*N)); |
550 | 0 | float *X=alloca(n*sizeof(*N)); |
551 | 0 | float *XX=alloca(n*sizeof(*N)); |
552 | 0 | float *Y=alloca(n*sizeof(*N)); |
553 | 0 | float *XY=alloca(n*sizeof(*N)); |
554 | |
|
555 | 0 | float tN, tX, tXX, tY, tXY; |
556 | 0 | int i; |
557 | |
|
558 | 0 | int lo, hi; |
559 | 0 | float R=0.f; |
560 | 0 | float A=0.f; |
561 | 0 | float B=0.f; |
562 | 0 | float D=1.f; |
563 | 0 | float w, x, y; |
564 | |
|
565 | 0 | tN = tX = tXX = tY = tXY = 0.f; |
566 | |
|
567 | 0 | y = f[0] + offset; |
568 | 0 | if (y < 1.f) y = 1.f; |
569 | |
|
570 | 0 | w = y * y * .5; |
571 | |
|
572 | 0 | tN += w; |
573 | 0 | tX += w; |
574 | 0 | tY += w * y; |
575 | |
|
576 | 0 | N[0] = tN; |
577 | 0 | X[0] = tX; |
578 | 0 | XX[0] = tXX; |
579 | 0 | Y[0] = tY; |
580 | 0 | XY[0] = tXY; |
581 | |
|
582 | 0 | for (i = 1, x = 1.f; i < n; i++, x += 1.f) { |
583 | |
|
584 | 0 | y = f[i] + offset; |
585 | 0 | if (y < 1.f) y = 1.f; |
586 | |
|
587 | 0 | w = y * y; |
588 | |
|
589 | 0 | tN += w; |
590 | 0 | tX += w * x; |
591 | 0 | tXX += w * x * x; |
592 | 0 | tY += w * y; |
593 | 0 | tXY += w * x * y; |
594 | |
|
595 | 0 | N[i] = tN; |
596 | 0 | X[i] = tX; |
597 | 0 | XX[i] = tXX; |
598 | 0 | Y[i] = tY; |
599 | 0 | XY[i] = tXY; |
600 | 0 | } |
601 | |
|
602 | 0 | for (i = 0, x = 0.f; i < n; i++, x += 1.f) { |
603 | |
|
604 | 0 | lo = b[i] >> 16; |
605 | 0 | hi = b[i] & 0xffff; |
606 | 0 | if( lo>=0 || -lo>=n ) break; |
607 | 0 | if( hi>=n ) break; |
608 | | |
609 | 0 | tN = N[hi] + N[-lo]; |
610 | 0 | tX = X[hi] - X[-lo]; |
611 | 0 | tXX = XX[hi] + XX[-lo]; |
612 | 0 | tY = Y[hi] + Y[-lo]; |
613 | 0 | tXY = XY[hi] - XY[-lo]; |
614 | |
|
615 | 0 | A = tY * tXX - tX * tXY; |
616 | 0 | B = tN * tXY - tX * tY; |
617 | 0 | D = tN * tXX - tX * tX; |
618 | 0 | R = (A + x * B) / D; |
619 | 0 | if (R < 0.f) R = 0.f; |
620 | |
|
621 | 0 | noise[i] = R - offset; |
622 | 0 | } |
623 | |
|
624 | 0 | for ( ; i < n; i++, x += 1.f) { |
625 | |
|
626 | 0 | lo = b[i] >> 16; |
627 | 0 | hi = b[i] & 0xffff; |
628 | 0 | if( lo<0 || lo>=n ) break; |
629 | 0 | if( hi>=n ) break; |
630 | | |
631 | 0 | tN = N[hi] - N[lo]; |
632 | 0 | tX = X[hi] - X[lo]; |
633 | 0 | tXX = XX[hi] - XX[lo]; |
634 | 0 | tY = Y[hi] - Y[lo]; |
635 | 0 | tXY = XY[hi] - XY[lo]; |
636 | |
|
637 | 0 | A = tY * tXX - tX * tXY; |
638 | 0 | B = tN * tXY - tX * tY; |
639 | 0 | D = tN * tXX - tX * tX; |
640 | 0 | R = (A + x * B) / D; |
641 | 0 | if (R < 0.f) R = 0.f; |
642 | |
|
643 | 0 | noise[i] = R - offset; |
644 | 0 | } |
645 | |
|
646 | 0 | for ( ; i < n; i++, x += 1.f) { |
647 | |
|
648 | 0 | R = (A + x * B) / D; |
649 | 0 | if (R < 0.f) R = 0.f; |
650 | |
|
651 | 0 | noise[i] = R - offset; |
652 | 0 | } |
653 | |
|
654 | 0 | if (fixed <= 0) return; |
655 | | |
656 | 0 | for (i = 0, x = 0.f; i < n; i++, x += 1.f) { |
657 | 0 | hi = i + fixed / 2; |
658 | 0 | lo = hi - fixed; |
659 | 0 | if ( hi>=n ) break; |
660 | 0 | if ( lo>=0 ) break; |
661 | | |
662 | 0 | tN = N[hi] + N[-lo]; |
663 | 0 | tX = X[hi] - X[-lo]; |
664 | 0 | tXX = XX[hi] + XX[-lo]; |
665 | 0 | tY = Y[hi] + Y[-lo]; |
666 | 0 | tXY = XY[hi] - XY[-lo]; |
667 | | |
668 | |
|
669 | 0 | A = tY * tXX - tX * tXY; |
670 | 0 | B = tN * tXY - tX * tY; |
671 | 0 | D = tN * tXX - tX * tX; |
672 | 0 | R = (A + x * B) / D; |
673 | |
|
674 | 0 | if (R - offset < noise[i]) noise[i] = R - offset; |
675 | 0 | } |
676 | 0 | for ( ; i < n; i++, x += 1.f) { |
677 | |
|
678 | 0 | hi = i + fixed / 2; |
679 | 0 | lo = hi - fixed; |
680 | 0 | if ( hi>=n ) break; |
681 | 0 | if ( lo<0 ) break; |
682 | | |
683 | 0 | tN = N[hi] - N[lo]; |
684 | 0 | tX = X[hi] - X[lo]; |
685 | 0 | tXX = XX[hi] - XX[lo]; |
686 | 0 | tY = Y[hi] - Y[lo]; |
687 | 0 | tXY = XY[hi] - XY[lo]; |
688 | |
|
689 | 0 | A = tY * tXX - tX * tXY; |
690 | 0 | B = tN * tXY - tX * tY; |
691 | 0 | D = tN * tXX - tX * tX; |
692 | 0 | R = (A + x * B) / D; |
693 | |
|
694 | 0 | if (R - offset < noise[i]) noise[i] = R - offset; |
695 | 0 | } |
696 | 0 | for ( ; i < n; i++, x += 1.f) { |
697 | 0 | R = (A + x * B) / D; |
698 | 0 | if (R - offset < noise[i]) noise[i] = R - offset; |
699 | 0 | } |
700 | 0 | } |
701 | | |
702 | | void _vp_noisemask(vorbis_look_psy *p, |
703 | | float *logmdct, |
704 | 0 | float *logmask){ |
705 | |
|
706 | 0 | int i,n=p->n; |
707 | 0 | float *work=alloca(n*sizeof(*work)); |
708 | |
|
709 | 0 | bark_noise_hybridmp(n,p->bark,logmdct,logmask, |
710 | 0 | 140.,-1); |
711 | |
|
712 | 0 | for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i]; |
713 | |
|
714 | 0 | bark_noise_hybridmp(n,p->bark,work,logmask,0., |
715 | 0 | p->vi->noisewindowfixed); |
716 | |
|
717 | 0 | for(i=0;i<n;i++)work[i]=logmdct[i]-work[i]; |
718 | |
|
719 | | #if 0 |
720 | | { |
721 | | static int seq=0; |
722 | | |
723 | | float work2[n]; |
724 | | for(i=0;i<n;i++){ |
725 | | work2[i]=logmask[i]+work[i]; |
726 | | } |
727 | | |
728 | | if(seq&1) |
729 | | _analysis_output("median2R",seq/2,work,n,1,0,0); |
730 | | else |
731 | | _analysis_output("median2L",seq/2,work,n,1,0,0); |
732 | | |
733 | | if(seq&1) |
734 | | _analysis_output("envelope2R",seq/2,work2,n,1,0,0); |
735 | | else |
736 | | _analysis_output("envelope2L",seq/2,work2,n,1,0,0); |
737 | | seq++; |
738 | | } |
739 | | #endif |
740 | |
|
741 | 0 | for(i=0;i<n;i++){ |
742 | 0 | int dB=logmask[i]+.5; |
743 | 0 | if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1; |
744 | 0 | if(dB<0)dB=0; |
745 | 0 | logmask[i]= work[i]+p->vi->noisecompand[dB]; |
746 | 0 | } |
747 | |
|
748 | 0 | } |
749 | | |
750 | | void _vp_tonemask(vorbis_look_psy *p, |
751 | | float *logfft, |
752 | | float *logmask, |
753 | | float global_specmax, |
754 | 0 | float local_specmax){ |
755 | |
|
756 | 0 | int i,n=p->n; |
757 | |
|
758 | 0 | float *seed=alloca(sizeof(*seed)*p->total_octave_lines); |
759 | 0 | float att=local_specmax+p->vi->ath_adjatt; |
760 | 0 | for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF; |
761 | | |
762 | | /* set the ATH (floating below localmax, not global max by a |
763 | | specified att) */ |
764 | 0 | if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt; |
765 | |
|
766 | 0 | for(i=0;i<n;i++) |
767 | 0 | logmask[i]=p->ath[i]+att; |
768 | | |
769 | | /* tone masking */ |
770 | 0 | seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax); |
771 | 0 | max_seeds(p,seed,logmask); |
772 | |
|
773 | 0 | } |
774 | | |
775 | | void _vp_offset_and_mix(vorbis_look_psy *p, |
776 | | float *noise, |
777 | | float *tone, |
778 | | int offset_select, |
779 | | float *logmask, |
780 | | float *mdct, |
781 | 0 | float *logmdct){ |
782 | 0 | int i,n=p->n; |
783 | 0 | float de, coeffi, cx;/* AoTuV */ |
784 | 0 | float toneatt=p->vi->tone_masteratt[offset_select]; |
785 | |
|
786 | 0 | cx = p->m_val; |
787 | |
|
788 | 0 | for(i=0;i<n;i++){ |
789 | 0 | float val= noise[i]+p->noiseoffset[offset_select][i]; |
790 | 0 | if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp; |
791 | 0 | logmask[i]=max(val,tone[i]+toneatt); |
792 | | |
793 | | |
794 | | /* AoTuV */ |
795 | | /** @ M1 ** |
796 | | The following codes improve a noise problem. |
797 | | A fundamental idea uses the value of masking and carries out |
798 | | the relative compensation of the MDCT. |
799 | | However, this code is not perfect and all noise problems cannot be solved. |
800 | | by Aoyumi @ 2004/04/18 |
801 | | */ |
802 | |
|
803 | 0 | if(offset_select == 1) { |
804 | 0 | coeffi = -17.2; /* coeffi is a -17.2dB threshold */ |
805 | 0 | val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */ |
806 | |
|
807 | 0 | if(val > coeffi){ |
808 | | /* mdct value is > -17.2 dB below floor */ |
809 | |
|
810 | 0 | de = 1.0-((val-coeffi)*0.005*cx); |
811 | | /* pro-rated attenuation: |
812 | | -0.00 dB boost if mdct value is -17.2dB (relative to floor) |
813 | | -0.77 dB boost if mdct value is 0dB (relative to floor) |
814 | | -1.64 dB boost if mdct value is +17.2dB (relative to floor) |
815 | | etc... */ |
816 | |
|
817 | 0 | if(de < 0) de = 0.0001; |
818 | 0 | }else |
819 | | /* mdct value is <= -17.2 dB below floor */ |
820 | | |
821 | 0 | de = 1.0-((val-coeffi)*0.0003*cx); |
822 | | /* pro-rated attenuation: |
823 | | +0.00 dB atten if mdct value is -17.2dB (relative to floor) |
824 | | +0.45 dB atten if mdct value is -34.4dB (relative to floor) |
825 | | etc... */ |
826 | |
|
827 | 0 | mdct[i] *= de; |
828 | |
|
829 | 0 | } |
830 | 0 | } |
831 | 0 | } |
832 | | |
833 | 0 | float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){ |
834 | 0 | vorbis_info *vi=vd->vi; |
835 | 0 | codec_setup_info *ci=vi->codec_setup; |
836 | 0 | vorbis_info_psy_global *gi=&ci->psy_g_param; |
837 | |
|
838 | 0 | int n=ci->blocksizes[vd->W]/2; |
839 | 0 | float secs=(float)n/vi->rate; |
840 | |
|
841 | 0 | amp+=secs*gi->ampmax_att_per_sec; |
842 | 0 | if(amp<-9999)amp=-9999; |
843 | 0 | return(amp); |
844 | 0 | } |
845 | | |
846 | | static float FLOOR1_fromdB_LOOKUP[256]={ |
847 | | 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F, |
848 | | 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F, |
849 | | 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F, |
850 | | 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F, |
851 | | 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F, |
852 | | 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F, |
853 | | 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F, |
854 | | 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F, |
855 | | 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F, |
856 | | 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F, |
857 | | 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F, |
858 | | 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F, |
859 | | 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F, |
860 | | 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F, |
861 | | 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F, |
862 | | 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F, |
863 | | 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F, |
864 | | 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F, |
865 | | 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F, |
866 | | 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F, |
867 | | 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F, |
868 | | 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F, |
869 | | 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F, |
870 | | 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F, |
871 | | 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F, |
872 | | 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F, |
873 | | 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F, |
874 | | 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F, |
875 | | 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F, |
876 | | 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F, |
877 | | 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F, |
878 | | 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F, |
879 | | 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F, |
880 | | 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F, |
881 | | 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F, |
882 | | 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F, |
883 | | 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F, |
884 | | 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F, |
885 | | 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F, |
886 | | 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F, |
887 | | 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F, |
888 | | 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F, |
889 | | 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F, |
890 | | 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F, |
891 | | 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F, |
892 | | 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F, |
893 | | 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F, |
894 | | 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F, |
895 | | 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F, |
896 | | 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F, |
897 | | 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F, |
898 | | 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F, |
899 | | 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F, |
900 | | 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F, |
901 | | 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F, |
902 | | 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F, |
903 | | 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F, |
904 | | 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F, |
905 | | 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F, |
906 | | 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F, |
907 | | 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F, |
908 | | 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F, |
909 | | 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F, |
910 | | 0.82788260F, 0.88168307F, 0.9389798F, 1.F, |
911 | | }; |
912 | | |
913 | | /* this is for per-channel noise normalization */ |
914 | 0 | static int apsort(const void *a, const void *b){ |
915 | 0 | float f1=**(float**)a; |
916 | 0 | float f2=**(float**)b; |
917 | 0 | return (f1<f2)-(f1>f2); |
918 | 0 | } |
919 | | |
920 | | static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct, |
921 | 0 | float *floor, int *flag, int i, int jn){ |
922 | 0 | int j; |
923 | 0 | for(j=0;j<jn;j++){ |
924 | 0 | float point = j>=limit-i ? postpoint : prepoint; |
925 | 0 | float r = fabs(mdct[j])/floor[j]; |
926 | 0 | if(r<point) |
927 | 0 | flag[j]=0; |
928 | 0 | else |
929 | 0 | flag[j]=1; |
930 | 0 | } |
931 | 0 | } |
932 | | |
933 | | /* Overload/Side effect: On input, the *q vector holds either the |
934 | | quantized energy (for elements with the flag set) or the absolute |
935 | | values of the *r vector (for elements with flag unset). On output, |
936 | | *q holds the quantized energy for all elements */ |
937 | 0 | static float noise_normalize(vorbis_look_psy *p, int limit, float *r, float *q, float *f, int *flags, float acc, int i, int n, int *out){ |
938 | |
|
939 | 0 | vorbis_info_psy *vi=p->vi; |
940 | 0 | float **sort = alloca(n*sizeof(*sort)); |
941 | 0 | int j,count=0; |
942 | 0 | int start = (vi->normal_p ? vi->normal_start-i : n); |
943 | 0 | if(start>n)start=n; |
944 | | |
945 | | /* force classic behavior where only energy in the current band is considered */ |
946 | 0 | acc=0.f; |
947 | | |
948 | | /* still responsible for populating *out where noise norm not in |
949 | | effect. There's no need to [re]populate *q in these areas */ |
950 | 0 | for(j=0;j<start;j++){ |
951 | 0 | if(!flags || !flags[j]){ /* lossless coupling already quantized. |
952 | | Don't touch; requantizing based on |
953 | | energy would be incorrect. */ |
954 | 0 | float ve = q[j]/f[j]; |
955 | 0 | if(r[j]<0) |
956 | 0 | out[j] = -rint(sqrt(ve)); |
957 | 0 | else |
958 | 0 | out[j] = rint(sqrt(ve)); |
959 | 0 | } |
960 | 0 | } |
961 | | |
962 | | /* sort magnitudes for noise norm portion of partition */ |
963 | 0 | for(;j<n;j++){ |
964 | 0 | if(!flags || !flags[j]){ /* can't noise norm elements that have |
965 | | already been loslessly coupled; we can |
966 | | only account for their energy error */ |
967 | 0 | float ve = q[j]/f[j]; |
968 | | /* Despite all the new, more capable coupling code, for now we |
969 | | implement noise norm as it has been up to this point. Only |
970 | | consider promotions to unit magnitude from 0. In addition |
971 | | the only energy error counted is quantizations to zero. */ |
972 | | /* also-- the original point code only applied noise norm at > pointlimit */ |
973 | 0 | if(ve<.25f && (!flags || j>=limit-i)){ |
974 | 0 | acc += ve; |
975 | 0 | sort[count++]=q+j; /* q is fabs(r) for unflagged element */ |
976 | 0 | }else{ |
977 | | /* For now: no acc adjustment for nonzero quantization. populate *out and q as this value is final. */ |
978 | 0 | if(r[j]<0) |
979 | 0 | out[j] = -rint(sqrt(ve)); |
980 | 0 | else |
981 | 0 | out[j] = rint(sqrt(ve)); |
982 | 0 | q[j] = out[j]*out[j]*f[j]; |
983 | 0 | } |
984 | 0 | }/* else{ |
985 | | again, no energy adjustment for error in nonzero quant-- for now |
986 | | }*/ |
987 | 0 | } |
988 | |
|
989 | 0 | if(count){ |
990 | | /* noise norm to do */ |
991 | 0 | qsort(sort,count,sizeof(*sort),apsort); |
992 | 0 | for(j=0;j<count;j++){ |
993 | 0 | int k=sort[j]-q; |
994 | 0 | if(acc>=vi->normal_thresh){ |
995 | 0 | out[k]=unitnorm(r[k]); |
996 | 0 | acc-=1.f; |
997 | 0 | q[k]=f[k]; |
998 | 0 | }else{ |
999 | 0 | out[k]=0; |
1000 | 0 | q[k]=0.f; |
1001 | 0 | } |
1002 | 0 | } |
1003 | 0 | } |
1004 | |
|
1005 | 0 | return acc; |
1006 | 0 | } |
1007 | | |
1008 | | /* Noise normalization, quantization and coupling are not wholly |
1009 | | seperable processes in depth>1 coupling. */ |
1010 | | void _vp_couple_quantize_normalize(int blobno, |
1011 | | vorbis_info_psy_global *g, |
1012 | | vorbis_look_psy *p, |
1013 | | vorbis_info_mapping0 *vi, |
1014 | | float **mdct, |
1015 | | int **iwork, |
1016 | | int *nonzero, |
1017 | | int sliding_lowpass, |
1018 | 0 | int ch){ |
1019 | |
|
1020 | 0 | int i; |
1021 | 0 | int n = p->n; |
1022 | 0 | int partition=(p->vi->normal_p ? p->vi->normal_partition : 16); |
1023 | 0 | int limit = g->coupling_pointlimit[p->vi->blockflag][blobno]; |
1024 | 0 | float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]]; |
1025 | 0 | float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]]; |
1026 | | #if 0 |
1027 | | float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */ |
1028 | | #endif |
1029 | | |
1030 | | /* mdct is our raw mdct output, floor not removed. */ |
1031 | | /* inout passes in the ifloor, passes back quantized result */ |
1032 | | |
1033 | | /* unquantized energy (negative indicates amplitude has negative sign) */ |
1034 | 0 | float **raw = alloca(ch*sizeof(*raw)); |
1035 | | |
1036 | | /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */ |
1037 | 0 | float **quant = alloca(ch*sizeof(*quant)); |
1038 | | |
1039 | | /* floor energy */ |
1040 | 0 | float **floor = alloca(ch*sizeof(*floor)); |
1041 | | |
1042 | | /* flags indicating raw/quantized status of elements in raw vector */ |
1043 | 0 | int **flag = alloca(ch*sizeof(*flag)); |
1044 | | |
1045 | | /* non-zero flag working vector */ |
1046 | 0 | int *nz = alloca(ch*sizeof(*nz)); |
1047 | | |
1048 | | /* energy surplus/defecit tracking */ |
1049 | 0 | float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc)); |
1050 | | |
1051 | | /* The threshold of a stereo is changed with the size of n */ |
1052 | 0 | if(n > 1000) |
1053 | 0 | postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]]; |
1054 | |
|
1055 | 0 | raw[0] = alloca(ch*partition*sizeof(**raw)); |
1056 | 0 | quant[0] = alloca(ch*partition*sizeof(**quant)); |
1057 | 0 | floor[0] = alloca(ch*partition*sizeof(**floor)); |
1058 | 0 | flag[0] = alloca(ch*partition*sizeof(**flag)); |
1059 | |
|
1060 | 0 | for(i=1;i<ch;i++){ |
1061 | 0 | raw[i] = &raw[0][partition*i]; |
1062 | 0 | quant[i] = &quant[0][partition*i]; |
1063 | 0 | floor[i] = &floor[0][partition*i]; |
1064 | 0 | flag[i] = &flag[0][partition*i]; |
1065 | 0 | } |
1066 | 0 | for(i=0;i<ch+vi->coupling_steps;i++) |
1067 | 0 | acc[i]=0.f; |
1068 | |
|
1069 | 0 | for(i=0;i<n;i+=partition){ |
1070 | 0 | int k,j,jn = partition > n-i ? n-i : partition; |
1071 | 0 | int step,track = 0; |
1072 | |
|
1073 | 0 | memcpy(nz,nonzero,sizeof(*nz)*ch); |
1074 | | |
1075 | | /* prefill */ |
1076 | 0 | memset(flag[0],0,ch*partition*sizeof(**flag)); |
1077 | 0 | for(k=0;k<ch;k++){ |
1078 | 0 | int *iout = &iwork[k][i]; |
1079 | 0 | if(nz[k]){ |
1080 | |
|
1081 | 0 | for(j=0;j<jn;j++) |
1082 | 0 | floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]]; |
1083 | |
|
1084 | 0 | flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn); |
1085 | |
|
1086 | 0 | for(j=0;j<jn;j++){ |
1087 | 0 | quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j]; |
1088 | 0 | if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f; |
1089 | 0 | floor[k][j]*=floor[k][j]; |
1090 | 0 | } |
1091 | |
|
1092 | 0 | acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout); |
1093 | |
|
1094 | 0 | }else{ |
1095 | 0 | for(j=0;j<jn;j++){ |
1096 | 0 | floor[k][j] = 1e-10f; |
1097 | 0 | raw[k][j] = 0.f; |
1098 | 0 | quant[k][j] = 0.f; |
1099 | 0 | flag[k][j] = 0; |
1100 | 0 | iout[j]=0; |
1101 | 0 | } |
1102 | 0 | acc[track]=0.f; |
1103 | 0 | } |
1104 | 0 | track++; |
1105 | 0 | } |
1106 | | |
1107 | | /* coupling */ |
1108 | 0 | for(step=0;step<vi->coupling_steps;step++){ |
1109 | 0 | int Mi = vi->coupling_mag[step]; |
1110 | 0 | int Ai = vi->coupling_ang[step]; |
1111 | 0 | int *iM = &iwork[Mi][i]; |
1112 | 0 | int *iA = &iwork[Ai][i]; |
1113 | 0 | float *reM = raw[Mi]; |
1114 | 0 | float *reA = raw[Ai]; |
1115 | 0 | float *qeM = quant[Mi]; |
1116 | 0 | float *qeA = quant[Ai]; |
1117 | 0 | float *floorM = floor[Mi]; |
1118 | 0 | float *floorA = floor[Ai]; |
1119 | 0 | int *fM = flag[Mi]; |
1120 | 0 | int *fA = flag[Ai]; |
1121 | |
|
1122 | 0 | if(nz[Mi] || nz[Ai]){ |
1123 | 0 | nz[Mi] = nz[Ai] = 1; |
1124 | |
|
1125 | 0 | for(j=0;j<jn;j++){ |
1126 | |
|
1127 | 0 | if(j<sliding_lowpass-i){ |
1128 | 0 | if(fM[j] || fA[j]){ |
1129 | | /* lossless coupling */ |
1130 | |
|
1131 | 0 | reM[j] = fabs(reM[j])+fabs(reA[j]); |
1132 | 0 | qeM[j] = qeM[j]+qeA[j]; |
1133 | 0 | fM[j]=fA[j]=1; |
1134 | | |
1135 | | /* couple iM/iA */ |
1136 | 0 | { |
1137 | 0 | int A = iM[j]; |
1138 | 0 | int B = iA[j]; |
1139 | |
|
1140 | 0 | if(abs(A)>abs(B)){ |
1141 | 0 | iA[j]=(A>0?A-B:B-A); |
1142 | 0 | }else{ |
1143 | 0 | iA[j]=(B>0?A-B:B-A); |
1144 | 0 | iM[j]=B; |
1145 | 0 | } |
1146 | | |
1147 | | /* collapse two equivalent tuples to one */ |
1148 | 0 | if(iA[j]>=abs(iM[j])*2){ |
1149 | 0 | iA[j]= -iA[j]; |
1150 | 0 | iM[j]= -iM[j]; |
1151 | 0 | } |
1152 | |
|
1153 | 0 | } |
1154 | |
|
1155 | 0 | }else{ |
1156 | | /* lossy (point) coupling */ |
1157 | 0 | if(j<limit-i){ |
1158 | | /* dipole */ |
1159 | 0 | reM[j] += reA[j]; |
1160 | 0 | qeM[j] = fabs(reM[j]); |
1161 | 0 | }else{ |
1162 | | #if 0 |
1163 | | /* AoTuV */ |
1164 | | /** @ M2 ** |
1165 | | The boost problem by the combination of noise normalization and point stereo is eased. |
1166 | | However, this is a temporary patch. |
1167 | | by Aoyumi @ 2004/04/18 |
1168 | | */ |
1169 | | float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit))); |
1170 | | /* elliptical */ |
1171 | | if(reM[j]+reA[j]<0){ |
1172 | | reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); |
1173 | | }else{ |
1174 | | reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); |
1175 | | } |
1176 | | #else |
1177 | | /* elliptical */ |
1178 | 0 | if(reM[j]+reA[j]<0){ |
1179 | 0 | reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j])); |
1180 | 0 | }else{ |
1181 | 0 | reM[j] = (qeM[j] = fabs(reM[j])+fabs(reA[j])); |
1182 | 0 | } |
1183 | 0 | #endif |
1184 | |
|
1185 | 0 | } |
1186 | 0 | reA[j]=qeA[j]=0.f; |
1187 | 0 | fA[j]=1; |
1188 | 0 | iA[j]=0; |
1189 | 0 | } |
1190 | 0 | } |
1191 | 0 | floorM[j]=floorA[j]=floorM[j]+floorA[j]; |
1192 | 0 | } |
1193 | | /* normalize the resulting mag vector */ |
1194 | 0 | acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM); |
1195 | 0 | track++; |
1196 | 0 | } |
1197 | 0 | } |
1198 | 0 | } |
1199 | |
|
1200 | 0 | for(i=0;i<vi->coupling_steps;i++){ |
1201 | | /* make sure coupling a zero and a nonzero channel results in two |
1202 | | nonzero channels. */ |
1203 | 0 | if(nonzero[vi->coupling_mag[i]] || |
1204 | 0 | nonzero[vi->coupling_ang[i]]){ |
1205 | 0 | nonzero[vi->coupling_mag[i]]=1; |
1206 | 0 | nonzero[vi->coupling_ang[i]]=1; |
1207 | 0 | } |
1208 | 0 | } |
1209 | 0 | } |