b85e0c9c2dbc9f27536ad6451d7122a81b398e37
[WebKit-https.git] / LayoutTests / webaudio / resources / biquad-testing.js
1 // Globals, to make testing and debugging easier.
2 var context;
3 var filter;
4 var signal;
5 var renderedBuffer;
6 var renderedData;
7
8 var sampleRate = 44100.0;
9 var pulseLengthFrames = .1 * sampleRate;
10
11 // Maximum allowed error for the test to succeed.  Experimentally determined. 
12 var maxAllowedError = 5.9e-8;
13
14 // This must be large enough so that the filtered result is
15 // essentially zero.  See comments for createTestAndRun.
16 var timeStep = .1;
17
18 // Maximum number of filters we can process (mostly for setting the
19 // render length correctly.)
20 var maxFilters = 5;
21
22 // How long to render.  Must be long enough for all of the filters we
23 // want to test.
24 var renderLengthSeconds = timeStep * (maxFilters + 1) ;
25
26 var renderLengthSamples = Math.round(renderLengthSeconds * sampleRate);
27
28 // Number of filters that will be processed.
29 var nFilters;
30
31 // A biquad filter has a z-transform of
32 // H(z) = (b0 + b1 / z + b2 / z^2) / (1 + a1 / z + a2 / z^2)
33 //
34 // The formulas for the various filters were taken from
35 // http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt.
36
37
38 // Lowpass filter.
39 function createLowpassFilter(freq, q, gain) {
40     var b0;
41     var b1;
42     var b2;
43     var a1;
44     var a2;
45
46     if (freq == 1) {
47         // The formula below works, except for roundoff.  When freq = 1,
48         // the filter is just a wire, so hardwire the coefficients.
49         b0 = 1;
50         b1 = 0;
51         b2 = 0;
52         a1 = 0;
53         a2 = 0;
54     } else {
55         var g = Math.pow(10, q / 20);
56         var d = Math.sqrt((4 - Math.sqrt(16 - 16 / (g * g))) / 2);
57         var theta = Math.PI * freq;
58         var sn = d * Math.sin(theta) / 2;
59         var beta = 0.5 * (1 - sn) / (1 + sn);
60         var gamma = (0.5 + beta) * Math.cos(theta);
61         var alpha = 0.25 * (0.5 + beta - gamma);
62
63         b0 = 2 * alpha;
64         b1 = 4 * alpha;
65         b2 = 2 * alpha;
66         a1 = 2 * (-gamma);
67         a2 = 2 * beta;
68     }
69
70     return {b0 : b0, b1 : b1, b2 : b2, a1 : a1, a2 : a2};
71 }
72
73 function createHighpassFilter(freq, q, gain) {
74     var b0;
75     var b1;
76     var b2;
77     var a1;
78     var a2;
79
80     if (freq == 1) {
81         // The filter is 0
82         b0 = 0;
83         b1 = 0;
84         b2 = 0;
85         a1 = 0;
86         a2 = 0;
87     } else if (freq == 0) {
88         // The filter is 1.  Computation of coefficients below is ok, but
89         // there's a pole at 1 and a zero at 1, so round-off could make
90         // the filter unstable.
91         b0 = 1;
92         b1 = 0;
93         b2 = 0;
94         a1 = 0;
95         a2 = 0;
96     } else {
97         var g = Math.pow(10, q / 20);
98         var d = Math.sqrt((4 - Math.sqrt(16 - 16 / (g * g))) / 2);
99         var theta = Math.PI * freq;
100         var sn = d * Math.sin(theta) / 2;
101         var beta = 0.5 * (1 - sn) / (1 + sn);
102         var gamma = (0.5 + beta) * Math.cos(theta);
103         var alpha = 0.25 * (0.5 + beta + gamma);
104
105         b0 = 2 * alpha;
106         b1 = -4 * alpha;
107         b2 = 2 * alpha;
108         a1 = 2 * (-gamma);
109         a2 = 2 * beta;
110     }
111
112     return {b0 : b0, b1 : b1, b2 : b2, a1 : a1, a2 : a2};
113 }
114
115 function normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2) {
116     var scale = 1 / a0;
117
118     return {b0 : b0 * scale,
119             b1 : b1 * scale,
120             b2 : b2 * scale,
121             a1 : a1 * scale,
122             a2 : a2 * scale};
123 }
124
125 function createBandpassFilter(freq, q, gain) {
126     var b0;
127     var b1;
128     var b2;
129     var a0;
130     var a1;
131     var a2;
132     var coef;
133
134     if (freq > 0 && freq < 1) {
135         var w0 = Math.PI * freq;
136         if (q > 0) {
137             var alpha = Math.sin(w0) / (2 * q);
138             var k = Math.cos(w0);
139
140             b0 = alpha;
141             b1 = 0;
142             b2 = -alpha;
143             a0 = 1 + alpha;
144             a1 = -2 * k;
145             a2 = 1 - alpha;
146
147             coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
148         } else {
149             // q = 0, and frequency is not 0 or 1.  The above formula has a
150             // divide by zero problem.  The limit of the z-transform as q
151             // approaches 0 is 1, so set the filter that way.
152             coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
153         }
154     } else {
155         // When freq = 0 or 1, the z-transform is identically 0,
156         // independent of q.
157         coef = {b0 : 0, b1 : 0, b2 : 0, a1 : 0, a2 : 0}
158     }
159   
160     return coef;
161 }
162
163 function createLowShelfFilter(freq, q, gain) {
164     // q not used
165     var b0;
166     var b1;
167     var b2;
168     var a0;
169     var a1;
170     var a2;
171     var coef;
172   
173     var S = 1;
174     var A = Math.pow(10, gain / 40);
175
176     if (freq == 1) {
177         // The filter is just a constant gain
178         coef = {b0 : A * A, b1 : 0, b2 : 0, a1 : 0, a2 : 0};  
179     } else if (freq == 0) {
180         // The filter is 1
181         coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};  
182     } else {
183         var w0 = Math.PI * freq;
184         var alpha = 1 / 2 * Math.sin(w0) * Math.sqrt((A + 1 / A) * (1 / S - 1) + 2);
185         var k = Math.cos(w0);
186         var k2 = 2 * Math.sqrt(A) * alpha;
187         var Ap1 = A + 1;
188         var Am1 = A - 1;
189
190         b0 = A * (Ap1 - Am1 * k + k2);
191         b1 = 2 * A * (Am1 - Ap1 * k);
192         b2 = A * (Ap1 - Am1 * k - k2);
193         a0 = Ap1 + Am1 * k + k2;
194         a1 = -2 * (Am1 + Ap1 * k);
195         a2 = Ap1 + Am1 * k - k2;
196         coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
197     }
198
199     return coef;
200 }
201
202 function createHighShelfFilter(freq, q, gain) {
203     // q not used
204     var b0;
205     var b1;
206     var b2;
207     var a0;
208     var a1;
209     var a2;
210     var coef;
211
212     var A = Math.pow(10, gain / 40);
213
214     if (freq == 1) {
215         // When freq = 1, the z-transform is 1
216         coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
217     } else if (freq > 0) {
218         var w0 = Math.PI * freq;
219         var S = 1;
220         var alpha = 0.5 * Math.sin(w0) * Math.sqrt((A + 1 / A) * (1 / S - 1) + 2);
221         var k = Math.cos(w0);
222         var k2 = 2 * Math.sqrt(A) * alpha;
223         var Ap1 = A + 1;
224         var Am1 = A - 1;
225
226         b0 = A * (Ap1 + Am1 * k + k2);
227         b1 = -2 * A * (Am1 + Ap1 * k);
228         b2 = A * (Ap1 + Am1 * k - k2);
229         a0 = Ap1 - Am1 * k + k2;
230         a1 = 2 * (Am1 - Ap1*k);
231         a2 = Ap1 - Am1 * k-k2;
232
233         coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
234     } else {
235         // When freq = 0, the filter is just a gain
236         coef = {b0 : A * A, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
237     }
238
239     return coef;
240 }
241
242 function createPeakingFilter(freq, q, gain) {
243     var b0;
244     var b1;
245     var b2;
246     var a0;
247     var a1;
248     var a2;
249     var coef;
250
251     var A = Math.pow(10, gain / 40);
252
253     if (freq > 0 && freq < 1) {
254         if (q > 0) {
255             var w0 = Math.PI * freq;
256             var alpha = Math.sin(w0) / (2 * q);
257             var k = Math.cos(w0);
258
259             b0 = 1 + alpha * A;
260             b1 = -2 * k;
261             b2 = 1 - alpha * A;
262             a0 = 1 + alpha / A;
263             a1 = -2 * k;
264             a2 = 1 - alpha / A;
265   
266             coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
267         } else {
268             // q = 0, we have a divide by zero problem in the formulas
269             // above.  But if we look at the z-transform, we see that the
270             // limit as q approaches 0 is A^2.
271             coef = {b0 : A * A, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
272         }
273     } else {
274         // freq = 0 or 1, the z-transform is 1
275         coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
276     }
277
278     return coef;
279 }
280
281 function createNotchFilter(freq, q, gain) {
282     var b0;
283     var b1;
284     var b2;
285     var a0;
286     var a1;
287     var a2;
288     var coef;
289
290     if (freq > 0 && freq < 1) {
291         if (q > 0) {
292             var w0 = Math.PI * freq;
293             var alpha = Math.sin(w0) / (2 * q);
294             var k = Math.cos(w0);
295
296             b0 = 1;
297             b1 = -2 * k;
298             b2 = 1;
299             a0 = 1 + alpha;
300             a1 = -2 * k;
301             a2 = 1 - alpha;
302             coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
303         } else {
304             // When q = 0, we get a divide by zero above.  The limit of the
305             // z-transform as q approaches 0 is 0, so set the coefficients
306             // appropriately.
307             coef = {b0 : 0, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
308         }
309     } else {
310         // When freq = 0 or 1, the z-transform is 1
311         coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
312     }
313
314     return coef;
315 }
316
317 function createAllpassFilter(freq, q, gain) {
318     var b0;
319     var b1;
320     var b2;
321     var a0;
322     var a1;
323     var a2;
324     var coef;
325
326     if (freq > 0 && freq < 1) {
327         if (q > 0) {
328             var w0 = Math.PI * freq;
329             var alpha = Math.sin(w0) / (2 * q);
330             var k = Math.cos(w0);
331
332             b0 = 1 - alpha;
333             b1 = -2 * k;
334             b2 = 1 + alpha;
335             a0 = 1 + alpha;
336             a1 = -2 * k;
337             a2 = 1 - alpha;
338             coef = normalizeFilterCoefficients(b0, b1, b2, a0, a1, a2);
339         } else {
340             // q = 0
341             coef = {b0 : -1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
342         }
343     } else {
344         coef = {b0 : 1, b1 : 0, b2 : 0, a1 : 0, a2 : 0};
345     }
346
347     return coef;
348 }
349
350 // Array of functions to compute the filter coefficients.  This must
351 // be arraned in the same order as the filter types in the idl file.
352 var filterCreatorFunction = [createLowpassFilter,
353                              createHighpassFilter,
354                              createBandpassFilter,
355                              createLowShelfFilter,
356                              createHighShelfFilter,
357                              createPeakingFilter,
358                              createNotchFilter,
359                              createAllpassFilter];
360
361 var filterTypeName = ["Lowpass filter",
362                       "Highpass filter",
363                       "Bandpass filter",
364                       "Lowshelf filter",
365                       "Highshelf filter",
366                       "Peaking filter",
367                       "Notch filter",
368                       "Allpass filter"];
369
370 function createFilter(filterType, freq, q, gain) {
371     return filterCreatorFunction[filterType](freq, q, gain);
372 }
373
374 function filterData(filterCoef, signal, len) {
375     var y = new Array(len);
376     var b0 = filterCoef.b0;
377     var b1 = filterCoef.b1;
378     var b2 = filterCoef.b2;
379     var a1 = filterCoef.a1;
380     var a2 = filterCoef. a2;
381
382     // Prime the pump. (Assumes the signal has length >= 2!)
383     y[0] = b0 * signal[0];
384     y[1] = b0 * signal[1] + b1 * signal[0] - a1 * y[0];
385
386     // Filter all of the signal that we have.
387     for (var k = 2; k < Math.min(signal.length, len); ++k) {
388         y[k] = b0 * signal[k] + b1 * signal[k-1] + b2 * signal[k-2] - a1 * y[k-1] - a2 * y[k-2];
389     }
390
391     // If we need to filter more, but don't have any signal left,
392     // assume the signal is zero.
393     for (var k = signal.length; k < len; ++k) {
394         y[k] = - a1 * y[k-1] - a2 * y[k-2];
395     }
396
397     return y;
398 }
399
400 function createImpulseBuffer(context, length) {
401     var impulse = context.createBuffer(1, length, context.sampleRate);
402     var data = impulse.getChannelData(0);
403     for (var k = 1; k < data.length; ++k) {
404         data[k] = 0;
405     }
406     data[0] = 1;
407
408     return impulse;
409 }
410
411
412 function createTestAndRun(context, filterType, filterParameters) {
413     // To test the filters, we apply a signal (an impulse) to each of
414     // the specified filters, with each signal starting at a different
415     // time.  The output of the filters is summed together at the
416     // output.  Thus for filter k, the signal input to the filter
417     // starts at time k * timeStep.  For this to work well, timeStep
418     // must be large enough for the output of each filter to have
419     // decayed to zero with timeStep seconds.  That way the filter
420     // outputs don't interfere with each other.
421     
422     nFilters = Math.min(filterParameters.length, maxFilters);
423
424     signal = new Array(nFilters);
425     filter = new Array(nFilters);
426
427     impulse = createImpulseBuffer(context, pulseLengthFrames);
428
429     // Create all of the signal sources and filters that we need.
430     for (var k = 0; k < nFilters; ++k) {
431         signal[k] = context.createBufferSource();
432         signal[k].buffer = impulse;
433
434         filter[k] = context.createBiquadFilter();
435         filter[k].type = filterType;
436         filter[k].frequency.value = context.sampleRate / 2 * filterParameters[k].cutoff;
437         filter[k].Q.value = filterParameters[k].q;
438         filter[k].gain.value = filterParameters[k].gain;
439
440         signal[k].connect(filter[k]);
441         filter[k].connect(context.destination);
442
443         signal[k].noteOn(timeStep * k);
444     }
445
446     context.oncomplete = checkFilterResponse(filterType, filterParameters);
447     context.startRendering();
448 }
449
450 function addSignal(dest, src, destOffset) {
451     // Add src to dest at the given dest offset.
452     for (var k = destOffset, j = 0; k < dest.length, j < src.length; ++k, ++j) {
453         dest[k] += src[j];
454     }
455 }
456
457 function generateReference(filterType, filterParameters) {
458     var result = new Array(renderLengthSamples);
459     var data = new Array(renderLengthSamples);
460     // Initialize the result array and data.
461     for (var k = 0; k < result.length; ++k) {
462         result[k] = 0;
463         data[k] = 0;
464     }
465     // Make data an impulse.
466     data[0] = 1;
467     
468     for (var k = 0; k < nFilters; ++k) {
469         // Filter an impulse
470         var filterCoef = createFilter(filterType,
471                                       filterParameters[k].cutoff,
472                                       filterParameters[k].q,
473                                       filterParameters[k].gain);
474         var y = filterData(filterCoef, data, renderLengthSamples);
475
476         // Accumulate this filtered data into the final output at the desired offset.
477         addSignal(result, y, timeToSampleFrame(timeStep * k, sampleRate));
478     }
479
480     return result;
481 }
482
483 // True if the number is not an infinity or NaN
484 function isValidNumber(x) {
485     return !isNaN(x) && (x != Infinity) && (x != -Infinity);
486 }
487
488 function checkFilterResponse(filterType, filterParameters) {
489     return function(event) {
490         renderedBuffer = event.renderedBuffer;
491         renderedData = renderedBuffer.getChannelData(0);
492
493         reference = generateReference(filterType, filterParameters);
494         
495         var len = Math.min(renderedData.length, reference.length);
496
497         var success = true;
498
499         // Maximum error between rendered data and expected data
500         var maxError = 0;
501
502         // Sample offset where the maximum error occurred.
503         var maxPosition = 0;
504
505         // Number of infinities or NaNs that occurred in the rendered data.
506         var invalidNumberCount = 0;
507
508         if (nFilters != filterParameters.length) {
509             testFailed("Test wanted " + filterParameters.length + " filters but only " + maxFilters + " allowed.");
510             success = false;
511         }
512
513         // Compare the rendered signal with our reference, keeping
514         // track of the maximum difference (and the offset of the max
515         // difference.)  Check for bad numbers in the rendered output
516         // too.  There shouldn't be any.
517         for (var k = 0; k < len; ++k) {
518             var err = Math.abs(renderedData[k] - reference[k]);
519             if (err > maxError) {
520                 maxError = err;
521                 maxPosition = k;
522             }
523             if (!isValidNumber(renderedData[k])) {
524                 ++invalidNumberCount;
525             }
526         }
527
528         if (invalidNumberCount > 0) {
529             testFailed("Rendered output has " + invalidNumberCount + " infinities or NaNs.");
530             success = false;
531         } else {
532             testPassed("Rendered output did not have infinities or NaNs.");
533         }
534         
535         if (maxError <= maxAllowedError) {
536             testPassed(filterTypeName[filterType] + " response is correct.");
537         } else {
538             testFailed(filterTypeName[filterType] + " response is incorrect.  Max err = " + maxError + " at " + maxPosition + ".  Threshold = " + maxAllowedError);
539             success = false;
540         }
541         
542         if (success) {
543             testPassed("Test signal was correctly filtered.");
544         } else {
545             testFailed("Test signal was not correctly filtered.");
546         }
547         finishJSTest();
548     }
549 }