Mercurial > mplayer.hg
annotate libaf/af_hrtf.c @ 14495:bfe061a6ae4c
Sync with 1.858:
1.858: disambiguate wording of -aspect
1.857: more on H.264's quantization parameter
1.856: (was already corrected on previous version)
author | gpoirier |
---|---|
date | Fri, 14 Jan 2005 16:58:39 +0000 |
parents | 0a22a046f0d3 |
children | 71570687c1a3 |
rev | line source |
---|---|
13996 | 1 /* Experimental audio filter that mixes 5.1 and 5.1 with matrix |
2 encoded rear channels into headphone signal using FIR filtering | |
3 with HRTF. | |
4 */ | |
5 //#include <stdio.h> | |
6 #include <stdlib.h> | |
7 #include <string.h> | |
8 #include <unistd.h> | |
9 #include <inttypes.h> | |
10 | |
11 #include <math.h> | |
12 | |
13 #include "af.h" | |
14 #include "dsp.h" | |
15 | |
16 /* HRTF filter coefficients and adjustable parameters */ | |
17 #include "af_hrtf.h" | |
18 | |
19 typedef struct af_hrtf_s { | |
20 /* Lengths */ | |
21 int dlbuflen, hrflen, basslen; | |
22 /* L, C, R, Ls, Rs channels */ | |
23 float *lf, *rf, *lr, *rr, *cf, *cr; | |
24 float *cf_ir, *af_ir, *of_ir, *ar_ir, *or_ir, *cr_ir; | |
25 int cf_o, af_o, of_o, ar_o, or_o, cr_o; | |
26 /* Bass */ | |
27 float *ba_l, *ba_r; | |
28 float *ba_ir; | |
29 /* Whether to matrix decode the rear center channel */ | |
30 int matrix_mode; | |
31 /* Full wave rectified amplitude used to steer the active matrix | |
32 decoding of center rear channel */ | |
33 float lr_fwr, rr_fwr; | |
34 /* Cyclic position on the ring buffer */ | |
35 int cyc_pos; | |
36 } af_hrtf_t; | |
37 | |
38 /* Convolution on a ring buffer | |
39 * nx: length of the ring buffer | |
40 * nk: length of the convolution kernel | |
41 * sx: ring buffer | |
42 * sk: convolution kernel | |
43 * offset: offset on the ring buffer, can be | |
44 */ | |
45 static float conv(const int nx, const int nk, float *sx, float *sk, | |
46 const int offset) | |
47 { | |
48 /* k = reminder of offset / nx */ | |
49 int k = offset >= 0 ? offset % nx : nx + (offset % nx); | |
50 | |
51 if(nk + k <= nx) | |
14275
de13fd557440
less namespace pollution #2 (prefixed globals in filter.c with af_filter_)
alex
parents:
14245
diff
changeset
|
52 return af_filter_fir(nk, sx + k, sk); |
13996 | 53 else |
14275
de13fd557440
less namespace pollution #2 (prefixed globals in filter.c with af_filter_)
alex
parents:
14245
diff
changeset
|
54 return af_filter_fir(nk + k - nx, sx, sk + nx - k) + |
de13fd557440
less namespace pollution #2 (prefixed globals in filter.c with af_filter_)
alex
parents:
14245
diff
changeset
|
55 af_filter_fir(nx - k, sx + k, sk); |
13996 | 56 } |
57 | |
58 /* Detect when the impulse response starts (significantly) */ | |
59 int pulse_detect(float *sx) | |
60 { | |
61 /* nmax must be the reference impulse response length (128) minus | |
62 s->hrflen */ | |
63 const int nmax = 128 - HRTFFILTLEN; | |
64 const float thresh = IRTHRESH; | |
65 int i; | |
66 | |
67 for(i = 0; i < nmax; i++) | |
68 if(fabs(sx[i]) > thresh) | |
69 return i; | |
70 return 0; | |
71 } | |
72 | |
73 inline void update_ch(af_hrtf_t *s, short *in, const int k) | |
74 { | |
75 /* Update the full wave rectified total amplutude */ | |
76 s->lr_fwr += abs(in[2]) - fabs(s->lr[k]); | |
77 s->rr_fwr += abs(in[3]) - fabs(s->rr[k]); | |
78 | |
79 s->lf[k] = in[0]; | |
80 s->cf[k] = in[4]; | |
81 s->rf[k] = in[1]; | |
82 s->lr[k] = in[2]; | |
83 s->rr[k] = in[3]; | |
84 | |
85 s->ba_l[k] = in[0] + in[4] + in[2]; | |
86 s->ba_r[k] = in[4] + in[1] + in[3]; | |
87 } | |
88 | |
89 inline void matrix_decode_cr(af_hrtf_t *s, short *in, const int k) | |
90 { | |
91 /* Active matrix decoding of the center rear channel, 1 in the | |
92 denominator is to prevent singularity */ | |
93 float lr_agc = in[2] * (s->lr_fwr + s->rr_fwr) / | |
94 (1 + s->lr_fwr + s->lr_fwr); | |
95 float rr_agc = in[3] * (s->lr_fwr + s->rr_fwr) / | |
96 (1 + s->rr_fwr + s->rr_fwr); | |
97 | |
98 s->cr[k] = (lr_agc + rr_agc) * M_SQRT1_2; | |
99 } | |
100 | |
101 /* Initialization and runtime control */ | |
102 static int control(struct af_instance_s *af, int cmd, void* arg) | |
103 { | |
104 af_hrtf_t *s = af->setup; | |
105 char mode; | |
106 | |
107 switch(cmd) { | |
108 case AF_CONTROL_REINIT: | |
109 af->data->rate = ((af_data_t*)arg)->rate; | |
110 if(af->data->rate != 48000) { | |
14213
3c56b18bbb0c
Make filters request a supported input format instead of failing.
reimar
parents:
13996
diff
changeset
|
111 // automatic samplerate adjustment in the filter chain |
3c56b18bbb0c
Make filters request a supported input format instead of failing.
reimar
parents:
13996
diff
changeset
|
112 // is not yet supported. |
13996 | 113 af_msg(AF_MSG_ERROR, |
114 "[hrtf] ERROR: Sampling rate is not 48000 Hz (%d)!\n", | |
115 af->data->rate); | |
116 return AF_ERROR; | |
117 } | |
118 af->data->nch = ((af_data_t*)arg)->nch; | |
119 if(af->data->nch < 5) { | |
14213
3c56b18bbb0c
Make filters request a supported input format instead of failing.
reimar
parents:
13996
diff
changeset
|
120 af->data->nch = 5; |
13996 | 121 } |
14245 | 122 af->data->format = AF_FORMAT_S16_NE; |
13996 | 123 af->data->bps = 2; |
14213
3c56b18bbb0c
Make filters request a supported input format instead of failing.
reimar
parents:
13996
diff
changeset
|
124 return af_test_output(af, (af_data_t*)arg); |
13996 | 125 case AF_CONTROL_COMMAND_LINE: |
126 sscanf((char*)arg, "%c", &mode); | |
127 switch(mode) { | |
128 case 'm': | |
129 s->matrix_mode = 1; | |
130 break; | |
131 case '0': | |
132 s->matrix_mode = 0; | |
133 break; | |
134 default: | |
135 af_msg(AF_MSG_ERROR, | |
136 "[hrtf] Mode is neither 'm', nor '0' (%c).\n", | |
137 mode); | |
138 return AF_ERROR; | |
139 } | |
140 return AF_OK; | |
141 } | |
142 | |
143 af_msg(AF_MSG_INFO, | |
144 "[hrtf] Using HRTF to mix %s discrete surround into " | |
145 "L, R channels\n", s->matrix_mode ? "5" : "5+1"); | |
146 if(s->matrix_mode) | |
147 af_msg(AF_MSG_INFO, | |
148 "[hrtf] Using active matrix to decode rear center " | |
149 "channel\n"); | |
150 | |
151 return AF_UNKNOWN; | |
152 } | |
153 | |
154 /* Deallocate memory */ | |
155 static void uninit(struct af_instance_s *af) | |
156 { | |
157 if(af->setup) { | |
158 af_hrtf_t *s = af->setup; | |
159 | |
160 if(s->lf) | |
161 free(s->lf); | |
162 if(s->rf) | |
163 free(s->rf); | |
164 if(s->lr) | |
165 free(s->lr); | |
166 if(s->rr) | |
167 free(s->rr); | |
168 if(s->cf) | |
169 free(s->cf); | |
170 if(s->cr) | |
171 free(s->cr); | |
172 if(s->ba_l) | |
173 free(s->ba_l); | |
174 if(s->ba_r) | |
175 free(s->ba_r); | |
176 if(s->ba_ir) | |
177 free(s->ba_ir); | |
178 free(af->setup); | |
179 } | |
180 if(af->data) | |
181 free(af->data); | |
182 } | |
183 | |
184 /* Filter data through filter | |
185 | |
186 Two "tricks" are used to compensate the "color" of the KEMAR data: | |
187 | |
188 1. The KEMAR data is refiltered to ensure that the front L, R channels | |
189 on the same side of the ear are equalized (especially in the high | |
190 frequencies). | |
191 | |
192 2. A bass compensation is introduced to ensure that 0-200 Hz are not | |
193 damped (without any real 3D acoustical image, however). | |
194 */ | |
195 static af_data_t* play(struct af_instance_s *af, af_data_t *data) | |
196 { | |
197 af_hrtf_t *s = af->setup; | |
198 short *in = data->audio; // Input audio data | |
199 short *out = NULL; // Output audio data | |
200 short *end = in + data->len / sizeof(short); // Loop end | |
201 float common, left, right, diff, left_b, right_b; | |
202 const int dblen = s->dlbuflen, hlen = s->hrflen, blen = s->basslen; | |
203 | |
204 if(AF_OK != RESIZE_LOCAL_BUFFER(af, data)) | |
205 return NULL; | |
206 | |
207 out = af->data->audio; | |
208 | |
209 /* MPlayer's 5 channel layout (notation for the variable): | |
210 * | |
211 * 0: L (LF), 1: R (RF), 2: Ls (LR), 3: Rs (RR), 4: C (CF), matrix | |
212 * encoded: Cs (CR) | |
213 * | |
214 * or: L = left, C = center, R = right, F = front, R = rear | |
215 * | |
216 * Filter notation: | |
217 * | |
218 * CF | |
219 * OF AF | |
220 * Ear-> | |
221 * OR AR | |
222 * CR | |
223 * | |
224 * or: C = center, A = same side, O = opposite, F = front, R = rear | |
225 */ | |
226 | |
227 while(in < end) { | |
228 const int k = s->cyc_pos; | |
229 | |
230 update_ch(s, in, k); | |
231 | |
232 /* Simulate a 7.5 ms -20 dB echo of the center channel in the | |
233 front channels (like reflection from a room wall) - a kind of | |
234 psycho-acoustically "cheating" to focus the center front | |
235 channel, which is normally hard to be perceived as front */ | |
236 s->lf[k] += CFECHOAMPL * s->cf[(k + CFECHODELAY) % s->dlbuflen]; | |
237 s->rf[k] += CFECHOAMPL * s->cf[(k + CFECHODELAY) % s->dlbuflen]; | |
238 | |
239 /* Mixer filter matrix */ | |
240 common = conv(dblen, hlen, s->cf, s->cf_ir, k + s->cf_o); | |
241 if(s->matrix_mode) { | |
242 /* In matrix decoding mode, the rear channel gain must be | |
243 renormalized, as there is an additional channel. */ | |
244 matrix_decode_cr(s, in, k); | |
245 common += | |
246 conv(dblen, hlen, s->cr, s->cr_ir, k + s->cr_o) * | |
247 M1_76DB; | |
248 left = | |
249 ( conv(dblen, hlen, s->lf, s->af_ir, k + s->af_o) + | |
250 conv(dblen, hlen, s->rf, s->of_ir, k + s->of_o) + | |
251 (conv(dblen, hlen, s->lr, s->ar_ir, k + s->ar_o) + | |
252 conv(dblen, hlen, s->rr, s->or_ir, k + s->or_o)) * | |
253 M1_76DB + common); | |
254 right = | |
255 ( conv(dblen, hlen, s->rf, s->af_ir, k + s->af_o) + | |
256 conv(dblen, hlen, s->lf, s->of_ir, k + s->of_o) + | |
257 (conv(dblen, hlen, s->rr, s->ar_ir, k + s->ar_o) + | |
258 conv(dblen, hlen, s->lr, s->or_ir, k + s->or_o)) * | |
259 M1_76DB + common); | |
260 } | |
261 else { | |
262 left = | |
263 ( conv(dblen, hlen, s->lf, s->af_ir, k + s->af_o) + | |
264 conv(dblen, hlen, s->rf, s->of_ir, k + s->of_o) + | |
265 conv(dblen, hlen, s->lr, s->ar_ir, k + s->ar_o) + | |
266 conv(dblen, hlen, s->rr, s->or_ir, k + s->or_o) + | |
267 common); | |
268 right = | |
269 ( conv(dblen, hlen, s->rf, s->af_ir, k + s->af_o) + | |
270 conv(dblen, hlen, s->lf, s->of_ir, k + s->of_o) + | |
271 conv(dblen, hlen, s->rr, s->ar_ir, k + s->ar_o) + | |
272 conv(dblen, hlen, s->lr, s->or_ir, k + s->or_o) + | |
273 common); | |
274 } | |
275 | |
276 /* Bass compensation for the lower frequency cut of the HRTF. A | |
277 cross talk of the left and right channel is introduced to | |
278 match the directional characteristics of higher frequencies. | |
279 The bass will not have any real 3D perception, but that is | |
280 OK. */ | |
281 left_b = conv(dblen, blen, s->ba_l, s->ba_ir, k); | |
282 right_b = conv(dblen, blen, s->ba_r, s->ba_ir, k); | |
283 left += (1 - BASSCROSS) * left_b + BASSCROSS * right_b; | |
284 right += (1 - BASSCROSS) * right_b + BASSCROSS * left_b; | |
285 /* Also mix the LFE channel (if available) */ | |
286 if(af->data->nch >= 6) { | |
287 left += out[5] * M3_01DB; | |
288 right += out[5] * M3_01DB; | |
289 } | |
290 | |
291 /* Amplitude renormalization. */ | |
292 left *= AMPLNORM; | |
293 right *= AMPLNORM; | |
294 | |
295 /* "Cheating": linear stereo expansion to amplify the 3D | |
296 perception. Note: Too much will destroy the acoustic space | |
297 and may even result in headaches. */ | |
298 diff = STEXPAND2 * (left - right); | |
299 out[0] = (int16_t)(left + diff); | |
300 out[1] = (int16_t)(right - diff); | |
301 | |
302 /* The remaining channels are not needed any more */ | |
303 out[2] = out[3] = out[4] = 0; | |
304 if(af->data->nch >= 6) | |
305 out[5] = 0; | |
306 | |
307 /* Next sample... */ | |
308 in = &in[data->nch]; | |
309 out = &out[af->data->nch]; | |
310 (s->cyc_pos)--; | |
311 if(s->cyc_pos < 0) | |
312 s->cyc_pos += dblen; | |
313 } | |
314 | |
315 /* Set output data */ | |
316 data->audio = af->data->audio; | |
317 data->len = (data->len * af->mul.n) / af->mul.d; | |
318 data->nch = af->data->nch; | |
319 | |
320 return data; | |
321 } | |
322 | |
323 static int allocate(af_hrtf_t *s) | |
324 { | |
325 if ((s->lf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
326 if ((s->rf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
327 if ((s->lr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
328 if ((s->rr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
329 if ((s->cf = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
330 if ((s->cr = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
331 if ((s->ba_l = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
332 if ((s->ba_r = malloc(s->dlbuflen * sizeof(float))) == NULL) return -1; | |
333 return 0; | |
334 } | |
335 | |
336 /* Allocate memory and set function pointers */ | |
337 static int open(af_instance_t* af) | |
338 { | |
339 int i; | |
340 af_hrtf_t *s; | |
341 float fc; | |
342 | |
343 af->control = control; | |
344 af->uninit = uninit; | |
345 af->play = play; | |
346 af->mul.n = 1; | |
347 af->mul.d = 1; | |
348 af->data = calloc(1, sizeof(af_data_t)); | |
349 af->setup = calloc(1, sizeof(af_hrtf_t)); | |
350 if((af->data == NULL) || (af->setup == NULL)) | |
351 return AF_ERROR; | |
352 | |
353 s = af->setup; | |
354 | |
355 s->dlbuflen = DELAYBUFLEN; | |
356 s->hrflen = HRTFFILTLEN; | |
357 s->basslen = BASSFILTLEN; | |
358 | |
359 s->cyc_pos = s->dlbuflen - 1; | |
360 s->matrix_mode = 1; | |
361 | |
362 if (allocate(s) != 0) { | |
363 af_msg(AF_MSG_ERROR, "[hrtf] Memory allocation error.\n"); | |
364 return AF_ERROR; | |
365 } | |
366 | |
367 for(i = 0; i < s->dlbuflen; i++) | |
368 s->lf[i] = s->rf[i] = s->lr[i] = s->rr[i] = s->cf[i] = | |
369 s->cr[i] = 0; | |
370 | |
371 s->lr_fwr = | |
372 s->rr_fwr = 0; | |
373 | |
374 s->cf_ir = cf_filt + (s->cf_o = pulse_detect(cf_filt)); | |
375 s->af_ir = af_filt + (s->af_o = pulse_detect(af_filt)); | |
376 s->of_ir = of_filt + (s->of_o = pulse_detect(of_filt)); | |
377 s->ar_ir = ar_filt + (s->ar_o = pulse_detect(ar_filt)); | |
378 s->or_ir = or_filt + (s->or_o = pulse_detect(or_filt)); | |
379 s->cr_ir = cr_filt + (s->cr_o = pulse_detect(cr_filt)); | |
380 | |
381 if((s->ba_ir = malloc(s->basslen * sizeof(float))) == NULL) { | |
382 af_msg(AF_MSG_ERROR, "[hrtf] Memory allocation error.\n"); | |
383 return AF_ERROR; | |
384 } | |
385 fc = 2.0 * BASSFILTFREQ / (float)af->data->rate; | |
14275
de13fd557440
less namespace pollution #2 (prefixed globals in filter.c with af_filter_)
alex
parents:
14245
diff
changeset
|
386 if(af_filter_design_fir(s->basslen, s->ba_ir, &fc, LP | KAISER, 4 * M_PI) == |
13996 | 387 -1) { |
388 af_msg(AF_MSG_ERROR, "[hrtf] Unable to design low-pass " | |
389 "filter.\n"); | |
390 return AF_ERROR; | |
391 } | |
392 for(i = 0; i < s->basslen; i++) | |
393 s->ba_ir[i] *= BASSGAIN; | |
394 | |
395 return AF_OK; | |
396 } | |
397 | |
398 /* Description of this filter */ | |
399 af_info_t af_info_hrtf = { | |
400 "HRTF Headphone", | |
401 "hrtf", | |
402 "ylai", | |
403 "", | |
404 AF_FLAGS_REENTRANT, | |
405 open | |
406 }; |