1/*
2 * Copyright (C) 2012 Google Inc. All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. Neither the name of Apple Inc. ("Apple") nor the names of
14 * its contributors may be used to endorse or promote products derived
15 * from this software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
18 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20 * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
21 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29#include "config.h"
30
31#if ENABLE(WEB_AUDIO)
32
33#include "PeriodicWave.h"
34
35#include "FFTFrame.h"
36#include "VectorMath.h"
37#include <algorithm>
38
39const unsigned PeriodicWaveSize = 4096; // This must be a power of two.
40const unsigned NumberOfRanges = 36; // There should be 3 * log2(PeriodicWaveSize) 1/3 octave ranges.
41const float CentsPerRange = 1200 / 3; // 1/3 Octave.
42
43namespace WebCore {
44
45using namespace VectorMath;
46
47Ref<PeriodicWave> PeriodicWave::create(float sampleRate, Float32Array& real, Float32Array& imaginary)
48{
49 ASSERT(real.length() == imaginary.length());
50
51 auto waveTable = adoptRef(*new PeriodicWave(sampleRate));
52 waveTable->createBandLimitedTables(real.data(), imaginary.data(), real.length());
53 return waveTable;
54}
55
56Ref<PeriodicWave> PeriodicWave::createSine(float sampleRate)
57{
58 Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
59 waveTable->generateBasicWaveform(Type::Sine);
60 return waveTable;
61}
62
63Ref<PeriodicWave> PeriodicWave::createSquare(float sampleRate)
64{
65 Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
66 waveTable->generateBasicWaveform(Type::Square);
67 return waveTable;
68}
69
70Ref<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate)
71{
72 Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
73 waveTable->generateBasicWaveform(Type::Sawtooth);
74 return waveTable;
75}
76
77Ref<PeriodicWave> PeriodicWave::createTriangle(float sampleRate)
78{
79 Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
80 waveTable->generateBasicWaveform(Type::Triangle);
81 return waveTable;
82}
83
84PeriodicWave::PeriodicWave(float sampleRate)
85 : m_sampleRate(sampleRate)
86 , m_periodicWaveSize(PeriodicWaveSize)
87 , m_numberOfRanges(NumberOfRanges)
88 , m_centsPerRange(CentsPerRange)
89{
90 float nyquist = 0.5 * m_sampleRate;
91 m_lowestFundamentalFrequency = nyquist / maxNumberOfPartials();
92 m_rateScale = m_periodicWaveSize / m_sampleRate;
93}
94
95void PeriodicWave::waveDataForFundamentalFrequency(float fundamentalFrequency, float* &lowerWaveData, float* &higherWaveData, float& tableInterpolationFactor)
96{
97 // Negative frequencies are allowed, in which case we alias to the positive frequency.
98 fundamentalFrequency = fabsf(fundamentalFrequency);
99
100 // Calculate the pitch range.
101 float ratio = fundamentalFrequency > 0 ? fundamentalFrequency / m_lowestFundamentalFrequency : 0.5;
102 float centsAboveLowestFrequency = log2f(ratio) * 1200;
103
104 // Add one to round-up to the next range just in time to truncate partials before aliasing occurs.
105 float pitchRange = 1 + centsAboveLowestFrequency / m_centsPerRange;
106
107 pitchRange = std::max(pitchRange, 0.0f);
108 pitchRange = std::min(pitchRange, static_cast<float>(m_numberOfRanges - 1));
109
110 // The words "lower" and "higher" refer to the table data having the lower and higher numbers of partials.
111 // It's a little confusing since the range index gets larger the more partials we cull out.
112 // So the lower table data will have a larger range index.
113 unsigned rangeIndex1 = static_cast<unsigned>(pitchRange);
114 unsigned rangeIndex2 = rangeIndex1 < m_numberOfRanges - 1 ? rangeIndex1 + 1 : rangeIndex1;
115
116 lowerWaveData = m_bandLimitedTables[rangeIndex2]->data();
117 higherWaveData = m_bandLimitedTables[rangeIndex1]->data();
118
119 // Ranges from 0 -> 1 to interpolate between lower -> higher.
120 tableInterpolationFactor = pitchRange - rangeIndex1;
121}
122
123unsigned PeriodicWave::maxNumberOfPartials() const
124{
125 return m_periodicWaveSize / 2;
126}
127
128unsigned PeriodicWave::numberOfPartialsForRange(unsigned rangeIndex) const
129{
130 // Number of cents below nyquist where we cull partials.
131 float centsToCull = rangeIndex * m_centsPerRange;
132
133 // A value from 0 -> 1 representing what fraction of the partials to keep.
134 float cullingScale = pow(2, -centsToCull / 1200);
135
136 // The very top range will have all the partials culled.
137 unsigned numberOfPartials = cullingScale * maxNumberOfPartials();
138
139 return numberOfPartials;
140}
141
142// Convert into time-domain wave tables.
143// One table is created for each range for non-aliasing playback at different playback rates.
144// Thus, higher ranges have more high-frequency partials culled out.
145void PeriodicWave::createBandLimitedTables(const float* realData, const float* imagData, unsigned numberOfComponents)
146{
147 float normalizationScale = 1;
148
149 unsigned fftSize = m_periodicWaveSize;
150 unsigned halfSize = fftSize / 2;
151 unsigned i;
152
153 numberOfComponents = std::min(numberOfComponents, halfSize);
154
155 m_bandLimitedTables.reserveCapacity(m_numberOfRanges);
156
157 for (unsigned rangeIndex = 0; rangeIndex < m_numberOfRanges; ++rangeIndex) {
158 // This FFTFrame is used to cull partials (represented by frequency bins).
159 FFTFrame frame(fftSize);
160 float* realP = frame.realData();
161 float* imagP = frame.imagData();
162
163 // Copy from loaded frequency data and scale.
164 float scale = fftSize;
165 vsmul(realData, 1, &scale, realP, 1, numberOfComponents);
166 vsmul(imagData, 1, &scale, imagP, 1, numberOfComponents);
167
168 // If fewer components were provided than 1/2 FFT size, then clear the remaining bins.
169 for (i = numberOfComponents; i < halfSize; ++i) {
170 realP[i] = 0;
171 imagP[i] = 0;
172 }
173
174 // Generate complex conjugate because of the way the inverse FFT is defined.
175 float minusOne = -1;
176 vsmul(imagP, 1, &minusOne, imagP, 1, halfSize);
177
178 // Find the starting bin where we should start culling.
179 // We need to clear out the highest frequencies to band-limit the waveform.
180 unsigned numberOfPartials = numberOfPartialsForRange(rangeIndex);
181
182 // Cull the aliasing partials for this pitch range.
183 for (i = numberOfPartials + 1; i < halfSize; ++i) {
184 realP[i] = 0;
185 imagP[i] = 0;
186 }
187 // Clear packed-nyquist if necessary.
188 if (numberOfPartials < halfSize)
189 imagP[0] = 0;
190
191 // Clear any DC-offset.
192 realP[0] = 0;
193
194 // Create the band-limited table.
195 m_bandLimitedTables.append(std::make_unique<AudioFloatArray>(m_periodicWaveSize));
196
197 // Apply an inverse FFT to generate the time-domain table data.
198 float* data = m_bandLimitedTables[rangeIndex]->data();
199 frame.doInverseFFT(data);
200
201 // For the first range (which has the highest power), calculate its peak value then compute normalization scale.
202 if (!rangeIndex) {
203 float maxValue;
204 vmaxmgv(data, 1, &maxValue, m_periodicWaveSize);
205
206 if (maxValue)
207 normalizationScale = 1.0f / maxValue;
208 }
209
210 // Apply normalization scale.
211 vsmul(data, 1, &normalizationScale, data, 1, m_periodicWaveSize);
212 }
213}
214
215void PeriodicWave::generateBasicWaveform(Type shape)
216{
217 unsigned fftSize = periodicWaveSize();
218 unsigned halfSize = fftSize / 2;
219
220 AudioFloatArray real(halfSize);
221 AudioFloatArray imag(halfSize);
222 float* realP = real.data();
223 float* imagP = imag.data();
224
225 // Clear DC and Nyquist.
226 realP[0] = 0;
227 imagP[0] = 0;
228
229 for (unsigned n = 1; n < halfSize; ++n) {
230 float omega = 2 * piFloat * n;
231 float invOmega = 1 / omega;
232
233 // Fourier coefficients according to standard definition.
234 float a; // Coefficient for cos().
235 float b; // Coefficient for sin().
236
237 // Calculate Fourier coefficients depending on the shape.
238 // Note that the overall scaling (magnitude) of the waveforms is normalized in createBandLimitedTables().
239 switch (shape) {
240 case Type::Sine:
241 // Standard sine wave function.
242 a = 0;
243 b = (n == 1) ? 1 : 0;
244 break;
245 case Type::Square:
246 // Square-shaped waveform with the first half its maximum value and the second half its minimum value.
247 a = 0;
248 b = invOmega * ((n & 1) ? 2 : 0);
249 break;
250 case Type::Sawtooth:
251 // Sawtooth-shaped waveform with the first half ramping from zero to maximum and the second half from minimum to zero.
252 a = 0;
253 b = -invOmega * cos(0.5 * omega);
254 break;
255 case Type::Triangle:
256 // Triangle-shaped waveform going from its maximum value to its minimum value then back to the maximum value.
257 a = (4 - 4 * cos(0.5 * omega)) / (n * n * piFloat * piFloat);
258 b = 0;
259 break;
260 }
261
262 realP[n] = a;
263 imagP[n] = b;
264 }
265
266 createBandLimitedTables(realP, imagP, halfSize);
267}
268
269} // namespace WebCore
270
271#endif // ENABLE(WEB_AUDIO)
272