« back to all changes in this revision
Viewing changes to filter.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Matt Zimmerman
- Date: 2001-03-14 17:18:53 UTC
- Revision ID: james.westby@ubuntu.com-20010314171853-actxdsm1k2otn4ok
Tags: upstream-0.3.2
Import upstream version 0.3.2
- files added:
- raw
- xpm
added
removed
filter.c
1
/* Freebirth
2
* Copyright (C) 1999 topher lafata <topher@topher.com>,
3
* Jake Donham <jake@bitmechanic.com>
4
*
5
* This program is free software; you can redistribute it and/or
6
* modify it under the terms of the GNU General Public
7
* License as published by the Free Software Foundation; either
8
* version 2 of the License, or (at your option) any later version.
9
*
10
* This program is distributed in the hope that it will be useful,
11
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13
* General Public License for more details.
14
*
15
* You should have received a copy of the GNU General Public
16
* License along with this program (see COPYING); if not, write to the
17
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
18
* Boston, MA 02111-1307, USA.
19
*/
20
21
#include "filter.h"
22
#include <stdlib.h>
23
#include <stdio.h>
24
#include <math.h>
25
26
#define FILTER_SECTIONS 2 /* 2 filter sections for 24 db/oct filter */
27
#define MIN_GAIN .4
28
#define MAX_GAIN .7
29
30
#define SAMPLE_RATE 44100
31
32
33
void filter_next_buffer(filter *this);
34
35
typedef struct {
36
double a0, a1, a2; /* numerator coefficients */
37
double b0, b1, b2; /* denominator coefficients */
38
} BIQUAD;
39
40
BIQUAD ProtoCoef[FILTER_SECTIONS] = {
41
{ 1.0, 0, 0, 1.0, 0.765367, 1.0 },
42
{ 1.0, 0, 0, 1.0, 1.847759, 1.0 }
43
};
44
45
46
static float **coef_lookup_table;
47
static int looked_up = 0;
48
49
50
51
/*coeficient calculation function */
52
void szxform(
53
double *a0, double *a1, double *a2, /* numerator coefficients */
54
double *b0, double *b1, double *b2, /* denominator coefficients */
55
double fc, /* Filter cutoff frequency */
56
double fs, /* sampling rate */
57
double *k, /* overall gain factor */
58
float *coef); /* pointer to 4 iir coefficients */
59
60
float *filter_get_coef(int sections,double cutoff,double q, double gain)
61
{
62
63
int nInd;
64
65
double k = gain;
66
float *Coef = (float *)calloc(4 * sections + 1, sizeof(float));
67
float *coef = Coef + 1;
68
double Q = q;
69
double fc = cutoff;
70
double fs = SAMPLE_RATE;
71
72
double a0, a1, a2, b0, b1, b2;
73
74
75
/* printf("cutoff is %f\n",fc); */
76
/* printf("Q is %f\n",Q); */
77
78
for(nInd = 0; nInd < sections; nInd++)
79
{
80
a0 = ProtoCoef[nInd].a0;
81
a1 = ProtoCoef[nInd].a1;
82
a2 = ProtoCoef[nInd].a2;
83
84
b0 = ProtoCoef[nInd].b0;
85
b1 = ProtoCoef[nInd].b1 / Q; /* Divide by resonance or Q */
86
b2 = ProtoCoef[nInd].b2;
87
szxform(&a0, &a1, &a2, &b0, &b1, &b2, fc, fs, &k, coef);
88
coef += 4; /* Point to next filter section */
89
90
}
91
92
Coef[0] = k;
93
return Coef;
94
}
95
96
97
void filter_init_lookup_table()
98
{
99
int index;
100
double q,fc;
101
int freq_range, q_range, table_length;
102
float q_scale ;
103
float gain_scale;
104
105
q_scale = MAX_Q/MAX_FC;
106
gain_scale = - (MAX_GAIN - MIN_GAIN) / MAX_FC;
107
freq_range = MAX_FC - MIN_FC;
108
q_range = (MAX_Q - MIN_Q) / (double)Q_INCREMENT;
109
110
table_length = freq_range * q_range;
111
112
printf("Calculating coefficient table.");
113
114
index = 0;
115
for(fc = MIN_FC; fc < MAX_FC;fc++)
116
for(q = MIN_Q;q < MAX_Q; q += (double)Q_INCREMENT){
117
/*PATCH BY: Harold Gutch <logix@foobar.franken.de>
118
FreeBSD can't handle devide by zero and it was possible
119
for scaled_q to go to zero*/
120
double scaled_q = (((q_scale * (fc)) / MAX_Q)* q) + MIN_Q;
121
coef_lookup_table[index] = filter_get_coef(FILTER_SECTIONS,
122
fc,
123
scaled_q,
124
gain_scale * fc + MAX_GAIN);
125
126
index++;
127
if(index % 1000 == 0){
128
printf(".");
129
fflush(0);
130
}
131
}
132
looked_up = 1;
133
printf("\nDone\n");
134
}
135
136
void filter_set_coeffs(filter *this, int fc, float q)
137
{
138
int index;
139
index = ((fc - MIN_FC) * (MAX_Q - MIN_Q) / Q_INCREMENT)
140
+ (q / Q_INCREMENT);
141
142
this->coef = coef_lookup_table[index];
143
144
}
145
146
147
148
149
sample filter_process_sample(filter *this, sample input)
150
{
151
152
unsigned int i;
153
float *hist1_ptr,*hist2_ptr,*coef_ptr;
154
float output,new_hist,history1,history2;
155
156
157
/* allocate history array if different size than last call */
158
if(!this->history)
159
this->history = (float *) calloc(2*this->length,sizeof(float));
160
161
coef_ptr = this->coef; /* coefficient pointer */
162
163
hist1_ptr = this->history; /* first history */
164
hist2_ptr = hist1_ptr + 1; /* next history */
165
166
167
output = (float)input * (*coef_ptr++);
168
169
for (i = 0 ; i < this->length; i++)
170
{
171
history1 = *hist1_ptr; /* history values */
172
history2 = *hist2_ptr;
173
174
output = output - history1 * (*coef_ptr++);
175
new_hist = output - history2 * (*coef_ptr++); /* poles */
176
177
output = new_hist + history1 * (*coef_ptr++);
178
output = output + history2 * (*coef_ptr++); /* zeros */
179
180
*hist2_ptr++ = *hist1_ptr;
181
*hist1_ptr++ = new_hist;
182
hist1_ptr++;
183
hist2_ptr++;
184
}
185
186
return (sample)output;
187
}
188
189
void filter_next_buffer(filter *this) {
190
this->next = 1;
191
}
192
193
sample *filter_get_buffer(filter *this)
194
{
195
int i;
196
sample *filt_e_buffer;
197
sample *buffer;
198
199
if (this->next) {
200
this->next = 0;
201
buffer = this->sp->get_buffer(this->sp);
202
filt_e_buffer = this->filt_e->get_buffer(this->filt_e);
203
for(i = 0; i < TBASS_BUFF_SIZE; i++)
204
{
205
float s;
206
207
/*
208
we're scaling the 0-SAMPLE_MAX range of the envelope into
209
a LIMIT-SAMPLE_MAX range 'cause the filter doesn't like the
210
cutoff to be 0
211
*/
212
int fe = filt_e_buffer[i] * (SAMPLE_MAX - LIMIT) / SAMPLE_MAX + LIMIT;
213
int fc = this->fc * fe / SAMPLE_MAX;
214
filter_set_coeffs(this,
215
fc,
216
this->Q);
217
s = filter_process_sample(this,buffer[i]);
218
219
/*
220
PATCH: Benno Senoner <benno@gardena.net>
221
clip on underflow as well as overflow
222
*/
223
if(s > SAMPLE_MAX)
224
s = SAMPLE_MAX;
225
else if (s < SAMPLE_MIN)
226
s = SAMPLE_MIN;
227
228
this->buffer[i] = s;
229
230
}
231
}
232
return this->buffer;
233
}
234
235
236
void filter_set_cutoff(filter *this,double cutoff)
237
{
238
this->fc = cutoff;
239
filter_set_coeffs(this,this->fc,this->Q);
240
241
}
242
243
void filter_set_gain(filter *this,double gain)
244
{
245
this->gain = gain;
246
247
}
248
249
void filter_set_q(filter *this,double q)
250
{
251
this->Q = q;
252
filter_set_coeffs(this,this->fc,this->Q);
253
}
254
255
256
void filter_delete(filter *this)
257
{
258
free(this->history);
259
free(this->coef);
260
free(this);
261
262
}
263
264
265
266
/*
267
* ----------------------------------------------------------
268
* bilinear.c
269
*
270
* Perform bilinear transformation on s-domain coefficients
271
* of 2nd order biquad section.
272
* First design an analog filter and use s-domain coefficients
273
* as input to szxform() to convert them to z-domain.
274
*
275
* Here's the butterworth polinomials for 2nd, 4th and 6th order sections.
276
* When we construct a 24 db/oct filter, we take to 2nd order
277
* sections and compute the coefficients separately for each section.
278
*
279
* n Polinomials
280
* --------------------------------------------------------------------
281
* 2 s^2 + 1.4142s +1
282
* 4 (s^2 + 0.765367s + 1) (s^2 + 1.847759s + 1)
283
* 6 (s^2 + 0.5176387s + 1) (s^2 + 1.414214 + 1) (s^2 + 1.931852s + 1)
284
*
285
* Where n is a filter order.
286
* For n=4, or two second order sections, we have following equasions for each
287
* 2nd order stage:
288
*
289
* (1 / (s^2 + (1/Q) * 0.765367s + 1)) * (1 / (s^2 + (1/Q) * 1.847759s + 1))
290
*
291
* Where Q is filter quality factor in the range of
292
* 1 to 1000. The overall filter Q is a product of all
293
* 2nd order stages. For example, the 6th order filter
294
* (3 stages, or biquads) with individual Q of 2 will
295
* have filter Q = 2 * 2 * 2 = 8.
296
*
297
* The nominator part is just 1.
298
* The denominator coefficients for stage 1 of filter are:
299
* b2 = 1; b1 = 0.765367; b0 = 1;
300
* numerator is
301
* a2 = 0; a1 = 0; a0 = 1;
302
*
303
* The denominator coefficients for stage 1 of filter are:
304
* b2 = 1; b1 = 1.847759; b0 = 1;
305
* numerator is
306
* a2 = 0; a1 = 0; a0 = 1;
307
*
308
* These coefficients are used directly by the szxform()
309
* and bilinear() functions. For all stages the numerator
310
* is the same and the only thing that is different between
311
* different stages is 1st order coefficient. The rest of
312
* coefficients are the same for any stage and equal to 1.
313
*
314
* Any filter could be constructed using this approach.
315
*
316
* References:
317
* Van Valkenburg, "Analog Filter Design"
318
* Oxford University Press 1982
319
* ISBN 0-19-510734-9
320
*
321
* C Language Algorithms for Digital Signal Processing
322
* Paul Embree, Bruce Kimble
323
* Prentice Hall, 1991
324
* ISBN 0-13-133406-9
325
*
326
* Digital Filter Designer's Handbook
327
* With C++ Algorithms
328
* Britton Rorabaugh
329
* McGraw Hill, 1997
330
* ISBN 0-07-053806-9
331
* ----------------------------------------------------------
332
*/
333
334
335
void prewarp(double *a0, double *a1, double *a2, double fc, double fs);
336
void bilinear(
337
double a0, double a1, double a2, /* numerator coefficients */
338
double b0, double b1, double b2, /* denominator coefficients */
339
double *k, /* overall gain factor */
340
double fs, /* sampling rate */
341
float *coef); /* pointer to 4 iir coefficients */
342
343
344
/*
345
* ----------------------------------------------------------
346
* Pre-warp the coefficients of a numerator or denominator.
347
* Note that a0 is assumed to be 1, so there is no wrapping
348
* of it.
349
* ----------------------------------------------------------
350
*/
351
void prewarp(
352
double *a0, double *a1, double *a2,
353
double fc, double fs)
354
{
355
double wp, pi;
356
357
pi = 4.0 * atan(1.0);
358
wp = 2.0 * fs * tan(pi * fc / fs);
359
360
if(wp == 0) /*PATCH BY: free bsd from Hans Huebner<hans@huebner.org> */
361
wp = 0.0001;
362
*a2 = (*a2) / (wp * wp);
363
*a1 = (*a1) / wp;
364
}
365
366
367
368
369
370
371
/*
372
* ----------------------------------------------------------
373
* bilinear()
374
*
375
* Transform the numerator and denominator coefficients
376
* of s-domain biquad section into corresponding
377
* z-domain coefficients.
378
*
379
* Store the 4 IIR coefficients in array pointed by coef
380
* in following order:
381
* beta1, beta2 (denominator)
382
* alpha1, alpha2 (numerator)
383
*
384
* Arguments:
385
* a0-a2 - s-domain numerator coefficients
386
* b0-b2 - s-domain denominator coefficients
387
* k - filter gain factor. initially set to 1
388
* and modified by each biquad section in such
389
* a way, as to make it the coefficient by
390
* which to multiply the overall filter gain
391
* in order to achieve a desired overall filter gain,
392
* specified in initial value of k.
393
* fs - sampling rate (Hz)
394
* coef - array of z-domain coefficients to be filled in.
395
*
396
* Return:
397
* On return, set coef z-domain coefficients
398
* ----------------------------------------------------------
399
*/
400
void bilinear(
401
double a0, double a1, double a2, /* numerator coefficients */
402
double b0, double b1, double b2, /* denominator coefficients */
403
double *k, /* overall gain factor */
404
double fs, /* sampling rate */
405
float *coef /* pointer to 4 iir coefficients */
406
)
407
{
408
double ad, bd;
409
410
/* alpha (Numerator in s-domain) */
411
ad = 4. * a2 * fs * fs + 2. * a1 * fs + a0;
412
/* beta (Denominator in s-domain) */
413
bd = 4. * b2 * fs * fs + 2. * b1* fs + b0;
414
415
/* update gain constant for this section */
416
*k *= ad/bd;
417
418
/* Denominator */
419
*coef++ = (2. * b0 - 8. * b2 * fs * fs)
420
/ bd; /* beta1 */
421
*coef++ = (4. * b2 * fs * fs - 2. * b1 * fs + b0)
422
/ bd; /* beta2 */
423
424
/* Nominator */
425
*coef++ = (2. * a0 - 8. * a2 * fs * fs)
426
/ ad; /* alpha1 */
427
*coef = (4. * a2 * fs * fs - 2. * a1 * fs + a0)
428
/ ad; /* alpha2 */
429
}
430
431
432
/*
433
* ----------------------------------------------------------
434
* Transform from s to z domain using bilinear transform
435
* with prewarp.
436
*
437
* Arguments:
438
* For argument description look at bilinear()
439
*
440
* coef - pointer to array of floating point coefficients,
441
* corresponding to output of bilinear transofrm
442
* (z domain).
443
*
444
* Note: frequencies are in Hz.
445
* ----------------------------------------------------------
446
*/
447
void szxform(
448
double *a0, double *a1, double *a2, /* numerator coefficients */
449
double *b0, double *b1, double *b2, /* denominator coefficients */
450
double fc, /* Filter cutoff frequency */
451
double fs, /* sampling rate */
452
double *k, /* overall gain factor */
453
float *coef) /* pointer to 4 iir coefficients */
454
{
455
/* Calculate a1 and a2 and overwrite the original values */
456
prewarp(a0, a1, a2, fc, fs);
457
prewarp(b0, b1, b2, fc, fs);
458
bilinear(*a0, *a1, *a2, *b0, *b1, *b2, k, fs, coef);
459
}
460
461
462
463
static sample_producer **get_children(filter *this)
464
{
465
return &(this->sp);
466
}
467
468
static char **get_header(filter *this)
469
{
470
static char *header[] = {
471
"int $n_fc = ((filter *)$t)->fc;",
472
"float $n_Q = ((filter *)$t)->Q;",
473
"float $n_s;",
474
NULL
475
};
476
return header;
477
}
478
479
static char **get_code(filter *this)
480
{
481
static char *code[] = {
482
"filter_set_coeffs((filter *)$t, $n_fc *",
483
" ($i1 * (SAMPLE_MAX - LIMIT) / SAMPLE_MAX + LIMIT) /",
484
" SAMPLE_MAX, $n_Q);",
485
"$n_s = filter_process_sample((filter *)$t, $i0);",
486
"if ($n_s > SAMPLE_MAX) $n_s = SAMPLE_MAX;",
487
"else if ($n_s < SAMPLE_MIN) $n_s = SAMPLE_MIN;",
488
"$o = $n_s;",
489
NULL
490
};
491
return code;
492
}
493
494
static char **get_footer(filter *this)
495
{
496
static char *footer[] = { NULL };
497
return footer;
498
}
499
500
filter *filter_new(double cutoff, double q, sample_producer *filt_e,
501
sample_producer* sp)
502
{
503
int i;
504
filter *out = malloc(sizeof(filter));
505
out->get_buffer = filter_get_buffer;
506
out->next_buffer = filter_next_buffer;
507
out->get_children = get_children;
508
out->get_header = get_header;
509
out->get_code = get_code;
510
out->get_footer = get_footer;
511
out->next = 0;
512
out->schedule = sp_schedule_event;
513
out->trigger = NULL; /*what does it mean to trigger the filter? */
514
515
/*we aren't responding to events yet in the filter */
516
out->events = (event **)malloc(sizeof(event *)
517
* TBASS_BUFF_SIZE);
518
519
for(i = 0; i < TBASS_BUFF_SIZE; i++)
520
out->events[i] = NULL;
521
522
523
out->length = FILTER_SECTIONS;
524
out->gain = MAX_GAIN;
525
out->coef = (float *) calloc(4 * out->length + 1, sizeof(float));
526
out->fc = cutoff;
527
out->Q = q;
528
out->sp = sp;
529
out->filt_e = filt_e;
530
out->unused = NULL;
531
out->buffer = (sample *)malloc(sizeof(sample) * TBASS_BUFF_SIZE);
532
533
/* we're doing this here because we don't call filter_init_lookup_table
534
in fusebirth, and we need to create the filter but not have it
535
blow up in filter_set_coeffs because the table doesn't exist. */
536
coef_lookup_table = (float **)
537
malloc(sizeof(float*) *
538
(MAX_FC - MIN_FC) * (MAX_Q - MIN_Q) / (double)Q_INCREMENT);
539
540
filter_set_coeffs(out,out->fc,out->Q);
541
return out;
542
543
}
544
545
/*
546
Local Variables:
547
mode: font-lock
548
End:
549
*/
Loggerhead is a web-based interface for Breezy
Version: 2.0.1