1 /*
2     wdatutil.c - Code to generate wavedata for bandlimited waveforms.
3 
4     Copyright (C) 2003  Mike Rawes
5 
6     This program is free software; you can redistribute it and/or modify
7     it under the terms of the GNU General Public License as published by
8     the Free Software Foundation; either version 2 of the License, or
9     (at your option) any later version.
10 
11     This program is distributed in the hope that it will be useful,
12     but WITHOUT ANY WARRANTY; without even the implied warranty of
13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14     GNU General Public License for more details.
15 
16     You should have received a copy of the GNU General Public License
17     along with this program; if not, write to the Free Software
18     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19 */
20 
21 #include <string.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <math.h>
25 #include <ladspa.h>
26 #include "common.h"
27 #include "wavedata.h"
28 #include "wdatutil.h"
29 
30 #ifdef __cplusplus
31 extern "C" {
32 #endif
33 
34 void
35 generate_sine (LADSPA_Data * samples,
36                unsigned long sample_count);
37 void
38 generate_sawtooth (LADSPA_Data * samples,
39                    unsigned long sample_count,
40                    unsigned long harmonics,
41                    float gibbs_comp);
42 void
43 generate_square (LADSPA_Data * samples,
44                  unsigned long sample_count,
45                  unsigned long harmonics,
46                  float gibbs_comp);
47 void
48 generate_parabola (LADSPA_Data * samples,
49                    unsigned long sample_count,
50                    unsigned long harmonics,
51                    float gibbs_comp);
52 
53 #ifdef __cplusplus
54 } /* extern "C" { */
55 #endif
56 
57 char *wave_names[] = {
58 	"Saw",
59 	"Square",
60 	"Parabola"
61 };
62 
63 char *wave_descriptions[] = {
64 	"Sawtooth Wave",
65 	"Square Wave",
66 	"Parabola Wave"
67 };
68 
69 unsigned long wave_first_harmonics[] = {
70 	1,
71 	1,
72 	1
73 };
74 
75 unsigned long wave_harmonic_intervals[] = {
76 	1,
77 	2,
78 	1
79 };
80 
81 Wavedata *
wavedata_new(unsigned long sample_rate)82 wavedata_new (unsigned long sample_rate)
83 {
84 	Wavedata * w;
85 
86 	w = (Wavedata *) malloc (sizeof (Wavedata));
87 
88 	if (!w)
89 		return 0;
90 
91 	w->data_handle = 0;
92 	w->table_count = 0;
93 	w->tables = 0;
94 	w->lookup = 0;
95 	w->lookup_max = 0;
96 	w->sample_rate = (LADSPA_Data) sample_rate;
97 	w->nyquist = w->sample_rate * 0.5f;
98 
99 	return w;
100 }
101 
102 void
wavedata_cleanup(Wavedata * w)103 wavedata_cleanup (Wavedata * w)
104 {
105 	unsigned long ti;
106 	Wavetable * t;
107 
108 	for (ti = 0; ti < w->table_count; ti++)
109 	{
110 		t = w->tables[ti];
111 		if (t)
112 		{
113 			if (t->samples_hf)
114 				free (t->samples_hf);
115 
116 			if (t->samples_lf)
117 				free (t->samples_lf);
118 
119 			free (t);
120 		}
121 	}
122 
123 	free (w);
124 }
125 
126 int
wavedata_add_table(Wavedata * w,unsigned long sample_count,unsigned long harmonics)127 wavedata_add_table (Wavedata * w,
128                     unsigned long sample_count,
129                     unsigned long harmonics)
130 {
131 	Wavetable ** tables;
132 	Wavetable * t;
133 	size_t bytes;
134 
135 	t = (Wavetable *) malloc (sizeof (Wavetable));
136 
137 	if (!t)
138 		return -1;
139 
140 /* Extra 3 samples for interpolation */
141 	bytes = (sample_count + 3) * sizeof (LADSPA_Data);
142 
143 	t->samples_lf = (LADSPA_Data *) malloc (bytes);
144 
145 	if (!t->samples_lf)
146 	{
147 		free (t);
148 		return -1;
149 	}
150 
151 	t->samples_hf = (LADSPA_Data *) malloc (bytes);
152 
153 	if (!t->samples_hf)
154 	{
155 		free (t->samples_lf);
156 		free (t);
157 		return -1;
158 	}
159 
160 	bytes = (w->table_count + 1) * sizeof (Wavetable *);
161 	if (w->table_count == 0)
162 		tables = (Wavetable **) malloc (bytes);
163 	else
164 		tables = (Wavetable **) realloc (w->tables, bytes);
165 
166 	if (!tables)
167 	{
168 		free (t);
169 		return -1;
170 	}
171 
172 	t->sample_count = sample_count;
173 	t->harmonics = harmonics;
174 
175 	if (w->lookup_max < harmonics)
176 		w->lookup_max = harmonics;
177 
178 	tables[w->table_count] = t;
179 	w->tables = tables;
180 	w->table_count++;
181 
182 	return 0;
183 }
184 
185 void
wavedata_generate_tables(Wavedata * w,Wavetype wavetype,float gibbs_comp)186 wavedata_generate_tables (Wavedata * w,
187                           Wavetype wavetype,
188                           float gibbs_comp)
189 {
190 	Wavetable * t;
191 	LADSPA_Data * samples_lf;
192 	LADSPA_Data * samples_hf;
193 	unsigned long h_lf;
194 	unsigned long h_hf;
195 	unsigned long s;
196 	unsigned long i;
197 
198 	for (i = 0; i < w->table_count; i++)
199 	{
200 		t = w->tables[i];
201 
202 		h_lf = t->harmonics;
203 
204 		if (i < w->table_count - 1)
205 			h_hf = w->tables[i+1]->harmonics;
206 		else
207 			h_hf = 1;
208 
209 		samples_lf = t->samples_lf;
210 		samples_hf = t->samples_hf;
211 		samples_lf++;
212 		samples_hf++;
213 
214 		switch (wavetype)
215 		{
216 			case SAW:
217 				generate_sawtooth (samples_lf, t->sample_count, h_lf, gibbs_comp);
218 				generate_sawtooth (samples_hf, t->sample_count, h_hf, gibbs_comp);
219 				break;
220 			case SQUARE:
221 				generate_square (samples_lf, t->sample_count, h_lf, gibbs_comp);
222 				generate_square (samples_hf, t->sample_count, h_hf, gibbs_comp);
223 				break;
224 			case PARABOLA:
225 				generate_parabola (samples_lf, t->sample_count, h_lf, gibbs_comp);
226 				generate_parabola (samples_hf, t->sample_count, h_hf, gibbs_comp);
227 				break;
228 		}
229 
230 	/* Basic denormalization */
231 		for (s = 0; s < t->sample_count; s++)
232 			samples_lf[s] = FABSF (samples_lf[s]) < SMALLEST_FLOAT ? 0.0 : samples_lf[s];
233 
234 		samples_lf--;
235 		samples_lf[0] = samples_lf[t->sample_count];
236 		samples_lf[t->sample_count + 1] = samples_hf[1];
237 		samples_lf[t->sample_count + 2] = samples_hf[2];
238 
239 		for (s = 0; s < t->sample_count; s++)
240 			samples_hf[s] = FABSF (samples_hf[s]) < SMALLEST_FLOAT ? 0.0 : samples_hf[s];
241 
242 		samples_hf--;
243 		samples_hf[0] = samples_hf[t->sample_count];
244 		samples_hf[t->sample_count + 1] = samples_hf[1];
245 		samples_hf[t->sample_count + 2] = samples_hf[2];
246 	}
247 }
248 
249 int
wavedata_write(Wavedata * w,FILE * wdat_fp,char * data_name)250 wavedata_write (Wavedata * w,
251                 FILE * wdat_fp,
252                 char * data_name)
253 {
254 	Wavetable * t = 0;
255 	unsigned long table_count;
256 	unsigned long i;
257 	unsigned long j;
258 	unsigned long s;
259 	int column;
260 /*
261  * Extra table at end
262  */
263 	table_count = w->table_count + 1;
264 
265 	fprintf (wdat_fp, "#include <ladspa.h>\n");
266 	fprintf (wdat_fp, "#include <stdio.h>\n");
267 	fprintf (wdat_fp, "#include \"wavedata.h\"\n");
268 	fprintf (wdat_fp, "\n");
269 /*
270  * Fixed data and tables
271  */
272  	fprintf (wdat_fp, "unsigned long ref_count = 0;\n");
273 	fprintf (wdat_fp, "unsigned long first_sample_rate = 0;\n");
274 	fprintf (wdat_fp, "unsigned long table_count = %ld;\n", table_count);
275 	fprintf (wdat_fp, "Wavetable tables[%ld];\n", table_count);
276 	fprintf (wdat_fp, "Wavetable * ptables[%ld];\n", table_count);
277 	fprintf (wdat_fp, "unsigned long lookup[%ld];\n", w->lookup_max + 1);
278 	fprintf (wdat_fp, "unsigned long lookup_max = %ld;\n", w->lookup_max);
279 	fprintf (wdat_fp, "\n");
280 /*
281  * Sample data
282  * Each table has an extra 3 samples for interpolation
283  */
284 	for (i = 0; i < w->table_count; i++)
285 	{
286 		t = w->tables[i];
287 
288 		fprintf(wdat_fp, "static LADSPA_Data samples_lf_%ld[%ld] = {\n", i, t->sample_count + 3);
289 
290 		column = 0;
291 		for (s = 0; s < t->sample_count + 3 - 1; s++, column++)
292 		{
293 			if (column == 5)
294 			{
295 				fprintf (wdat_fp, "\n");
296 				column = 0;
297 			}
298 			fprintf (wdat_fp, "%+.8ef,", t->samples_lf[s]);
299 		}
300 
301 		if (column == 5)
302 			fprintf (wdat_fp, "\n");
303 
304 		fprintf (wdat_fp, "%+.8ef\n", t->samples_lf[s]);
305 		fprintf (wdat_fp, "};\n");
306 		fprintf (wdat_fp, "\n");
307 
308 		fprintf(wdat_fp, "static LADSPA_Data samples_hf_%ld[%ld] = {\n", i, t->sample_count + 3);
309 
310 		column = 0;
311 		for (s = 0; s < t->sample_count + 3 - 1; s++, column++)
312 		{
313 			if (column == 5)
314 			{
315 				fprintf (wdat_fp, "\n");
316 				column = 0;
317 			}
318 			fprintf (wdat_fp, "%+.8ef,", t->samples_hf[s]);
319 		}
320 
321 		if (column == 5)
322 			fprintf (wdat_fp, "\n");
323 
324 		fprintf (wdat_fp, "%+.8ef\n", t->samples_hf[s]);
325 		fprintf (wdat_fp, "};\n");
326 		fprintf (wdat_fp, "\n");
327 	}
328 
329 	fprintf (wdat_fp, "LADSPA_Data samples_zero[%ld];\n", t->sample_count + 3);
330 	fprintf (wdat_fp, "\n");
331 /*
332  * Function to get Wavedata - the sample rate is needed to calculate
333  * frequencies and related things
334  */
335 	fprintf (wdat_fp, "int\n");
336 	fprintf (wdat_fp, "blop_get_%s (Wavedata * w, unsigned long sample_rate)\n", data_name);
337 	fprintf (wdat_fp, "{\n");
338 	fprintf (wdat_fp, "\tWavetable * t;\n");
339 	fprintf (wdat_fp, "\tunsigned long ti;\n");
340 	fprintf (wdat_fp, "\n");
341 /*
342  * Sample rate must be > 0
343  */
344 	fprintf (wdat_fp, "\tif (sample_rate == 0)\n");
345 	fprintf (wdat_fp, "\t\treturn -1;\n");
346 	fprintf (wdat_fp, "\n");
347 /*
348  * First time call - set up all sample rate dependent data
349  */
350 	fprintf (wdat_fp, "\tif (first_sample_rate == 0)\n");
351 	fprintf (wdat_fp, "\t{\n");
352 	fprintf (wdat_fp, "\t\tfirst_sample_rate = sample_rate;\n");
353 	fprintf (wdat_fp, "\t\tw->sample_rate = (LADSPA_Data) sample_rate;\n");
354 	fprintf (wdat_fp, "\t\tw->nyquist = w->sample_rate / 2.0f;\n");
355 	fprintf (wdat_fp, "\t\tw->table_count = table_count;\n");
356 	fprintf (wdat_fp, "\t\tw->tables = ptables;\n");
357 	fprintf (wdat_fp, "\t\tw->lookup = lookup;\n");
358 	fprintf (wdat_fp, "\t\tw->lookup_max = lookup_max;\n");
359 	fprintf (wdat_fp, "\n");
360 	fprintf (wdat_fp, "\t\tfor (ti = 1; ti < table_count - 1; ti++)\n");
361 	fprintf (wdat_fp, "\t\t{\n");
362 	fprintf (wdat_fp, "\t\t\tt = ptables[ti];\n");
363 	fprintf (wdat_fp, "\t\t\tt->min_frequency = w->nyquist / (LADSPA_Data) (ptables[ti - 1]->harmonics);\n");
364 	fprintf (wdat_fp, "\t\t\tt->max_frequency = w->nyquist / (LADSPA_Data) (t->harmonics);\n");
365 	fprintf (wdat_fp, "\t\t}\n");
366 	fprintf (wdat_fp, "\n");
367 	fprintf (wdat_fp, "\t\tt = w->tables[0];\n");
368 	fprintf (wdat_fp, "\t\tt->min_frequency = 0.0f;\n");
369 	fprintf (wdat_fp, "\t\tt->max_frequency = ptables[1]->min_frequency;\n");
370 	fprintf (wdat_fp, "\n");
371 	fprintf (wdat_fp, "\t\tt = ptables[table_count - 1];\n");
372 	fprintf (wdat_fp, "\t\tt->min_frequency = ptables[w->table_count - 2]->max_frequency;\n");
373 	fprintf (wdat_fp, "\t\tt->max_frequency = w->nyquist;\n");
374 	fprintf (wdat_fp, "\t\n");
375 	fprintf (wdat_fp, "\t\tfor (ti = 0; ti < w->table_count; ti++)\n");
376 	fprintf (wdat_fp, "\t\t{\n");
377 	fprintf (wdat_fp, "\t\t\tt = w->tables[ti];\n");
378 	fprintf (wdat_fp, "\t\t\tt->phase_scale_factor = (LADSPA_Data) (t->sample_count) / w->sample_rate;\n");
379 	fprintf (wdat_fp, "\t\t\tt->range_scale_factor = 1.0f / (t->max_frequency - t->min_frequency);\n");
380 	fprintf (wdat_fp, "\t\t}\n");
381 	fprintf (wdat_fp, "\n");
382 	fprintf (wdat_fp, "\t\treturn 0;\n");
383 	fprintf (wdat_fp, "\t}\n");
384 /*
385  * Already called at least once, so just set up wavedata
386  */
387 	fprintf (wdat_fp, "\telse if (sample_rate == first_sample_rate)\n");
388 	fprintf (wdat_fp, "\t{\n");
389 	fprintf (wdat_fp, "\t\tw->sample_rate = (LADSPA_Data) sample_rate;\n");
390 	fprintf (wdat_fp, "\t\tw->nyquist = w->sample_rate / 2.0f;\n");
391 	fprintf (wdat_fp, "\t\tw->table_count = table_count;\n");
392 	fprintf (wdat_fp, "\t\tw->tables = ptables;\n");
393 	fprintf (wdat_fp, "\t\tw->lookup = lookup;\n");
394 	fprintf (wdat_fp, "\t\tw->lookup_max = lookup_max;\n");
395 	fprintf (wdat_fp, "\n");
396 	fprintf (wdat_fp, "\t\treturn 0;\n");
397 	fprintf (wdat_fp, "\t}\n");
398 /*
399  * Sample rate does not match, so fail
400  *
401  * NOTE: This means multiple sample rates are not supported
402  *       This should not present any problems
403  */
404 	fprintf (wdat_fp, "\telse\n");
405 	fprintf (wdat_fp, "\t{\n");
406 	fprintf (wdat_fp, "\t\treturn -1;\n");
407 	fprintf (wdat_fp, "\t}\n");
408 	fprintf (wdat_fp, "}\n");
409 	fprintf (wdat_fp, "\n");
410 /*
411  * _init()
412  * Assemble tables and lookup
413  */
414 	fprintf (wdat_fp, "void\n");
415 	fprintf (wdat_fp, "_init (void)\n");
416 	fprintf (wdat_fp, "{\n");
417 	fprintf (wdat_fp, "\tunsigned long max_harmonic;\n");
418 	fprintf (wdat_fp, "\tunsigned long ti;\n");
419 	fprintf (wdat_fp, "\tunsigned long li;\n");
420 	fprintf (wdat_fp, "\tunsigned long s;\n");
421 	fprintf (wdat_fp, "\n");
422 
423 	for (i = 0; i < w->table_count; i++)
424 	{
425 		t = w->tables[i];
426 
427 		fprintf (wdat_fp, "\ttables[%ld].sample_count = %ld;\n", i, t->sample_count);
428 		fprintf (wdat_fp, "\ttables[%ld].samples_lf = samples_lf_%ld;\n", i, i);
429 		fprintf (wdat_fp, "\ttables[%ld].samples_hf = samples_hf_%ld;\n", i, i);
430 		fprintf (wdat_fp, "\ttables[%ld].harmonics = %ld;\n", i, t->harmonics);
431 		fprintf (wdat_fp, "\n");
432 	}
433 /*
434  * Last table - uses same sample data as previous table for lf data,
435  * and zeroes for hf data
436  */
437 	i = w->table_count - 1;
438 	j = i + 1;
439 	t = w->tables[i];
440 /*
441  * Zero silent samples
442  */
443 	fprintf (wdat_fp, "\tfor (s = 0; s < %ld; s++)\n", t->sample_count + 3);
444 	fprintf (wdat_fp, "\t\tsamples_zero[s] = 0.0f;\n");
445 	fprintf (wdat_fp, "\n");
446 
447 	fprintf (wdat_fp, "\ttables[%ld].sample_count = %ld;\n", j, t->sample_count);
448 	fprintf (wdat_fp, "\ttables[%ld].samples_lf = samples_hf_%ld;\n", j, i);
449 	fprintf (wdat_fp, "\ttables[%ld].samples_hf = samples_zero;\n", j);
450 	fprintf (wdat_fp, "\ttables[%ld].harmonics = 1;\n", j);
451 	fprintf (wdat_fp, "\n");
452 /*
453  * Get pointers to each wavetable and put them in the pointer array
454  */
455 	fprintf (wdat_fp, "\tfor (ti = 0; ti < table_count; ti++)\n");
456 	fprintf (wdat_fp, "\t\tptables[ti] = &tables[ti];\n");
457 	fprintf (wdat_fp, "\n");
458 /*
459  * Shift all sample offsets forward by one sample
460  * !!! NO! Don't!
461 	fprintf (wdat_fp, "\tfor (ti = 0; ti < table_count; ti++)\n");
462 	fprintf (wdat_fp, "\t{\n");
463 	fprintf (wdat_fp, "\t\tptables[ti]->samples_lf++;\n");
464 	fprintf (wdat_fp, "\t\tptables[ti]->samples_hf++;\n");
465 	fprintf (wdat_fp, "\t}\n");
466 	fprintf (wdat_fp, "\n");
467  */
468 /*
469  * Table lookup vector indexed by harmonic
470  * Add lookup data to vector
471  */
472 	fprintf (wdat_fp, "\tli = 0;");
473 	fprintf (wdat_fp, "\n");
474 	fprintf (wdat_fp, "\tfor (ti = table_count - 1; ti > 0; ti--)\n");
475 	fprintf (wdat_fp, "\t{\n");
476 	fprintf (wdat_fp, "\t\tmax_harmonic = ptables[ti]->harmonics;\n");
477 	fprintf (wdat_fp, "\n");
478 	fprintf (wdat_fp, "\t\tfor ( ; li <= max_harmonic; li++)\n");
479 	fprintf (wdat_fp, "\t\t\tlookup[li] = ti;\n");
480 	fprintf (wdat_fp, "\t}\n");
481 	fprintf (wdat_fp, "\n");
482 	fprintf (wdat_fp, "\tfor ( ; li <= lookup_max; li++)\n");
483 	fprintf (wdat_fp, "\t\tlookup[li] = 0;\n");
484 	fprintf (wdat_fp, "}\n");
485 
486 	return 0;
487 }
488 
489 void
generate_sawtooth(LADSPA_Data * samples,unsigned long sample_count,unsigned long harmonics,float gibbs_comp)490 generate_sawtooth (LADSPA_Data * samples,
491                    unsigned long sample_count,
492                    unsigned long harmonics,
493                    float gibbs_comp)
494 {
495 	double phase_scale = 2.0 * M_PI / (double) sample_count;
496 	LADSPA_Data scale = 2.0f / M_PI;
497 	unsigned long i;
498 	unsigned long h;
499 	double mhf;
500 	double hf;
501 	double k;
502 	double m;
503 	double phase;
504 	double partial;
505 
506 	if (gibbs_comp > 0.0f)
507 	{
508 	/* Degree of Gibbs Effect compensation */
509 		mhf = (double) harmonics;
510 		k = M_PI * (double) gibbs_comp / mhf;
511 
512 		for (i = 0; i < sample_count; i++)
513 			samples[i] = 0.0f;
514 
515 		for (h = 1; h <= harmonics; h++)
516 		{
517 			hf = (double) h;
518 
519 		/* Gibbs Effect compensation - Hamming window */
520 		/* Modified slightly for smoother fade at highest frequencies */
521 			m = 0.54 + 0.46 * cos ((hf - 0.5 / mhf) * k);
522 
523 			for (i = 0; i < sample_count; i++)
524 			{
525 				phase = (double) i * phase_scale;
526 				partial = (m / hf) * sin (phase * hf);
527 				samples[i] += (LADSPA_Data) partial;
528 			}
529 		}
530 
531 		for (i = 0; i < sample_count; i++)
532 			samples[i] *= scale;
533 	}
534 	else
535 	{
536 	/* Allow overshoot */
537 		for (i = 0; i < sample_count; i++)
538 		{
539 			phase = (double) i * phase_scale;
540 			samples[i] = 0.0f;
541 
542 			for (h = 1; h <= harmonics; h++)
543 			{
544 				hf = (double) h;
545 				partial = (1.0 / hf) * sin (phase * hf);
546 				samples[i] += (LADSPA_Data) partial;
547 			}
548 			samples[i] *= scale;
549 		}
550 	}
551 }
552 
553 void
generate_square(LADSPA_Data * samples,unsigned long sample_count,unsigned long harmonics,float gibbs_comp)554 generate_square (LADSPA_Data * samples,
555                  unsigned long sample_count,
556                  unsigned long harmonics,
557                  float gibbs_comp)
558 {
559 	double phase_scale = 2.0 * M_PI / (double) sample_count;
560 	LADSPA_Data scale = 4.0f / M_PI;
561 	unsigned long i;
562 	unsigned long h;
563 	double mhf;
564 	double hf;
565 	double k;
566 	double m;
567 	double phase;
568 	double partial;
569 
570 	if (gibbs_comp > 0.0f)
571 	{
572 	/* Degree of Gibbs Effect compensation */
573 		mhf = (double) harmonics;
574 		k = M_PI * (double) gibbs_comp / mhf;
575 
576 		for (i = 0; i < sample_count; i++)
577 			samples[i] = 0.0f;
578 
579 		for (h = 1; h <= harmonics; h += 2)
580 		{
581 			hf = (double) h;
582 
583 		/* Gibbs Effect compensation - Hamming window */
584 		/* Modified slightly for smoother fade at highest frequencies */
585 			m = 0.54 + 0.46 * cos((hf - 0.5 / pow (mhf, 2.2)) * k);
586 
587 			for (i = 0; i < sample_count; i++)
588 			{
589 				phase = (double) i * phase_scale;
590 				partial = (m / hf) * sin (phase * hf);
591 				samples[i] += (LADSPA_Data) partial;
592 			}
593 		}
594 
595 		for (i = 0; i < sample_count; i++)
596 			samples[i] *= scale;
597 	}
598 	else
599 	{
600 	/* Allow overshoot */
601 		for (i = 0; i < sample_count; i++)
602 		{
603 			phase = (double) i * phase_scale;
604 			samples[i] = 0.0f;
605 
606 			for (h = 1; h <= harmonics; h += 2)
607 			{
608 				hf = (double) h;
609 				partial = (1.0 / hf) * sin (phase * hf);
610 				samples[i] += (LADSPA_Data) partial;
611 			}
612 			samples[i] *= scale;
613 		}
614 	}
615 }
616 
617 void
generate_parabola(LADSPA_Data * samples,unsigned long sample_count,unsigned long harmonics,float gibbs_comp)618 generate_parabola (LADSPA_Data * samples,
619                    unsigned long sample_count,
620                    unsigned long harmonics,
621                    float gibbs_comp)
622 {
623 	double phase_scale = 2.0 * M_PI / (double) sample_count;
624 	LADSPA_Data scale = 2.0f / (M_PI * M_PI);
625 	unsigned long i;
626 	unsigned long h;
627 	double mhf;
628 	double hf;
629 	double k;
630 	double m;
631 	double phase;
632 	double partial;
633     double sign;
634 
635 	if (gibbs_comp > 0.0f)
636 	{
637 	/* Degree of Gibbs Effect compensation */
638 		mhf = (double) harmonics;
639 		k = M_PI * (double) gibbs_comp / mhf;
640 
641 		for (i = 0; i < sample_count; i++)
642 			samples[i] = 0.0f;
643 
644 		sign = -1.0;
645 
646 		for (h = 1; h <= harmonics; h++)
647 		{
648 			hf = (double) h;
649 
650 		/* Gibbs Effect compensation - Hamming window */
651 		/* Modified slightly for smoother fade at highest frequencies */
652 			m = 0.54 + 0.46 * cos ((hf - 0.5 / mhf) * k);
653 
654 			for (i = 0; i < sample_count; i++)
655 			{
656 				phase = (double) i * phase_scale;
657 				partial = (sign * 4.0 / (hf * hf)) * cos (phase * hf);
658 				samples[i] += (LADSPA_Data) partial;
659 			}
660 			sign = -sign;
661 		}
662 
663 		for (i = 0; i < sample_count; i++)
664 			samples[i] *= scale;
665 	}
666 	else
667 	{
668 	/* Allow overshoot */
669 		for (i = 0; i < sample_count; i++)
670 		{
671 			phase = (double) i * phase_scale;
672 			samples[i] = 0.0f;
673 			sign = -1.0;
674 
675 			for (h = 1; h <= harmonics; h++)
676 			{
677 				hf = (double) h;
678 				partial = (sign * 4.0 / (hf * hf)) * cos (phase * hf);
679 				samples[i] += (LADSPA_Data) partial;
680 				sign = -sign;
681 			}
682 			samples[i] *= scale;
683 		}
684 	}
685 }
686