Commit 65d50038 authored by Raymond Toy's avatar Raymond Toy Committed by Commit Bot

Better interpolation of oscillator at low frequencies

For low frequencies, use a 3-point or 5-point Lagrange interpolation
scheme to improve the accuracy of the oscillator instead of using
simple linear interpolation.

For a typical context, the linear interpolator is used if the
frequency >= 3.2 Hz; the 3-point, for 3.2 > frequency > 1.7 Hz; and
the 5-point for all other cases.  Thus, typical use-cases are not
impaired by this change.

Some benchmark results from running an oscillator for 1000 seconds in an
offline context.  This is the time measured from startRendering()
until the promise is delivered.  The SNR is the measured SNR between
the test signal and a reference 440 Hz signal.  The linear time was
the same before and after this CL.

linear:
time: 1844.6 ms, stddev: 154.4
SNR: 27.7 dB

3-point:
time: 2426.5 ms, stddev: 265.0
SNR: 37.2 dB

5-point:
time: 2752.6 ms, stddev: 215.1
SNR: 59.3 dB

7-point (for reference, not used now):
time: 3239.4 ms, stddev: 263.8
SNR: 82.5 dB

Bug: 773507
Test: Oscillator/osc-440hz.html
Change-Id: Ib49fab4fb27738f4f02bd41f8adcc2e3b2dd2137
Reviewed-on: https://chromium-review.googlesource.com/713576
Commit-Queue: Raymond Toy <rtoy@chromium.org>
Reviewed-by: default avatarHongchan Choi <hongchan@chromium.org>
Cr-Commit-Position: refs/heads/master@{#510133}
parent 98422f9b
......@@ -38,7 +38,7 @@
minDecibels: -50,
floatRelError: 9.6549e-7,
},
{order: 6, floatRelError: 7.3494e-6},
{order: 6, floatRelError: 1.0084e-5},
{order: 7, floatRelError: 1.5645e-6},
{order: 8, floatRelError: 8.4828e-7},
{order: 9, floatRelError: 2.3906e-5},
......
......@@ -56,8 +56,8 @@
// experiments.
compareBuffersWithConstraints(should, actual, expected, {
prefix: '',
thresholdSNR: 93.31,
thresholdDiffULP: 1.01,
thresholdSNR: 93.337,
thresholdDiffULP: 1.0141,
thresholdDiffCount: 0,
bitDepth: 16
});
......
......@@ -55,8 +55,8 @@
// experiments.
compareBuffersWithConstraints(should, actual, expected, {
prefix: 'Verify',
thresholdSNR: 92.72,
thresholdDiffULP: 0.985,
thresholdSNR: 93.274,
thresholdDiffULP: 1.0519,
thresholdDiffCount: 0,
bitDepth: 16
});
......
<!DOCTYPE html>
<html>
<head>
<title>
Test Custom 440 Hz Oscillator
</title>
<script src="../../resources/testharness.js"></script>
<script src="../../resources/testharnessreport.js"></script>
<script src="../resources/audit-util.js"></script>
<script src="../resources/audit.js"></script>
</head>
<body>
<script id = 'layout-test-code'>let sampleRate = 48000;
// Minimum allowed SNR between the actual oscillator and the expected
// result. Experimentally determined.
let snrThreshold = 59.287;
// Max absolute difference between actual and expected oscillator
// outputs. Experimentally determined.
let maxDiffThreshold = 1.5632e-3;
let audit = Audit.createTaskRunner();
// Test that two ways of creating a 440-Hz sine wave produce the same
// results, within some tolerance.
audit.define('440 Hz Osc', (task, should) => {
// 3-channel context: channel 0 = actual output; channel 1 = expected
// output, channel 2 = difference between channels 0 and 1.
let context = new OfflineAudioContext(3, sampleRate, sampleRate);
let merger = new ChannelMergerNode(
context, {numberOfInputs: context.destination.channelCount});
merger.connect(context.destination);
// Create a 440 Hz oscillator in following way. Use a PeriodicWave with
// a single harmonic value of 440 and use that in an oscillator with
// fundamental frequeqncy of 1 Hz. This should produce a 440 Hz sine
// wave. This is the test signal.
let imag = new Float32Array(512);
imag[440] = 1;
let wave440 =
new PeriodicWave(context, {imag: imag, disableNormaliztion: true});
let oscTest =
new OscillatorNode(context, {frequency: 1, periodicWave: wave440});
// Create a 440 Hz oscillator by specifying a frequency of 440 Hz. We
// use a periodicWave for this so we can disable normalization.
let oscRef = new OscillatorNode(context, {
frequency: 440,
periodocWave: new PeriodicWave(
context, {imag: [0, 1], disableNormaliztion: true})
});
oscTest.connect(merger, 0, 0);
oscRef.connect(merger, 0, 1);
// Compute the difference of the two signals
let neg = new GainNode(context, {gain: -1});
neg.connect(merger, 0, 2);
oscTest.connect(merger, 0, 2);
oscRef.connect(neg);
oscTest.start();
oscRef.start();
context.startRendering()
.then(result => {
// Extract the actual output, the reference output and the
// difference signal from the rendered data.
actual = result.getChannelData(0);
ref = result.getChannelData(1);
diff = result.getChannelData(2);
// The actual and reference signals should be close, and the SNR
// should be high.
should(actual, 'Test oscillator output').beCloseToArray(
ref, {absoluteThreshold: maxDiffThreshold});
let snr = 10 * Math.log10(computeSNR(actual, ref));
should(snr, 'SNR').beGreaterThanOrEqualTo(snrThreshold);
})
.then(() => task.done());
});
audit.run();
</script>
</body>
</html>
......@@ -21,7 +21,7 @@
1, tester.sampleRate * tester.lengthInSeconds, tester.sampleRate);
// The thresholds are experimentally determined.
tester.setThresholds({snr: 140.073, maxDiff: 1.2815e-6});
tester.setThresholds({snr: 139.991, maxDiff: 1.2815e-6});
tester.runTest(
context, 'custom', 'Custom Oscillator with Exponential Sweep', task,
should);
......
......@@ -21,7 +21,7 @@
1, tester.sampleRate * tester.lengthInSeconds, tester.sampleRate);
// The thresholds are experimentally determined.
tester.setThresholds({snr: 140.51, maxDiff: 1.5386e-6});
tester.setThresholds({snr: 140.51, maxDiff: 1.5423e-6});
tester.runTest(
context, 'triangle', 'Triangle Oscillator with Exponential Sweep ',
task, should);
......
......@@ -229,6 +229,110 @@ bool OscillatorHandler::CalculateSampleAccuratePhaseIncrements(
return has_sample_accurate_values;
}
static float DoInterpolation(double virtual_read_index,
float incr,
unsigned read_index_mask,
float table_interpolation_factor,
const float* lower_wave_data,
const float* higher_wave_data) {
DCHECK_GE(incr, 0);
double sample_lower = 0;
double sample_higher = 0;
unsigned read_index_0 = static_cast<unsigned>(virtual_read_index);
// Consider a typical sample rate of 44100 Hz and max periodic wave
// size of 4096. The relationship between |incr| and the frequency
// of the oscillator is |incr| = freq * 4096/44100. Or freq =
// |incr|*44100/4096 = 10.8*|incr|.
//
// For the |incr| thresholds below, this means that we use linear
// interpolation for all freq >= 3.2 Hz, 3-point Lagrange
// for freq >= 1.7 Hz and 5-point Lagrange for every thing else.
//
// We use Lagrange interpolation because it's relatively simple to
// implement and fairly inexpensive, and the interpolator always
// passes through known points.
if (incr >= 0.3) {
// Increment is fairly large, so we're doing no more than about 3
// points between each wave table entry. Assume linear
// interpolation between points is good enough.
unsigned read_index2 = read_index_0 + 1;
// Contain within valid range.
read_index_0 = read_index_0 & read_index_mask;
read_index2 = read_index2 & read_index_mask;
float sample1_lower = lower_wave_data[read_index_0];
float sample2_lower = lower_wave_data[read_index2];
float sample1_higher = higher_wave_data[read_index_0];
float sample2_higher = higher_wave_data[read_index2];
// Linearly interpolate within each table (lower and higher).
double interpolation_factor =
static_cast<float>(virtual_read_index) - read_index_0;
sample_higher = (1 - interpolation_factor) * sample1_higher +
interpolation_factor * sample2_higher;
sample_lower = (1 - interpolation_factor) * sample1_lower +
interpolation_factor * sample2_lower;
} else if (incr >= .16) {
// We're doing about 6 interpolation values between each wave
// table sample. Just use a 3-point Lagrange interpolator to get a
// better estimate than just linear.
//
// See 3-point formula in http://dlmf.nist.gov/3.3#ii
unsigned read_index[3];
for (int k = -1; k <= 1; ++k) {
read_index[k + 1] = (read_index_0 + k) & read_index_mask;
}
double a[3];
double t = virtual_read_index - read_index_0;
a[0] = 0.5 * t * (t - 1);
a[1] = 1 - t * t;
a[2] = 0.5 * t * (t + 1);
for (int k = 0; k < 3; ++k) {
sample_lower += a[k] * lower_wave_data[read_index[k]];
sample_higher += a[k] * higher_wave_data[read_index[k]];
}
} else {
// For everything else (more than 6 points per entry), we'll do a
// 5-point Lagrange interpolator. This is a trade-off between
// quality and speed.
//
// See 5-point formula in http://dlmf.nist.gov/3.3#ii
unsigned read_index[5];
for (int k = -2; k <= 2; ++k) {
read_index[k + 2] = (read_index_0 + k) & read_index_mask;
}
double a[5];
double t = virtual_read_index - read_index_0;
double t2 = t * t;
a[0] = t * (t2 - 1) * (t - 2) / 24;
a[1] = -t * (t - 1) * (t2 - 4) / 6;
a[2] = (t2 - 1) * (t2 - 4) / 4;
a[3] = -t * (t + 1) * (t2 - 4) / 6;
a[4] = t * (t2 - 1) * (t + 2) / 24;
for (int k = 0; k < 5; ++k) {
sample_lower += a[k] * lower_wave_data[read_index[k]];
sample_higher += a[k] * higher_wave_data[read_index[k]];
}
}
// Then interpolate between the two tables.
float sample = (1 - table_interpolation_factor) * sample_higher +
table_interpolation_factor * sample_lower;
return sample;
}
void OscillatorHandler::Process(size_t frames_to_process) {
AudioBus* output_bus = Output(0).Bus();
......@@ -321,13 +425,6 @@ void OscillatorHandler::Process(size_t frames_to_process) {
}
while (n--) {
unsigned read_index = static_cast<unsigned>(virtual_read_index);
unsigned read_index2 = read_index + 1;
// Contain within valid range.
read_index = read_index & read_index_mask;
read_index2 = read_index2 & read_index_mask;
if (has_sample_accurate_values) {
incr = *phase_increments++;
......@@ -337,22 +434,9 @@ void OscillatorHandler::Process(size_t frames_to_process) {
table_interpolation_factor);
}
float sample1_lower = lower_wave_data[read_index];
float sample2_lower = lower_wave_data[read_index2];
float sample1_higher = higher_wave_data[read_index];
float sample2_higher = higher_wave_data[read_index2];
// Linearly interpolate within each table (lower and higher).
float interpolation_factor =
static_cast<float>(virtual_read_index) - read_index;
float sample_higher = (1 - interpolation_factor) * sample1_higher +
interpolation_factor * sample2_higher;
float sample_lower = (1 - interpolation_factor) * sample1_lower +
interpolation_factor * sample2_lower;
// Then interpolate between the two tables.
float sample = (1 - table_interpolation_factor) * sample_higher +
table_interpolation_factor * sample_lower;
float sample = DoInterpolation(virtual_read_index, fabs(incr),
read_index_mask, table_interpolation_factor,
lower_wave_data, higher_wave_data);
*dest_p++ = sample;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment