diff options
author | Sean Bartell <wingedtachikoma@gmail.com> | 2011-06-25 21:32:25 -0400 |
---|---|---|
committer | Nils Wallménius <nils@rockbox.org> | 2012-04-25 22:13:20 +0200 |
commit | f40bfc9267b13b54e6379dfe7539447662879d24 (patch) | |
tree | 9b20069d5e62809ff434061ad730096836f916f2 /lib/rbcodec/codecs/libspeex/vorbis_psy.c | |
parent | a0009907de7a0107d49040d8a180f140e2eff299 (diff) | |
download | rockbox-f40bfc9267b13b54e6379dfe7539447662879d24.tar.gz rockbox-f40bfc9267b13b54e6379dfe7539447662879d24.zip |
Add codecs to librbcodec.
Change-Id: Id7f4717d51ed02d67cb9f9cb3c0ada4a81843f97
Reviewed-on: http://gerrit.rockbox.org/137
Reviewed-by: Nils Wallménius <nils@rockbox.org>
Tested-by: Nils Wallménius <nils@rockbox.org>
Diffstat (limited to 'lib/rbcodec/codecs/libspeex/vorbis_psy.c')
-rw-r--r-- | lib/rbcodec/codecs/libspeex/vorbis_psy.c | 508 |
1 files changed, 508 insertions, 0 deletions
diff --git a/lib/rbcodec/codecs/libspeex/vorbis_psy.c b/lib/rbcodec/codecs/libspeex/vorbis_psy.c new file mode 100644 index 0000000000..2032bf63e2 --- /dev/null +++ b/lib/rbcodec/codecs/libspeex/vorbis_psy.c | |||
@@ -0,0 +1,508 @@ | |||
1 | /* Copyright (C) 2005 Jean-Marc Valin, CSIRO, Christopher Montgomery | ||
2 | File: vorbis_psy.c | ||
3 | |||
4 | Redistribution and use in source and binary forms, with or without | ||
5 | modification, are permitted provided that the following conditions | ||
6 | are met: | ||
7 | |||
8 | - Redistributions of source code must retain the above copyright | ||
9 | notice, this list of conditions and the following disclaimer. | ||
10 | |||
11 | - Redistributions in binary form must reproduce the above copyright | ||
12 | notice, this list of conditions and the following disclaimer in the | ||
13 | documentation and/or other materials provided with the distribution. | ||
14 | |||
15 | - Neither the name of the Xiph.org Foundation nor the names of its | ||
16 | contributors may be used to endorse or promote products derived from | ||
17 | this software without specific prior written permission. | ||
18 | |||
19 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | ||
20 | ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | ||
21 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | ||
22 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR | ||
23 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, | ||
24 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | ||
25 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | ||
26 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | ||
27 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | ||
28 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | ||
29 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | ||
30 | */ | ||
31 | |||
32 | #ifdef HAVE_CONFIG_H | ||
33 | #include "config-speex.h" | ||
34 | #endif | ||
35 | |||
36 | #ifdef VORBIS_PSYCHO | ||
37 | |||
38 | #include "arch.h" | ||
39 | #include "smallft.h" | ||
40 | #include "lpc.h" | ||
41 | #include "vorbis_psy.h" | ||
42 | |||
43 | #include <stdlib.h> | ||
44 | #include <math.h> | ||
45 | #include <string.h> | ||
46 | |||
47 | /* psychoacoustic setup ********************************************/ | ||
48 | |||
49 | static VorbisPsyInfo example_tuning = { | ||
50 | |||
51 | .5,.5, | ||
52 | 3,3,25, | ||
53 | |||
54 | /*63 125 250 500 1k 2k 4k 8k 16k*/ | ||
55 | // vorbis mode 4 style | ||
56 | //{-32,-32,-32,-32,-28,-24,-22,-20,-20, -20, -20, -8, -6, -6, -6, -6, -6}, | ||
57 | { -4, -6, -6, -6, -6, -6, -6, -6, -8, -8,-10,-10, -8, -6, -4, -4, -2}, | ||
58 | |||
59 | { | ||
60 | 0, 1, 2, 3, 4, 5, 5, 5, /* 7dB */ | ||
61 | 6, 6, 6, 5, 4, 4, 4, 4, /* 15dB */ | ||
62 | 4, 4, 5, 5, 5, 6, 6, 6, /* 23dB */ | ||
63 | 7, 7, 7, 8, 8, 8, 9, 10, /* 31dB */ | ||
64 | 11,12,13,14,15,16,17, 18, /* 39dB */ | ||
65 | } | ||
66 | |||
67 | }; | ||
68 | |||
69 | |||
70 | |||
71 | /* there was no great place to put this.... */ | ||
72 | #include <stdio.h> | ||
73 | static void _analysis_output(char *base,int i,float *v,int n,int bark,int dB){ | ||
74 | int j; | ||
75 | FILE *of; | ||
76 | char buffer[80]; | ||
77 | |||
78 | sprintf(buffer,"%s_%d.m",base,i); | ||
79 | of=fopen(buffer,"w"); | ||
80 | |||
81 | if(!of)perror("failed to open data dump file"); | ||
82 | |||
83 | for(j=0;j<n;j++){ | ||
84 | if(bark){ | ||
85 | float b=toBARK((4000.f*j/n)+.25); | ||
86 | fprintf(of,"%f ",b); | ||
87 | }else | ||
88 | fprintf(of,"%f ",(double)j); | ||
89 | |||
90 | if(dB){ | ||
91 | float val; | ||
92 | if(v[j]==0.) | ||
93 | val=-140.; | ||
94 | else | ||
95 | val=todB(v[j]); | ||
96 | fprintf(of,"%f\n",val); | ||
97 | }else{ | ||
98 | fprintf(of,"%f\n",v[j]); | ||
99 | } | ||
100 | } | ||
101 | fclose(of); | ||
102 | } | ||
103 | |||
104 | static void bark_noise_hybridmp(int n,const long *b, | ||
105 | const float *f, | ||
106 | float *noise, | ||
107 | const float offset, | ||
108 | const int fixed){ | ||
109 | |||
110 | float *N=alloca(n*sizeof(*N)); | ||
111 | float *X=alloca(n*sizeof(*N)); | ||
112 | float *XX=alloca(n*sizeof(*N)); | ||
113 | float *Y=alloca(n*sizeof(*N)); | ||
114 | float *XY=alloca(n*sizeof(*N)); | ||
115 | |||
116 | float tN, tX, tXX, tY, tXY; | ||
117 | int i; | ||
118 | |||
119 | int lo, hi; | ||
120 | float R, A, B, D; | ||
121 | float w, x, y; | ||
122 | |||
123 | tN = tX = tXX = tY = tXY = 0.f; | ||
124 | |||
125 | y = f[0] + offset; | ||
126 | if (y < 1.f) y = 1.f; | ||
127 | |||
128 | w = y * y * .5; | ||
129 | |||
130 | tN += w; | ||
131 | tX += w; | ||
132 | tY += w * y; | ||
133 | |||
134 | N[0] = tN; | ||
135 | X[0] = tX; | ||
136 | XX[0] = tXX; | ||
137 | Y[0] = tY; | ||
138 | XY[0] = tXY; | ||
139 | |||
140 | for (i = 1, x = 1.f; i < n; i++, x += 1.f) { | ||
141 | |||
142 | y = f[i] + offset; | ||
143 | if (y < 1.f) y = 1.f; | ||
144 | |||
145 | w = y * y; | ||
146 | |||
147 | tN += w; | ||
148 | tX += w * x; | ||
149 | tXX += w * x * x; | ||
150 | tY += w * y; | ||
151 | tXY += w * x * y; | ||
152 | |||
153 | N[i] = tN; | ||
154 | X[i] = tX; | ||
155 | XX[i] = tXX; | ||
156 | Y[i] = tY; | ||
157 | XY[i] = tXY; | ||
158 | } | ||
159 | |||
160 | for (i = 0, x = 0.f;; i++, x += 1.f) { | ||
161 | |||
162 | lo = b[i] >> 16; | ||
163 | if( lo>=0 ) break; | ||
164 | hi = b[i] & 0xffff; | ||
165 | |||
166 | tN = N[hi] + N[-lo]; | ||
167 | tX = X[hi] - X[-lo]; | ||
168 | tXX = XX[hi] + XX[-lo]; | ||
169 | tY = Y[hi] + Y[-lo]; | ||
170 | tXY = XY[hi] - XY[-lo]; | ||
171 | |||
172 | A = tY * tXX - tX * tXY; | ||
173 | B = tN * tXY - tX * tY; | ||
174 | D = tN * tXX - tX * tX; | ||
175 | R = (A + x * B) / D; | ||
176 | if (R < 0.f) | ||
177 | R = 0.f; | ||
178 | |||
179 | noise[i] = R - offset; | ||
180 | } | ||
181 | |||
182 | for ( ;; i++, x += 1.f) { | ||
183 | |||
184 | lo = b[i] >> 16; | ||
185 | hi = b[i] & 0xffff; | ||
186 | if(hi>=n)break; | ||
187 | |||
188 | tN = N[hi] - N[lo]; | ||
189 | tX = X[hi] - X[lo]; | ||
190 | tXX = XX[hi] - XX[lo]; | ||
191 | tY = Y[hi] - Y[lo]; | ||
192 | tXY = XY[hi] - XY[lo]; | ||
193 | |||
194 | A = tY * tXX - tX * tXY; | ||
195 | B = tN * tXY - tX * tY; | ||
196 | D = tN * tXX - tX * tX; | ||
197 | R = (A + x * B) / D; | ||
198 | if (R < 0.f) R = 0.f; | ||
199 | |||
200 | noise[i] = R - offset; | ||
201 | } | ||
202 | for ( ; i < n; i++, x += 1.f) { | ||
203 | |||
204 | R = (A + x * B) / D; | ||
205 | if (R < 0.f) R = 0.f; | ||
206 | |||
207 | noise[i] = R - offset; | ||
208 | } | ||
209 | |||
210 | if (fixed <= 0) return; | ||
211 | |||
212 | for (i = 0, x = 0.f;; i++, x += 1.f) { | ||
213 | hi = i + fixed / 2; | ||
214 | lo = hi - fixed; | ||
215 | if(lo>=0)break; | ||
216 | |||
217 | tN = N[hi] + N[-lo]; | ||
218 | tX = X[hi] - X[-lo]; | ||
219 | tXX = XX[hi] + XX[-lo]; | ||
220 | tY = Y[hi] + Y[-lo]; | ||
221 | tXY = XY[hi] - XY[-lo]; | ||
222 | |||
223 | |||
224 | A = tY * tXX - tX * tXY; | ||
225 | B = tN * tXY - tX * tY; | ||
226 | D = tN * tXX - tX * tX; | ||
227 | R = (A + x * B) / D; | ||
228 | |||
229 | if (R - offset < noise[i]) noise[i] = R - offset; | ||
230 | } | ||
231 | for ( ;; i++, x += 1.f) { | ||
232 | |||
233 | hi = i + fixed / 2; | ||
234 | lo = hi - fixed; | ||
235 | if(hi>=n)break; | ||
236 | |||
237 | tN = N[hi] - N[lo]; | ||
238 | tX = X[hi] - X[lo]; | ||
239 | tXX = XX[hi] - XX[lo]; | ||
240 | tY = Y[hi] - Y[lo]; | ||
241 | tXY = XY[hi] - XY[lo]; | ||
242 | |||
243 | A = tY * tXX - tX * tXY; | ||
244 | B = tN * tXY - tX * tY; | ||
245 | D = tN * tXX - tX * tX; | ||
246 | R = (A + x * B) / D; | ||
247 | |||
248 | if (R - offset < noise[i]) noise[i] = R - offset; | ||
249 | } | ||
250 | for ( ; i < n; i++, x += 1.f) { | ||
251 | R = (A + x * B) / D; | ||
252 | if (R - offset < noise[i]) noise[i] = R - offset; | ||
253 | } | ||
254 | } | ||
255 | |||
256 | static void _vp_noisemask(VorbisPsy *p, | ||
257 | float *logfreq, | ||
258 | float *logmask){ | ||
259 | |||
260 | int i,n=p->n/2; | ||
261 | float *work=alloca(n*sizeof(*work)); | ||
262 | |||
263 | bark_noise_hybridmp(n,p->bark,logfreq,logmask, | ||
264 | 140.,-1); | ||
265 | |||
266 | for(i=0;i<n;i++)work[i]=logfreq[i]-logmask[i]; | ||
267 | |||
268 | bark_noise_hybridmp(n,p->bark,work,logmask,0., | ||
269 | p->vi->noisewindowfixed); | ||
270 | |||
271 | for(i=0;i<n;i++)work[i]=logfreq[i]-work[i]; | ||
272 | |||
273 | { | ||
274 | static int seq=0; | ||
275 | |||
276 | float work2[n]; | ||
277 | for(i=0;i<n;i++){ | ||
278 | work2[i]=logmask[i]+work[i]; | ||
279 | } | ||
280 | |||
281 | //_analysis_output("logfreq",seq,logfreq,n,0,0); | ||
282 | //_analysis_output("median",seq,work,n,0,0); | ||
283 | //_analysis_output("envelope",seq,work2,n,0,0); | ||
284 | seq++; | ||
285 | } | ||
286 | |||
287 | for(i=0;i<n;i++){ | ||
288 | int dB=logmask[i]+.5; | ||
289 | if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1; | ||
290 | if(dB<0)dB=0; | ||
291 | logmask[i]= work[i]+p->vi->noisecompand[dB]+p->noiseoffset[i]; | ||
292 | } | ||
293 | |||
294 | } | ||
295 | |||
296 | VorbisPsy *vorbis_psy_init(int rate, int n) | ||
297 | { | ||
298 | long i,j,lo=-99,hi=1; | ||
299 | VorbisPsy *p = speex_alloc(sizeof(VorbisPsy)); | ||
300 | memset(p,0,sizeof(*p)); | ||
301 | |||
302 | p->n = n; | ||
303 | spx_drft_init(&p->lookup, n); | ||
304 | p->bark = speex_alloc(n*sizeof(*p->bark)); | ||
305 | p->rate=rate; | ||
306 | p->vi = &example_tuning; | ||
307 | |||
308 | /* BH4 window */ | ||
309 | p->window = speex_alloc(sizeof(*p->window)*n); | ||
310 | float a0 = .35875f; | ||
311 | float a1 = .48829f; | ||
312 | float a2 = .14128f; | ||
313 | float a3 = .01168f; | ||
314 | for(i=0;i<n;i++) | ||
315 | p->window[i] = //a0 - a1*cos(2.*M_PI/n*(i+.5)) + a2*cos(4.*M_PI/n*(i+.5)) - a3*cos(6.*M_PI/n*(i+.5)); | ||
316 | sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI); | ||
317 | /* bark scale lookups */ | ||
318 | for(i=0;i<n;i++){ | ||
319 | float bark=toBARK(rate/(2*n)*i); | ||
320 | |||
321 | for(;lo+p->vi->noisewindowlomin<i && | ||
322 | toBARK(rate/(2*n)*lo)<(bark-p->vi->noisewindowlo);lo++); | ||
323 | |||
324 | for(;hi<=n && (hi<i+p->vi->noisewindowhimin || | ||
325 | toBARK(rate/(2*n)*hi)<(bark+p->vi->noisewindowhi));hi++); | ||
326 | |||
327 | p->bark[i]=((lo-1)<<16)+(hi-1); | ||
328 | |||
329 | } | ||
330 | |||
331 | /* set up rolling noise median */ | ||
332 | p->noiseoffset=speex_alloc(n*sizeof(*p->noiseoffset)); | ||
333 | |||
334 | for(i=0;i<n;i++){ | ||
335 | float halfoc=toOC((i+.5)*rate/(2.*n))*2.; | ||
336 | int inthalfoc; | ||
337 | float del; | ||
338 | |||
339 | if(halfoc<0)halfoc=0; | ||
340 | if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1; | ||
341 | inthalfoc=(int)halfoc; | ||
342 | del=halfoc-inthalfoc; | ||
343 | |||
344 | p->noiseoffset[i]= | ||
345 | p->vi->noiseoff[inthalfoc]*(1.-del) + | ||
346 | p->vi->noiseoff[inthalfoc+1]*del; | ||
347 | |||
348 | } | ||
349 | #if 0 | ||
350 | _analysis_output_always("noiseoff0",ls,p->noiseoffset,n,1,0,0); | ||
351 | #endif | ||
352 | |||
353 | return p; | ||
354 | } | ||
355 | |||
356 | void vorbis_psy_destroy(VorbisPsy *p) | ||
357 | { | ||
358 | if(p){ | ||
359 | spx_drft_clear(&p->lookup); | ||
360 | if(p->bark) | ||
361 | speex_free(p->bark); | ||
362 | if(p->noiseoffset) | ||
363 | speex_free(p->noiseoffset); | ||
364 | if(p->window) | ||
365 | speex_free(p->window); | ||
366 | memset(p,0,sizeof(*p)); | ||
367 | speex_free(p); | ||
368 | } | ||
369 | } | ||
370 | |||
371 | void compute_curve(VorbisPsy *psy, float *audio, float *curve) | ||
372 | { | ||
373 | int i; | ||
374 | float work[psy->n]; | ||
375 | |||
376 | float scale=4.f/psy->n; | ||
377 | float scale_dB; | ||
378 | |||
379 | scale_dB=todB(scale); | ||
380 | |||
381 | /* window the PCM data; use a BH4 window, not vorbis */ | ||
382 | for(i=0;i<psy->n;i++) | ||
383 | work[i]=audio[i] * psy->window[i]; | ||
384 | |||
385 | { | ||
386 | static int seq=0; | ||
387 | |||
388 | //_analysis_output("win",seq,work,psy->n,0,0); | ||
389 | |||
390 | seq++; | ||
391 | } | ||
392 | |||
393 | /* FFT yields more accurate tonal estimation (not phase sensitive) */ | ||
394 | spx_drft_forward(&psy->lookup,work); | ||
395 | |||
396 | /* magnitudes */ | ||
397 | work[0]=scale_dB+todB(work[0]); | ||
398 | for(i=1;i<psy->n-1;i+=2){ | ||
399 | float temp = work[i]*work[i] + work[i+1]*work[i+1]; | ||
400 | work[(i+1)>>1] = scale_dB+.5f * todB(temp); | ||
401 | } | ||
402 | |||
403 | /* derive a noise curve */ | ||
404 | _vp_noisemask(psy,work,curve); | ||
405 | #define SIDEL 12 | ||
406 | for (i=0;i<SIDEL;i++) | ||
407 | { | ||
408 | curve[i]=curve[SIDEL]; | ||
409 | } | ||
410 | #define SIDEH 12 | ||
411 | for (i=0;i<SIDEH;i++) | ||
412 | { | ||
413 | curve[(psy->n>>1)-i-1]=curve[(psy->n>>1)-SIDEH]; | ||
414 | } | ||
415 | for(i=0;i<((psy->n)>>1);i++) | ||
416 | curve[i] = fromdB(1.2*curve[i]+.2*i); | ||
417 | //curve[i] = fromdB(0.8*curve[i]+.35*i); | ||
418 | //curve[i] = fromdB(0.9*curve[i])*pow(1.0*i+45,1.3); | ||
419 | } | ||
420 | |||
421 | /* Transform a masking curve (power spectrum) into a pole-zero filter */ | ||
422 | void curve_to_lpc(VorbisPsy *psy, float *curve, float *awk1, float *awk2, int ord) | ||
423 | { | ||
424 | int i; | ||
425 | float ac[psy->n]; | ||
426 | float tmp; | ||
427 | int len = psy->n >> 1; | ||
428 | for (i=0;i<2*len;i++) | ||
429 | ac[i] = 0; | ||
430 | for (i=1;i<len;i++) | ||
431 | ac[2*i-1] = curve[i]; | ||
432 | ac[0] = curve[0]; | ||
433 | ac[2*len-1] = curve[len-1]; | ||
434 | |||
435 | spx_drft_backward(&psy->lookup, ac); | ||
436 | _spx_lpc(awk1, ac, ord); | ||
437 | tmp = 1.; | ||
438 | for (i=0;i<ord;i++) | ||
439 | { | ||
440 | tmp *= .99; | ||
441 | awk1[i] *= tmp; | ||
442 | } | ||
443 | #if 0 | ||
444 | for (i=0;i<ord;i++) | ||
445 | awk2[i] = 0; | ||
446 | #else | ||
447 | /* Use the second (awk2) filter to correct the first one */ | ||
448 | for (i=0;i<2*len;i++) | ||
449 | ac[i] = 0; | ||
450 | for (i=0;i<ord;i++) | ||
451 | ac[i+1] = awk1[i]; | ||
452 | ac[0] = 1; | ||
453 | spx_drft_forward(&psy->lookup, ac); | ||
454 | /* Compute (power) response of awk1 (all zero) */ | ||
455 | ac[0] *= ac[0]; | ||
456 | for (i=1;i<len;i++) | ||
457 | ac[i] = ac[2*i-1]*ac[2*i-1] + ac[2*i]*ac[2*i]; | ||
458 | ac[len] = ac[2*len-1]*ac[2*len-1]; | ||
459 | /* Compute correction required */ | ||
460 | for (i=0;i<len;i++) | ||
461 | curve[i] = 1. / (1e-6f+curve[i]*ac[i]); | ||
462 | |||
463 | for (i=0;i<2*len;i++) | ||
464 | ac[i] = 0; | ||
465 | for (i=1;i<len;i++) | ||
466 | ac[2*i-1] = curve[i]; | ||
467 | ac[0] = curve[0]; | ||
468 | ac[2*len-1] = curve[len-1]; | ||
469 | |||
470 | spx_drft_backward(&psy->lookup, ac); | ||
471 | _spx_lpc(awk2, ac, ord); | ||
472 | tmp = 1; | ||
473 | for (i=0;i<ord;i++) | ||
474 | { | ||
475 | tmp *= .99; | ||
476 | awk2[i] *= tmp; | ||
477 | } | ||
478 | #endif | ||
479 | } | ||
480 | |||
481 | #if 0 | ||
482 | #include <stdio.h> | ||
483 | #include <math.h> | ||
484 | |||
485 | #define ORDER 10 | ||
486 | #define CURVE_SIZE 24 | ||
487 | |||
488 | int main() | ||
489 | { | ||
490 | int i; | ||
491 | float curve[CURVE_SIZE]; | ||
492 | float awk1[ORDER], awk2[ORDER]; | ||
493 | for (i=0;i<CURVE_SIZE;i++) | ||
494 | scanf("%f ", &curve[i]); | ||
495 | for (i=0;i<CURVE_SIZE;i++) | ||
496 | curve[i] = pow(10.f, .1*curve[i]); | ||
497 | curve_to_lpc(curve, CURVE_SIZE, awk1, awk2, ORDER); | ||
498 | for (i=0;i<ORDER;i++) | ||
499 | printf("%f ", awk1[i]); | ||
500 | printf ("\n"); | ||
501 | for (i=0;i<ORDER;i++) | ||
502 | printf("%f ", awk2[i]); | ||
503 | printf ("\n"); | ||
504 | return 0; | ||
505 | } | ||
506 | #endif | ||
507 | |||
508 | #endif | ||