1 /*
2 reverb.c
3
4 Midi Wavetable Processing library
5
6 Copyright (C) Chris Ison 2001-2011
7 Copyright (C) Bret Curtis 2013-2014
8
9 This file is part of WildMIDI.
10
11 WildMIDI is free software: you can redistribute and/or modify the player
12 under the terms of the GNU General Public License and you can redistribute
13 and/or modify the library under the terms of the GNU Lesser General Public
14 License as published by the Free Software Foundation, either version 3 of
15 the licenses, or(at your option) any later version.
16
17 WildMIDI is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License and
20 the GNU Lesser General Public License for more details.
21
22 You should have received a copy of the GNU General Public License and the
23 GNU Lesser General Public License along with WildMIDI. If not, see
24 <http://www.gnu.org/licenses/>.
25 */
26
27 //#include "config.h"
28
29 #include <math.h>
30 #include <stdlib.h>
31
32 #include "common.h"
33 #include "reverb.h"
34
35 namespace WildMidi
36 {
37 /*
38 reverb function
39 */
_WM_reset_reverb(struct _rvb * rvb)40 void _WM_reset_reverb(struct _rvb *rvb) {
41 int i, j, k;
42 for (i = 0; i < rvb->l_buf_size; i++) {
43 rvb->l_buf[i] = 0;
44 }
45 for (i = 0; i < rvb->r_buf_size; i++) {
46 rvb->r_buf[i] = 0;
47 }
48 for (k = 0; k < 8; k++) {
49 for (i = 0; i < 6; i++) {
50 for (j = 0; j < 2; j++) {
51 rvb->l_buf_flt_in[k][i][j] = 0;
52 rvb->l_buf_flt_out[k][i][j] = 0;
53 rvb->r_buf_flt_in[k][i][j] = 0;
54 rvb->r_buf_flt_out[k][i][j] = 0;
55 }
56 }
57 }
58 }
59
60 /*
61 _WM_init_reverb
62
63 =========================
64 Engine Description
65
66 8 reflective points around the room
67 2 speaker positions
68 1 listener position
69
70 Sounds come from the speakers to all points and to the listener.
71 Sound comes from the reflective points to the listener.
72 These sounds are combined, put through a filter that mimics surface absorbtion.
73 The combined sounds are also sent to the reflective points on the opposite side.
74
75 */
76 struct _rvb *
_WM_init_reverb(int rate,float room_x,float room_y,float listen_x,float listen_y)77 _WM_init_reverb(int rate, float room_x, float room_y, float listen_x,
78 float listen_y) {
79
80 /* filters set at 125Hz, 250Hz, 500Hz, 1000Hz, 2000Hz, 4000Hz */
81 double Freq[] = {125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0};
82
83 /* numbers calculated from
84 * 101.325 kPa, 20 deg C, 50% relative humidity */
85 double dbAirAbs[] = {-0.00044, -0.00131, -0.002728, -0.004665, -0.009887, -0.029665};
86
87 /* modify these to adjust the absorption qualities of the surface.
88 * Remember that lower frequencies are less effected by surfaces
89 * Note: I am currently playing with the values and finding the ideal surfaces
90 * for nice default reverb.
91 */
92 double dbAttn[8][6] = {
93 {-1.839, -6.205, -8.891, -12.059, -15.935, -20.942},
94 {-0.131, -6.205, -12.059, -20.933, -20.933, -15.944},
95 {-0.131, -6.205, -12.059, -20.933, -20.933, -15.944},
96 {-1.839, -6.205, -8.891, -12.059, -15.935, -20.942},
97 {-1.839, -6.205, -8.891, -12.059, -15.935, -20.942},
98 {-0.131, -6.205, -12.059, -20.933, -20.933, -15.944},
99 {-0.131, -6.205, -12.059, -20.933, -20.933, -15.944},
100 {-1.839, -6.205, -8.891, -12.059, -15.935, -20.942}
101 };
102 /*
103 double dbAttn[6] = {
104 // concrete covered in carpet
105 // -0.175, -0.537, -1.412, -4.437, -7.959, -7.959
106 // pleated drapes
107 -0.630, -3.223, -5.849, -12.041, -10.458, -7.959
108 };
109 */
110
111 /* distance */
112 double SPL_DST[8] = {0.0};
113 double SPR_DST[8] = {0.0};
114 double RFN_DST[8] = {0.0};
115
116 double MAXL_DST = 0.0;
117 double MAXR_DST = 0.0;
118
119 double SPL_LSN_XOFS = 0.0;
120 double SPL_LSN_YOFS = 0.0;
121 double SPL_LSN_DST = 0.0;
122
123 double SPR_LSN_XOFS = 0.0;
124 double SPR_LSN_YOFS = 0.0;
125 double SPR_LSN_DST = 0.0;
126
127
128 struct _rvb *rtn_rvb = (struct _rvb*)malloc(sizeof(struct _rvb));
129 int j = 0;
130 int i = 0;
131
132 struct _coord {
133 double x;
134 double y;
135 };
136
137 #if 0
138 struct _coord SPL = {2.5, 5.0}; /* Left Speaker Position */
139 struct _coord SPR = {7.5, 5.0}; /* Right Speaker Position */
140 /* position of the reflective points */
141 struct _coord RFN[] = {
142 { 5.0, 0.0},
143 { 0.0, 6.66666},
144 { 0.0, 13.3333},
145 { 5.0, 20.0},
146 { 10.0, 20.0},
147 { 15.0, 13.3333},
148 { 15.0, 6.66666},
149 { 10.0, 0.0}
150 };
151 #else
152 struct _coord SPL; /* Left Speaker Position */
153 struct _coord SPR; /* Right Speaker Position */
154 /* position of the reflective points */
155 struct _coord RFN[8];
156
157 SPL.x = room_x / 4.0;
158 SPR.x = room_x / 4.0 * 3.0;
159 SPL.y = room_y / 10.0;
160 SPR.y = room_y / 10.0;
161
162 RFN[0].x = room_x / 3.0;
163 RFN[0].y = 0.0;
164 RFN[1].x = 0.0;
165 RFN[1].y = room_y / 3.0;
166 RFN[2].x = 0.0;
167 RFN[2].y = room_y / 3.0 * 2.0;
168 RFN[3].x = room_x / 3.0;
169 RFN[3].y = room_y;
170 RFN[4].x = room_x / 3.0 * 2.0;
171 RFN[4].y = room_y;
172 RFN[5].x = room_x;
173 RFN[5].y = room_y / 3.0 * 2.0;
174 RFN[6].x = room_x;
175 RFN[6].y = room_y / 3.0;
176 RFN[7].x = room_x / 3.0 * 2.0;
177 RFN[7].y = 0.0;
178 #endif
179
180 SPL_LSN_XOFS = SPL.x - listen_x;
181 SPL_LSN_YOFS = SPL.y - listen_y;
182 SPL_LSN_DST = sqrt((SPL_LSN_XOFS * SPL_LSN_XOFS) + (SPL_LSN_YOFS * SPL_LSN_YOFS));
183
184 if (SPL_LSN_DST > MAXL_DST)
185 MAXL_DST = SPL_LSN_DST;
186
187 SPR_LSN_XOFS = SPR.x - listen_x;
188 SPR_LSN_YOFS = SPR.y - listen_y;
189 SPR_LSN_DST = sqrt((SPR_LSN_XOFS * SPR_LSN_XOFS) + (SPR_LSN_YOFS * SPR_LSN_YOFS));
190
191 if (SPR_LSN_DST > MAXR_DST)
192 MAXR_DST = SPR_LSN_DST;
193
194 if (rtn_rvb == NULL) {
195 return NULL;
196 }
197
198 for (j = 0; j < 8; j++) {
199 double SPL_RFL_XOFS = 0;
200 double SPL_RFL_YOFS = 0;
201 double SPR_RFL_XOFS = 0;
202 double SPR_RFL_YOFS = 0;
203 double RFN_XOFS = listen_x - RFN[j].x;
204 double RFN_YOFS = listen_y - RFN[j].y;
205 RFN_DST[j] = sqrt((RFN_XOFS * RFN_XOFS) + (RFN_YOFS * RFN_YOFS));
206
207 SPL_RFL_XOFS = SPL.x - RFN[i].x;
208 SPL_RFL_YOFS = SPL.y - RFN[i].y;
209 SPR_RFL_XOFS = SPR.x - RFN[i].x;
210 SPR_RFL_YOFS = SPR.y - RFN[i].y;
211 SPL_DST[i] = sqrt(
212 (SPL_RFL_XOFS * SPL_RFL_XOFS) + (SPL_RFL_YOFS * SPL_RFL_YOFS));
213 SPR_DST[i] = sqrt(
214 (SPR_RFL_XOFS * SPR_RFL_XOFS) + (SPR_RFL_YOFS * SPR_RFL_YOFS));
215 /*
216 add the 2 distances together and remove the speaker to listener distance
217 so we dont have to delay the initial output
218 */
219 SPL_DST[i] += RFN_DST[i];
220
221 /* so i dont have to delay speaker output */
222 SPL_DST[i] -= SPL_LSN_DST;
223
224 if (i < 4) {
225 if (SPL_DST[i] > MAXL_DST)
226 MAXL_DST = SPL_DST[i];
227 } else {
228 if (SPL_DST[i] > MAXR_DST)
229 MAXR_DST = SPL_DST[i];
230 }
231
232 SPR_DST[i] += RFN_DST[i];
233
234 /* so i dont have to delay speaker output */
235 SPR_DST[i] -= SPR_LSN_DST;
236
237 if (i < 4) {
238 if (SPR_DST[i] > MAXL_DST)
239 MAXL_DST = SPR_DST[i];
240 } else {
241 if (SPR_DST[i] > MAXR_DST)
242 MAXR_DST = SPR_DST[i];
243 }
244
245 RFN_DST[j] *= 2.0;
246
247 if (j < 4) {
248 if (RFN_DST[j] > MAXL_DST)
249 MAXL_DST = RFN_DST[j];
250 } else {
251 if (RFN_DST[j] > MAXR_DST)
252 MAXR_DST = RFN_DST[j];
253 }
254
255 for (i = 0; i < 6; i++) {
256 double srate = (double) rate;
257 double bandwidth = 2.0;
258 double omega = 2.0 * M_PI * Freq[i] / srate;
259 double sn = sin(omega);
260 double cs = cos(omega);
261 double alpha = sn * sinh(M_LN2 / 2 * bandwidth * omega / sn);
262 double A = pow(10.0, ((/*dbAttn[i]*/dbAttn[j][i] +
263 (dbAirAbs[i] * RFN_DST[j])) / 40.0) );
264 /*
265 Peaking band EQ filter
266 */
267 double b0 = 1 + (alpha * A);
268 double b1 = -2 * cs;
269 double b2 = 1 - (alpha * A);
270 double a0 = 1 + (alpha / A);
271 double a1 = -2 * cs;
272 double a2 = 1 - (alpha / A);
273
274 rtn_rvb->coeff[j][i][0] = (signed int) ((b0 / a0) * 1024.0);
275 rtn_rvb->coeff[j][i][1] = (signed int) ((b1 / a0) * 1024.0);
276 rtn_rvb->coeff[j][i][2] = (signed int) ((b2 / a0) * 1024.0);
277 rtn_rvb->coeff[j][i][3] = (signed int) ((a1 / a0) * 1024.0);
278 rtn_rvb->coeff[j][i][4] = (signed int) ((a2 / a0) * 1024.0);
279 }
280 }
281
282 /* init the reverb buffers */
283 rtn_rvb->l_buf_size = (int) ((float) rate * (MAXL_DST / 340.29));
284 rtn_rvb->l_buf = (int*)malloc(
285 sizeof(signed int) * (rtn_rvb->l_buf_size + 1));
286 rtn_rvb->l_out = 0;
287
288 rtn_rvb->r_buf_size = (int) ((float) rate * (MAXR_DST / 340.29));
289 rtn_rvb->r_buf = (int*)malloc(
290 sizeof(signed int) * (rtn_rvb->r_buf_size + 1));
291 rtn_rvb->r_out = 0;
292
293 for (i = 0; i < 4; i++) {
294 rtn_rvb->l_sp_in[i] = (int) ((float) rate * (SPL_DST[i] / 340.29));
295 rtn_rvb->l_sp_in[i + 4] = (int) ((float) rate
296 * (SPL_DST[i + 4] / 340.29));
297 rtn_rvb->r_sp_in[i] = (int) ((float) rate * (SPR_DST[i] / 340.29));
298 rtn_rvb->r_sp_in[i + 4] = (int) ((float) rate
299 * (SPR_DST[i + 4] / 340.29));
300 rtn_rvb->l_in[i] = (int) ((float) rate * (RFN_DST[i] / 340.29));
301 rtn_rvb->r_in[i] = (int) ((float) rate * (RFN_DST[i + 4] / 340.29));
302 }
303
304 rtn_rvb->gain = 4;
305
306 _WM_reset_reverb(rtn_rvb);
307 return rtn_rvb;
308 }
309
310 /* _WM_free_reverb - free up memory used for reverb */
_WM_free_reverb(struct _rvb * rvb)311 void _WM_free_reverb(struct _rvb *rvb) {
312 if (!rvb) return;
313 free(rvb->l_buf);
314 free(rvb->r_buf);
315 free(rvb);
316 }
317
_WM_do_reverb(struct _rvb * rvb,signed int * buffer,int size)318 void _WM_do_reverb(struct _rvb *rvb, signed int *buffer, int size) {
319 int i, j, k;
320 signed int l_buf_flt = 0;
321 signed int r_buf_flt = 0;
322 signed int l_rfl = 0;
323 signed int r_rfl = 0;
324 int vol_div = 64;
325
326 for (i = 0; i < size; i += 2) {
327 signed int tmp_l_val = 0;
328 signed int tmp_r_val = 0;
329 /*
330 add the initial reflections
331 from each speaker, 4 to go the left, 4 go to the right buffers
332 */
333 tmp_l_val = buffer[i] / vol_div;
334 tmp_r_val = buffer[i + 1] / vol_div;
335 for (j = 0; j < 4; j++) {
336 rvb->l_buf[rvb->l_sp_in[j]] += tmp_l_val;
337 rvb->l_sp_in[j] = (rvb->l_sp_in[j] + 1) % rvb->l_buf_size;
338 rvb->l_buf[rvb->r_sp_in[j]] += tmp_r_val;
339 rvb->r_sp_in[j] = (rvb->r_sp_in[j] + 1) % rvb->l_buf_size;
340
341 rvb->r_buf[rvb->l_sp_in[j + 4]] += tmp_l_val;
342 rvb->l_sp_in[j + 4] = (rvb->l_sp_in[j + 4] + 1) % rvb->r_buf_size;
343 rvb->r_buf[rvb->r_sp_in[j + 4]] += tmp_r_val;
344 rvb->r_sp_in[j + 4] = (rvb->r_sp_in[j + 4] + 1) % rvb->r_buf_size;
345 }
346
347 /*
348 filter the reverb output and add to buffer
349 */
350 l_rfl = rvb->l_buf[rvb->l_out];
351 rvb->l_buf[rvb->l_out] = 0;
352 rvb->l_out = (rvb->l_out + 1) % rvb->l_buf_size;
353
354 r_rfl = rvb->r_buf[rvb->r_out];
355 rvb->r_buf[rvb->r_out] = 0;
356 rvb->r_out = (rvb->r_out + 1) % rvb->r_buf_size;
357
358 for (k = 0; k < 8; k++) {
359 for (j = 0; j < 6; j++) {
360 l_buf_flt = ((l_rfl * rvb->coeff[k][j][0])
361 + (rvb->l_buf_flt_in[k][j][0] * rvb->coeff[k][j][1])
362 + (rvb->l_buf_flt_in[k][j][1] * rvb->coeff[k][j][2])
363 - (rvb->l_buf_flt_out[k][j][0] * rvb->coeff[k][j][3])
364 - (rvb->l_buf_flt_out[k][j][1] * rvb->coeff[k][j][4]))
365 / 1024;
366 rvb->l_buf_flt_in[k][j][1] = rvb->l_buf_flt_in[k][j][0];
367 rvb->l_buf_flt_in[k][j][0] = l_rfl;
368 rvb->l_buf_flt_out[k][j][1] = rvb->l_buf_flt_out[k][j][0];
369 rvb->l_buf_flt_out[k][j][0] = l_buf_flt;
370 buffer[i] += l_buf_flt / 8;
371
372 r_buf_flt = ((r_rfl * rvb->coeff[k][j][0])
373 + (rvb->r_buf_flt_in[k][j][0] * rvb->coeff[k][j][1])
374 + (rvb->r_buf_flt_in[k][j][1] * rvb->coeff[k][j][2])
375 - (rvb->r_buf_flt_out[k][j][0] * rvb->coeff[k][j][3])
376 - (rvb->r_buf_flt_out[k][j][1] * rvb->coeff[k][j][4]))
377 / 1024;
378 rvb->r_buf_flt_in[k][j][1] = rvb->r_buf_flt_in[k][j][0];
379 rvb->r_buf_flt_in[k][j][0] = r_rfl;
380 rvb->r_buf_flt_out[k][j][1] = rvb->r_buf_flt_out[k][j][0];
381 rvb->r_buf_flt_out[k][j][0] = r_buf_flt;
382 buffer[i + 1] += r_buf_flt / 8;
383 }
384 }
385
386 /*
387 add filtered result back into the buffers but on the opposite side
388 */
389 tmp_l_val = buffer[i + 1] / vol_div;
390 tmp_r_val = buffer[i] / vol_div;
391 for (j = 0; j < 4; j++) {
392 rvb->l_buf[rvb->l_in[j]] += tmp_l_val;
393 rvb->l_in[j] = (rvb->l_in[j] + 1) % rvb->l_buf_size;
394
395 rvb->r_buf[rvb->r_in[j]] += tmp_r_val;
396 rvb->r_in[j] = (rvb->r_in[j] + 1) % rvb->r_buf_size;
397 }
398 }
399 }
400
401 }
402