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