1 /*
2 afilters.c:
3
4 Copyright (C) 1991 Barry Vercoe, John ffitch, Gabriel Maldonado
5
6 This file is part of Csound.
7
8 The Csound Library is free software; you can redistribute it
9 and/or modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12
13 Csound is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with Csound; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21 02110-1301 USA
22 */
23
24 #include "csoundCore.h" /* AFILTERS.C */
25 #include "ugens5.h"
26 #include <math.h>
27
28 extern int32_t rsnset(CSOUND *csound, RESON *p);
29
aresonaa(CSOUND * csound,RESON * p)30 static int32_t aresonaa(CSOUND *csound, RESON *p)
31 {
32 uint32_t offset = p->h.insdshead->ksmps_offset;
33 uint32_t early = p->h.insdshead->ksmps_no_end;
34 uint32_t flag = 0, n, nsmps = CS_KSMPS;
35 MYFLT *ar, *asig;
36 double c3p1, c3t4, omc3, c2sqr; /* 1/RMS = root2 (rand) */
37 /* or 1/.5 (sine) */
38 double yt1, yt2, c1, c2, c3;
39
40 asig = p->asig;
41 ar = p->ar;
42 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
43 if (UNLIKELY(early)) {
44 nsmps -= early;
45 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
46 }
47 c1 = p->c1; c2 = p->c2; c3 = p->c3; yt1 = p->yt1; yt2 = p->yt2;
48 if (p->scale == 1 || p->scale == 0) {
49 for (n=offset; n<nsmps; n++) {
50 double sig = (double)asig[n];
51 double ans;
52 if (p->kcf[n] != (MYFLT)p->prvcf) {
53 p->prvcf = (double)p->kcf[n];
54 p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
55 flag = 1;
56 }
57 if (p->kbw[n] != (MYFLT)p->prvbw) {
58 p->prvbw = (double)p->kbw[n];
59 p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
60 flag = 1;
61 }
62 if (flag) {
63 c3p1 = p->c3 + 1.0;
64 c3t4 = p->c3 * 4.0;
65 omc3 = 1.0 - p->c3;
66 p->c2 = c3t4 * p->cosf / c3p1;
67 c2sqr = p->c2 * p->c2;
68 if (p->scale == 1) /* i.e. 1 - A(reson) */
69 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
70 else if (p->scale == 2) /* i.e. D - A(reson) */
71 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
72 else p->c1 = 0.0; /* cannot tell */
73 }
74 ans = c1 * sig + c2 * yt1 - c3 * yt2;
75 yt2 = yt1;
76 yt1 = ans - sig; /* yt1 contains yt1-xt1 */
77 ar[n] = (MYFLT)ans;
78 }
79 }
80 else if (p->scale == 2) {
81 for (n=offset; n<nsmps; n++) {
82 double sig = (double)asig[n];
83 double ans;
84 if (p->kcf[n] != (MYFLT)p->prvcf) {
85 p->prvcf = (double)p->kcf[n];
86 p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
87 flag = 1;
88 }
89 if (p->kbw[n] != (MYFLT)p->prvbw) {
90 p->prvbw = (double)p->kbw[n];
91 p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
92 flag = 1;
93 }
94 if (flag) {
95 c3p1 = p->c3 + 1.0;
96 c3t4 = p->c3 * 4.0;
97 omc3 = 1.0 - p->c3;
98 p->c2 = c3t4 * p->cosf / c3p1;
99 c2sqr = p->c2 * p->c2;
100 if (p->scale == 1) /* i.e. 1 - A(reson) */
101 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
102 else if (p->scale == 2) /* i.e. D - A(reson) */
103 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
104 else p->c1 = 0.0; /* cannot tell */
105 }
106 ans = c1 * sig + c2 * yt1 - c3 * yt2;
107 yt2 = yt1;
108 yt1 = ans - 2.0 * sig; /* yt1 contains yt1-D*xt1 */
109 ar[n] = (MYFLT)ans;
110 }
111 }
112 p->yt1 = yt1; p->yt2 = yt2;
113 return OK;
114 }
115
aresonak(CSOUND * csound,RESON * p)116 static int32_t aresonak(CSOUND *csound, RESON *p)
117 {
118 uint32_t offset = p->h.insdshead->ksmps_offset;
119 uint32_t early = p->h.insdshead->ksmps_no_end;
120 uint32_t n, nsmps = CS_KSMPS;
121 MYFLT *ar, *asig;
122 double c3p1, c3t4, omc3, c2sqr; /* 1/RMS = root2 (rand) */
123 /* or 1/.5 (sine) */
124 double yt1, yt2, c1, c2, c3;
125
126 if (*p->kbw != (MYFLT)p->prvbw) {
127 p->prvbw = (double)*p->kbw;
128 p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
129 c3p1 = p->c3 + 1.0;
130 c3t4 = p->c3 * 4.0;
131 omc3 = 1.0 - p->c3;
132 p->c2 = c3t4 * p->cosf / c3p1;
133 c2sqr = p->c2 * p->c2;
134 if (p->scale == 1) /* i.e. 1 - A(reson) */
135 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
136 else if (p->scale == 2) /* i.e. D - A(reson) */
137 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
138 else p->c1 = 0.0; /* cannot tell */
139 }
140 asig = p->asig;
141 ar = p->ar;
142 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
143 if (UNLIKELY(early)) {
144 nsmps -= early;
145 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
146 }
147 c1 = p->c1; c2 = p->c2; c3 = p->c3; yt1 = p->yt1; yt2 = p->yt2;
148 if (p->scale == 1 || p->scale == 0) {
149 for (n=offset; n<nsmps; n++) {
150 double sig = (double)asig[n];
151 double ans;
152 if (p->kcf[n] != (MYFLT)p->prvcf) {
153 p->prvcf = (double)p->kcf[n];
154 p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
155 c3p1 = p->c3 + 1.0;
156 c3t4 = p->c3 * 4.0;
157 omc3 = 1.0 - p->c3;
158 p->c2 = c3t4 * p->cosf / c3p1;
159 c2sqr = p->c2 * p->c2;
160 if (p->scale == 1) /* i.e. 1 - A(reson) */
161 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
162 else if (p->scale == 2) /* i.e. D - A(reson) */
163 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
164 else p->c1 = 0.0; /* cannot tell */
165 }
166 ans = c1 * sig + c2 * yt1 - c3 * yt2;
167 yt2 = yt1;
168 yt1 = ans - sig; /* yt1 contains yt1-xt1 */
169 ar[n] = (MYFLT)ans;
170 }
171 }
172 else if (p->scale == 2) {
173 for (n=offset; n<nsmps; n++) {
174 double sig = (double)asig[n];
175 double ans;
176 if (p->kcf[n] != (MYFLT)p->prvcf) {
177 p->prvcf = (double)p->kcf[n];
178 p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
179 c3p1 = p->c3 + 1.0;
180 c3t4 = p->c3 * 4.0;
181 omc3 = 1.0 - p->c3;
182 p->c2 = c3t4 * p->cosf / c3p1;
183 c2sqr = p->c2 * p->c2;
184 if (p->scale == 1) /* i.e. 1 - A(reson) */
185 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
186 else if (p->scale == 2) /* i.e. D - A(reson) */
187 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
188 else p->c1 = 0.0; /* cannot tell */
189 }
190 ans = c1 * sig + c2 * yt1 - c3 * yt2;
191 yt2 = yt1;
192 yt1 = ans - 2.0 * sig; /* yt1 contains yt1-D*xt1 */
193 ar[n] = (MYFLT)ans;
194 }
195 }
196 p->yt1 = yt1; p->yt2 = yt2;
197 return OK;
198 }
199
aresonka(CSOUND * csound,RESON * p)200 static int32_t aresonka(CSOUND *csound, RESON *p)
201 {
202 uint32_t offset = p->h.insdshead->ksmps_offset;
203 uint32_t early = p->h.insdshead->ksmps_no_end;
204 uint32_t n, nsmps = CS_KSMPS;
205 MYFLT *ar, *asig;
206 double c3p1, c3t4, omc3, c2sqr; /* 1/RMS = root2 (rand) */
207 /* or 1/.5 (sine) */
208 double yt1, yt2, c1, c2, c3;
209
210 if (*p->kcf != (MYFLT)p->prvcf) {
211 p->prvcf = (double)*p->kcf;
212 p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
213 c3p1 = p->c3 + 1.0;
214 c3t4 = p->c3 * 4.0;
215 omc3 = 1.0 - p->c3;
216 p->c2 = c3t4 * p->cosf / c3p1;
217 c2sqr = p->c2 * p->c2;
218 if (p->scale == 1) /* i.e. 1 - A(reson) */
219 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
220 else if (p->scale == 2) /* i.e. D - A(reson) */
221 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
222 else p->c1 = 0.0; /* cannot tell */
223 }
224 asig = p->asig;
225 ar = p->ar;
226 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
227 if (UNLIKELY(early)) {
228 nsmps -= early;
229 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
230 }
231 c1 = p->c1; c2 = p->c2; c3 = p->c3; yt1 = p->yt1; yt2 = p->yt2;
232 if (p->scale == 1 || p->scale == 0) {
233 for (n=offset; n<nsmps; n++) {
234 double sig = (double)asig[n];
235 double ans;
236 if (p->kbw[n] != (MYFLT)p->prvbw) {
237 p->prvbw = (double)p->kbw[n];
238 p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
239 c3p1 = p->c3 + 1.0;
240 c3t4 = p->c3 * 4.0;
241 omc3 = 1.0 - p->c3;
242 p->c2 = c3t4 * p->cosf / c3p1;
243 c2sqr = p->c2 * p->c2;
244 if (p->scale == 1) /* i.e. 1 - A(reson) */
245 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
246 else if (p->scale == 2) /* i.e. D - A(reson) */
247 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
248 else p->c1 = 0.0; /* cannot tell */
249 }
250 ans = c1 * sig + c2 * yt1 - c3 * yt2;
251 yt2 = yt1;
252 yt1 = ans - sig; /* yt1 contains yt1-xt1 */
253 ar[n] = (MYFLT)ans;
254 }
255 }
256 else if (p->scale == 2) {
257 for (n=offset; n<nsmps; n++) {
258 double sig = (double)asig[n];
259 double ans;
260 if (p->kbw[n] != (MYFLT)p->prvbw) {
261 p->prvbw = (double)p->kbw[n];
262 p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
263 c3p1 = p->c3 + 1.0;
264 c3t4 = p->c3 * 4.0;
265 omc3 = 1.0 - p->c3;
266 p->c2 = c3t4 * p->cosf / c3p1;
267 c2sqr = p->c2 * p->c2;
268 if (p->scale == 1) /* i.e. 1 - A(reson) */
269 p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
270 else if (p->scale == 2) /* i.e. D - A(reson) */
271 p->c1 = 2.0 - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
272 else p->c1 = 0.0; /* cannot tell */
273 }
274 ans = c1 * sig + c2 * yt1 - c3 * yt2;
275 yt2 = yt1;
276 yt1 = ans - 2.0 * sig; /* yt1 contains yt1-D*xt1 */
277 ar[n] = (MYFLT)ans;
278 }
279 }
280 p->yt1 = yt1; p->yt2 = yt2;
281 return OK;
282 }
283
284 extern int32_t tonset(CSOUND*, TONE*);
285 extern int32_t tonsetx(CSOUND *csound, TONEX *p);
286
atonea(CSOUND * csound,TONE * p)287 static int32_t atonea(CSOUND *csound, TONE *p)
288 {
289 MYFLT *ar, *asig;
290 uint32_t offset = p->h.insdshead->ksmps_offset;
291 uint32_t early = p->h.insdshead->ksmps_no_end;
292 uint32_t n, nsmps = CS_KSMPS;
293 double c2 = p->c2, yt1 = p->yt1;
294
295 ar = p->ar;
296 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
297 if (UNLIKELY(early)) {
298 nsmps -= early;
299 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
300 }
301 asig = p->asig;
302 for (n=offset; n<nsmps; n++) {
303 double sig = (double)asig[n];
304 double x;
305 if (p->khp[n] != p->prvhp) {
306 double b;
307 p->prvhp = p->khp[n];
308 b = 2.0 - cos((double)(p->khp[n] * csound->tpidsr));
309 p->c2 = c2 = b - sqrt(b * b - 1.0);
310 /* p->c1 = c1 = 1.0 - c2; */
311 }
312 x = yt1 = c2 * (yt1 + sig);
313 ar[n] = (MYFLT)x;
314 yt1 -= sig; /* yt1 contains yt1-xt1 */
315 }
316 p->yt1 = yt1;
317 return OK;
318 }
319
tonea(CSOUND * csound,TONE * p)320 static int32_t tonea(CSOUND *csound, TONE *p)
321 {
322 MYFLT *ar, *asig;
323 uint32_t offset = p->h.insdshead->ksmps_offset;
324 uint32_t early = p->h.insdshead->ksmps_no_end;
325 uint32_t n, nsmps = CS_KSMPS;
326 double c1 = p->c1, c2 = p->c2;
327 double yt1 = p->yt1;
328 double prvhp = p->prvhp;
329
330 ar = p->ar;
331 asig = p->asig;
332 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
333 if (UNLIKELY(early)) {
334 nsmps -= early;
335 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
336 }
337 for (n=offset; n<nsmps; n++) {
338 if (p->khp[n] != prvhp) {
339 double b;
340 prvhp = (double)p->khp[n];
341 b = 2.0 - cos((double)(prvhp * csound->tpidsr));
342 c2 = b - sqrt(b * b - 1.0);
343 c1 = 1.0 - c2;
344 }
345 yt1 = c1 * (double)(asig[n]) + c2 * yt1;
346 ar[n] = (MYFLT)yt1;
347 }
348 p->yt1 = yt1;
349 p->prvhp = prvhp;
350 p->c1 = c1;
351 p->c2 = c2;
352 return OK;
353 }
354
tonexa(CSOUND * csound,TONEX * p)355 static int32_t tonexa(CSOUND *csound, TONEX *p) /* From G Maldonado, modified */
356 {
357 MYFLT *ar = p->ar;
358 double c2 = p->c2, *yt1 = p->yt1,c1 = p->c1;
359 uint32_t offset = p->h.insdshead->ksmps_offset;
360 uint32_t early = p->h.insdshead->ksmps_no_end;
361 uint32_t n, nsmps = CS_KSMPS;
362 int32_t j, lp = p->loop;
363
364 memmove(ar,p->asig,sizeof(MYFLT)*nsmps);
365 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
366 if (UNLIKELY(early)) {
367 nsmps -= early;
368 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
369 }
370 for (j=0; j< lp; j++) {
371 /* Should *yt1 be reset to something?? */
372 for (n=offset; n<nsmps; n++) {
373 double x;
374 if (p->khp[n] != p->prvhp) {
375 double b;
376 p->prvhp = (double)p->khp[n];
377 b = 2.0 - cos(p->prvhp * (double)csound->tpidsr);
378 p->c2 = b - sqrt(b * b - 1.0);
379 p->c1 = 1.0 - p->c2;
380 }
381 x = c1 * ar[n] + c2 * yt1[j];
382 yt1[j] = x;
383 ar[n] = (MYFLT)x;
384 }
385 }
386 return OK;
387 }
388
atonexa(CSOUND * csound,TONEX * p)389 static int32_t atonexa(CSOUND *csound, TONEX *p) /* Gabriel Maldonado, modified */
390 {
391 MYFLT *ar = p->ar;
392 double c2 = p->c2, *yt1 = p->yt1;
393 uint32_t offset = p->h.insdshead->ksmps_offset;
394 uint32_t early = p->h.insdshead->ksmps_no_end;
395 uint32_t n, nsmps = CS_KSMPS;
396 int32_t j, lp = p->loop;
397 MYFLT prvhp = p->prvhp;
398
399 memmove(ar,p->asig,sizeof(MYFLT)*nsmps);
400 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
401 if (UNLIKELY(early)) {
402 nsmps -= early;
403 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
404 }
405 for (j=1; j<lp; j++) {
406 for (n=offset; n<nsmps; n++) {
407 double sig = (double)ar[n];
408 double x;
409 if (p->khp[n] != prvhp) {
410 double b;
411 prvhp = p->khp[n];
412 b = 2.0 - cos((double)(p->prvhp * csound->tpidsr));
413 c2 = b - sqrt(b * b - 1.0);
414 }
415 x = c2 * (yt1[j] + sig);
416 yt1[j] = x - sig; /* yt1 contains yt1-xt1 */
417 ar[n] = (MYFLT)x;
418 }
419 }
420 p->c2 = c2;
421 p->prvhp = prvhp;
422 return OK;
423 }
424
425
426 typedef struct {
427 OPDS h;
428 MYFLT *sr, *ain, *afc, *istor;
429 MYFLT lkf;
430 double a[8];
431 } BFIL;
432
433 typedef struct {
434 OPDS h;
435 MYFLT *sr, *ain, *kfo, *kbw, *istor;
436 MYFLT lkf, lkb;
437 double a[8];
438 } BBFIL;
439
440 //#define ROOT2 (1.4142135623730950488)
441
442 extern int32_t butset(CSOUND *csound, BFIL *p);
443
bbutset(CSOUND * csound,BBFIL * p)444 static int32_t bbutset(CSOUND *csound, BBFIL *p) /* Band set-up */
445 {
446 IGN(csound);
447 if (*p->istor==FL(0.0)) {
448 p->a[6] = p->a[7] = 0.0;
449 p->lkb = FL(0.0);
450 p->lkf = FL(0.0);
451 }
452 return OK;
453 }
454
455
hibuta(CSOUND * csound,BFIL * p)456 static int32_t hibuta(CSOUND *csound, BFIL *p) /* Hipass filter */
457 {
458 MYFLT *out, *in;
459 uint32_t offset = p->h.insdshead->ksmps_offset;
460 uint32_t early = p->h.insdshead->ksmps_no_end;
461 uint32_t nsmps = CS_KSMPS;
462 double *a;
463 double t, y;
464 uint32_t nn;
465 a = p->a;
466
467 in = p->ain;
468 out = p->sr;
469 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
470 if (UNLIKELY(early)) {
471 nsmps -= early;
472 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
473 }
474
475 if (UNLIKELY(p->afc[0] <= FL(0.0))) {
476 memcpy(&out[offset], &in[offset], (nsmps-offset)*sizeof(MYFLT));
477 return OK;
478 }
479
480 /* if (p->afc[0] != p->lkf) { */
481 /* p->lkf = p->afc[0]; */
482 /* c = tan((double)(csound->pidsr * p->lkf)); */
483
484 /* a[1] = 1.0 / ( 1.0 + ROOT2 * c + c * c); */
485 /* a[2] = -(a[1] + a[1]); */
486 /* a[3] = a[1]; */
487 /* a[4] = 2.0 * ( c*c - 1.0) * a[1]; */
488 /* a[5] = ( 1.0 - ROOT2 * c + c * c) * a[1]; */
489 /* } */
490 for (nn=offset; nn<nsmps; nn++) {
491 if (p->afc[nn] != p->lkf) {
492 double c;
493 p->lkf = p->afc[nn];
494 c = tan((double)(csound->pidsr * p->lkf));
495
496 a[1] = 1.0 / ( 1.0 + ROOT2 * c + c * c);
497 a[2] = -(a[1] + a[1]);
498 a[3] = a[1];
499 a[4] = 2.0 * ( c*c - 1.0) * a[1];
500 a[5] = ( 1.0 - ROOT2 * c + c * c) * a[1];
501 }
502 t = (double)in[nn] - a[4] * a[6] - a[5] * a[7];
503 t = csoundUndenormalizeDouble(t); /* Not needed on AMD */
504 y = t * a[1] + a[2] * a[6] + a[3] * a[7];
505 a[7] = a[6];
506 a[6] = t;
507 out[nn] = (MYFLT)y;
508 }
509 return OK;
510 }
511
lobuta(CSOUND * csound,BFIL * p)512 static int32_t lobuta(CSOUND *csound, BFIL *p) /* Lopass filter */
513 {
514 MYFLT *out, *in;
515 uint32_t offset = p->h.insdshead->ksmps_offset;
516 uint32_t early = p->h.insdshead->ksmps_no_end;
517 uint32_t nsmps = CS_KSMPS;
518 double *a = p->a;
519 double t, y;
520 uint32_t nn;
521
522 in = p->ain;
523 out = p->sr;
524
525 if (UNLIKELY(*p->afc <= FL(0.0))) {
526 memset(out, 0, CS_KSMPS*sizeof(MYFLT));
527 return OK;
528 }
529
530 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
531 if (UNLIKELY(early)) {
532 nsmps -= early;
533 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
534 }
535
536 /* if (p->afc[0] != p->lkf) { */
537 /* p->lkf = p->afc[0]; */
538 /* c = 1.0 / tan((double)(csound->pidsr * p->lkf)); */
539 /* a[1] = 1.0 / ( 1.0 + ROOT2 * c + c * c); */
540 /* a[2] = a[1] + a[1]; */
541 /* a[3] = a[1]; */
542 /* a[4] = 2.0 * ( 1.0 - c*c) * a[1]; */
543 /* a[5] = ( 1.0 - ROOT2 * c + c * c) * a[1]; */
544 /* } */
545
546 for (nn=offset; nn<nsmps; nn++) {
547 if (p->afc[nn] != p->lkf) {
548 double c;
549 p->lkf = p->afc[nn];
550 c = 1.0 / tan((double)(csound->pidsr * p->lkf));
551 a[1] = 1.0 / ( 1.0 + ROOT2 * c + c * c);
552 a[2] = a[1] + a[1];
553 a[3] = a[1];
554 a[4] = 2.0 * ( 1.0 - c*c) * a[1];
555 a[5] = ( 1.0 - ROOT2 * c + c * c) * a[1];
556 }
557 t = (double)in[nn] - a[4] * a[6] - a[5] * a[7];
558 t = csoundUndenormalizeDouble(t); /* Not needed on AMD */
559 y = t * a[1] + a[2] * a[6] + a[3] * a[7];
560 a[7] = a[6];
561 a[6] = t;
562 out[nn] = (MYFLT)y;
563 }
564 return OK;
565 }
566
bppasxx(CSOUND * csound,BBFIL * p)567 static int32_t bppasxx(CSOUND *csound, BBFIL *p) /* Bandpass filter */
568 {
569 uint32_t offset = p->h.insdshead->ksmps_offset;
570 uint32_t early = p->h.insdshead->ksmps_no_end;
571 uint32_t nsmps = CS_KSMPS;
572 MYFLT *out, *in;
573 double *a = p->a;
574 double t, y;
575 uint32_t nn;
576 int32_t asgbw = IS_ASIG_ARG(p->kbw), asgfr = IS_ASIG_ARG(p->kfo);
577
578 in = p->ain;
579 out = p->sr;
580 if (UNLIKELY(p->kbw[0] <= FL(0.0))) {
581 memset(out, 0, CS_KSMPS*sizeof(MYFLT));
582 return OK;
583 }
584 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
585 if (UNLIKELY(early)) {
586 nsmps -= early;
587 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
588 }
589 /* if (p->kbw[0] != p->lkb || p->kfo[0] != p->lkf) { */
590 /* p->lkf = p->kfo[0]; */
591 /* p->lkb = p->kbw[0]; */
592 /* c = 1.0 / tan((double)(csound->pidsr * p->lkb)); */
593 /* d = 2.0 * cos((double)(csound->tpidsr * p->lkf)); */
594 /* a[1] = 1.0 / (1.0 + c); */
595 /* a[2] = 0.0; */
596 /* a[3] = -a[1]; */
597 /* a[4] = - c * d * a[1]; */
598 /* a[5] = (c - 1.0) * a[1]; */
599 /* } */
600 //butter_filter(nsmps, offset, in, out, p->a);
601 for (nn=offset; nn<nsmps; nn++) {
602 MYFLT bw, fr;
603 bw = (asgbw ? p->kbw[nn] : *p->kbw);
604 fr = (asgfr ? p->kfo[nn] : *p->kfo);
605 if (bw != p->lkb || fr != p->lkf) {
606 double c, d;
607 p->lkf = fr;
608 p->lkb = bw;
609 c = 1.0 / tan((double)(csound->pidsr * bw));
610 d = 2.0 * cos((double)(csound->tpidsr * fr));
611 a[1] = 1.0 / (1.0 + c);
612 a[2] = 0.0;
613 a[3] = -a[1];
614 a[4] = - c * d * a[1];
615 a[5] = (c - 1.0) * a[1];
616 }
617 t = (double)in[nn] - a[4] * a[6] - a[5] * a[7];
618 t = csoundUndenormalizeDouble(t); /* Not needed on AMD */
619 y = t * a[1] + a[2] * a[6] + a[3] * a[7];
620 a[7] = a[6];
621 a[6] = t;
622 out[nn] = (MYFLT)y;
623 }
624 return OK;
625 }
626
bpcutxx(CSOUND * csound,BBFIL * p)627 static int32_t bpcutxx(CSOUND *csound, BBFIL *p) /* Band reject filter */
628 {
629 uint32_t offset = p->h.insdshead->ksmps_offset;
630 uint32_t early = p->h.insdshead->ksmps_no_end;
631 uint32_t nsmps = CS_KSMPS;
632 MYFLT *out, *in;
633 double *a = p->a;
634 double t, y;
635 uint32_t nn;
636 int32_t
637 asgbw = IS_ASIG_ARG(p->kbw), asgfr = IS_ASIG_ARG(p->kfo);
638
639 in = p->ain;
640 out = p->sr;
641
642 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
643 if (UNLIKELY(early)) {
644 nsmps -= early;
645 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
646 }
647 if (UNLIKELY(p->kbw[0] <= FL(0.0))) {
648 memcpy(&out[offset], &in[offset], (nsmps-offset)*sizeof(MYFLT));
649 return OK;
650 }
651
652 for (nn=offset; nn<nsmps; nn++) {
653 MYFLT bw, fr;
654 bw = (asgbw ? p->kbw[nn] : *p->kbw);
655 fr = (asgfr ? p->kfo[nn] : *p->kfo);
656 if (bw != p->lkb || fr != p->lkf) {
657 double c, d;
658 p->lkf = fr;
659 p->lkb = bw;
660 c = tan((double)(csound->pidsr * bw));
661 d = 2.0 * cos((double)(csound->tpidsr * fr));
662 a[1] = 1.0 / (1.0 + c);
663 a[2] = - d * a[1];
664 a[3] = a[1];
665 a[4] = a[2];
666 a[5] = (1.0 - c) * a[1];
667 }
668 t = (double)in[nn] - a[4] * a[6] - a[5] * a[7];
669 t = csoundUndenormalizeDouble(t); /* Not needed on AMD */
670 y = t * a[1] + a[2] * a[6] + a[3] * a[7];
671 a[7] = a[6];
672 a[6] = t;
673 out[nn] = (MYFLT)y;
674 }
675 return OK;
676 }
677
678
679
680 static OENTRY afilts_localops[] =
681 {
682 { "areson.aa", sizeof(RESON), 0,3,"a","aaaoo",(SUBR)rsnset,(SUBR)aresonaa},
683 { "areson.ak", sizeof(RESON), 0,3,"a","aakoo",(SUBR)rsnset,(SUBR)aresonak},
684 { "areson.ka", sizeof(RESON), 0,3,"a","akaoo",(SUBR)rsnset,(SUBR)aresonka},
685 { "atone.a", sizeof(TONE), 0,3,"a","ako", (SUBR)tonset,(SUBR)atonea },
686 { "atonex.a", sizeof(TONEX), 0,3, "a","aaoo",(SUBR)tonsetx,(SUBR)atonexa},
687 { "tone.a", sizeof(TONE), 0,3,"a","aao", (SUBR)tonset,(SUBR)tonea },
688 { "tonex.a", sizeof(TONEX), 0,3,"a","aaoo", (SUBR)tonsetx,(SUBR)tonexa },
689 { "butterhp.a", sizeof(BFIL), 0,3,"a","aao", (SUBR)butset,(SUBR)hibuta },
690 { "butterlp.a", sizeof(BFIL), 0,3,"a","aao", (SUBR)butset,(SUBR)lobuta },
691 { "buthp.a", sizeof(BFIL), 0,3,"a","aao", (SUBR)butset,(SUBR)hibuta },
692 { "butlp.a", sizeof(BFIL), 0,3,"a","aao", (SUBR)butset,(SUBR)lobuta },
693 { "butterbp", sizeof(BBFIL),0,3,"a","axxo", (SUBR)bbutset,(SUBR)bppasxx},
694 { "butbp", sizeof(BBFIL),0,3,"a","axxo", (SUBR)bbutset,(SUBR)bppasxx},
695 { "butterbr", sizeof(BBFIL),0,3,"a","axxo", (SUBR)bbutset,(SUBR)bpcutxx},
696 { "butbr", sizeof(BBFIL),0,3,"a","axxo", (SUBR)bbutset,(SUBR)bpcutxx},
697 };
698
699 LINKAGE_BUILTIN(afilts_localops)
700
701
702