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 "OscillatorNode.h"
37#include "VectorMath.h"
38#include <algorithm>
39
40const unsigned PeriodicWaveSize = 4096; // This must be a power of two.
41const unsigned NumberOfRanges = 36; // There should be 3 * log2(PeriodicWaveSize) 1/3 octave ranges.
42const float CentsPerRange = 1200 / 3; // 1/3 Octave.
43
44namespace WebCore {
45
46using namespace VectorMath;
47
48PassRefPtr<PeriodicWave> PeriodicWave::create(float sampleRate, Float32Array* real, Float32Array* imag)
49{
50    bool isGood = real && imag && real->length() == imag->length();
51    ASSERT(isGood);
52    if (isGood) {
53        RefPtr<PeriodicWave> waveTable = adoptRef(new PeriodicWave(sampleRate));
54        size_t numberOfComponents = real->length();
55        waveTable->createBandLimitedTables(real->data(), imag->data(), numberOfComponents);
56        return waveTable;
57    }
58    return nullptr;
59}
60
61PassRefPtr<PeriodicWave> PeriodicWave::createSine(float sampleRate)
62{
63    RefPtr<PeriodicWave> waveTable = adoptRef(new PeriodicWave(sampleRate));
64    waveTable->generateBasicWaveform(OscillatorNode::SINE);
65    return waveTable;
66}
67
68PassRefPtr<PeriodicWave> PeriodicWave::createSquare(float sampleRate)
69{
70    RefPtr<PeriodicWave> waveTable = adoptRef(new PeriodicWave(sampleRate));
71    waveTable->generateBasicWaveform(OscillatorNode::SQUARE);
72    return waveTable;
73}
74
75PassRefPtr<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate)
76{
77    RefPtr<PeriodicWave> waveTable = adoptRef(new PeriodicWave(sampleRate));
78    waveTable->generateBasicWaveform(OscillatorNode::SAWTOOTH);
79    return waveTable;
80}
81
82PassRefPtr<PeriodicWave> PeriodicWave::createTriangle(float sampleRate)
83{
84    RefPtr<PeriodicWave> waveTable = adoptRef(new PeriodicWave(sampleRate));
85    waveTable->generateBasicWaveform(OscillatorNode::TRIANGLE);
86    return waveTable;
87}
88
89PeriodicWave::PeriodicWave(float sampleRate)
90    : m_sampleRate(sampleRate)
91    , m_periodicWaveSize(PeriodicWaveSize)
92    , m_numberOfRanges(NumberOfRanges)
93    , m_centsPerRange(CentsPerRange)
94{
95    float nyquist = 0.5 * m_sampleRate;
96    m_lowestFundamentalFrequency = nyquist / maxNumberOfPartials();
97    m_rateScale = m_periodicWaveSize / m_sampleRate;
98}
99
100void PeriodicWave::waveDataForFundamentalFrequency(float fundamentalFrequency, float* &lowerWaveData, float* &higherWaveData, float& tableInterpolationFactor)
101{
102    // Negative frequencies are allowed, in which case we alias to the positive frequency.
103    fundamentalFrequency = fabsf(fundamentalFrequency);
104
105    // Calculate the pitch range.
106    float ratio = fundamentalFrequency > 0 ? fundamentalFrequency / m_lowestFundamentalFrequency : 0.5;
107    float centsAboveLowestFrequency = log2f(ratio) * 1200;
108
109    // Add one to round-up to the next range just in time to truncate partials before aliasing occurs.
110    float pitchRange = 1 + centsAboveLowestFrequency / m_centsPerRange;
111
112    pitchRange = std::max(pitchRange, 0.0f);
113    pitchRange = std::min(pitchRange, static_cast<float>(m_numberOfRanges - 1));
114
115    // The words "lower" and "higher" refer to the table data having the lower and higher numbers of partials.
116    // It's a little confusing since the range index gets larger the more partials we cull out.
117    // So the lower table data will have a larger range index.
118    unsigned rangeIndex1 = static_cast<unsigned>(pitchRange);
119    unsigned rangeIndex2 = rangeIndex1 < m_numberOfRanges - 1 ? rangeIndex1 + 1 : rangeIndex1;
120
121    lowerWaveData = m_bandLimitedTables[rangeIndex2]->data();
122    higherWaveData = m_bandLimitedTables[rangeIndex1]->data();
123
124    // Ranges from 0 -> 1 to interpolate between lower -> higher.
125    tableInterpolationFactor = pitchRange - rangeIndex1;
126}
127
128unsigned PeriodicWave::maxNumberOfPartials() const
129{
130    return m_periodicWaveSize / 2;
131}
132
133unsigned PeriodicWave::numberOfPartialsForRange(unsigned rangeIndex) const
134{
135    // Number of cents below nyquist where we cull partials.
136    float centsToCull = rangeIndex * m_centsPerRange;
137
138    // A value from 0 -> 1 representing what fraction of the partials to keep.
139    float cullingScale = pow(2, -centsToCull / 1200);
140
141    // The very top range will have all the partials culled.
142    unsigned numberOfPartials = cullingScale * maxNumberOfPartials();
143
144    return numberOfPartials;
145}
146
147// Convert into time-domain wave tables.
148// One table is created for each range for non-aliasing playback at different playback rates.
149// Thus, higher ranges have more high-frequency partials culled out.
150void PeriodicWave::createBandLimitedTables(const float* realData, const float* imagData, unsigned numberOfComponents)
151{
152    float normalizationScale = 1;
153
154    unsigned fftSize = m_periodicWaveSize;
155    unsigned halfSize = fftSize / 2;
156    unsigned i;
157
158    numberOfComponents = std::min(numberOfComponents, halfSize);
159
160    m_bandLimitedTables.reserveCapacity(m_numberOfRanges);
161
162    for (unsigned rangeIndex = 0; rangeIndex < m_numberOfRanges; ++rangeIndex) {
163        // This FFTFrame is used to cull partials (represented by frequency bins).
164        FFTFrame frame(fftSize);
165        float* realP = frame.realData();
166        float* imagP = frame.imagData();
167
168        // Copy from loaded frequency data and scale.
169        float scale = fftSize;
170        vsmul(realData, 1, &scale, realP, 1, numberOfComponents);
171        vsmul(imagData, 1, &scale, imagP, 1, numberOfComponents);
172
173        // If fewer components were provided than 1/2 FFT size, then clear the remaining bins.
174        for (i = numberOfComponents; i < halfSize; ++i) {
175            realP[i] = 0;
176            imagP[i] = 0;
177        }
178
179        // Generate complex conjugate because of the way the inverse FFT is defined.
180        float minusOne = -1;
181        vsmul(imagP, 1, &minusOne, imagP, 1, halfSize);
182
183        // Find the starting bin where we should start culling.
184        // We need to clear out the highest frequencies to band-limit the waveform.
185        unsigned numberOfPartials = numberOfPartialsForRange(rangeIndex);
186
187        // Cull the aliasing partials for this pitch range.
188        for (i = numberOfPartials + 1; i < halfSize; ++i) {
189            realP[i] = 0;
190            imagP[i] = 0;
191        }
192        // Clear packed-nyquist if necessary.
193        if (numberOfPartials < halfSize)
194            imagP[0] = 0;
195
196        // Clear any DC-offset.
197        realP[0] = 0;
198
199        // Create the band-limited table.
200        m_bandLimitedTables.append(std::make_unique<AudioFloatArray>(m_periodicWaveSize));
201
202        // Apply an inverse FFT to generate the time-domain table data.
203        float* data = m_bandLimitedTables[rangeIndex]->data();
204        frame.doInverseFFT(data);
205
206        // For the first range (which has the highest power), calculate its peak value then compute normalization scale.
207        if (!rangeIndex) {
208            float maxValue;
209            vmaxmgv(data, 1, &maxValue, m_periodicWaveSize);
210
211            if (maxValue)
212                normalizationScale = 1.0f / maxValue;
213        }
214
215        // Apply normalization scale.
216        vsmul(data, 1, &normalizationScale, data, 1, m_periodicWaveSize);
217    }
218}
219
220void PeriodicWave::generateBasicWaveform(int shape)
221{
222    unsigned fftSize = periodicWaveSize();
223    unsigned halfSize = fftSize / 2;
224
225    AudioFloatArray real(halfSize);
226    AudioFloatArray imag(halfSize);
227    float* realP = real.data();
228    float* imagP = imag.data();
229
230    // Clear DC and Nyquist.
231    realP[0] = 0;
232    imagP[0] = 0;
233
234    for (unsigned n = 1; n < halfSize; ++n) {
235        float omega = 2 * piFloat * n;
236        float invOmega = 1 / omega;
237
238        // Fourier coefficients according to standard definition.
239        float a; // Coefficient for cos().
240        float b; // Coefficient for sin().
241
242        // Calculate Fourier coefficients depending on the shape.
243        // Note that the overall scaling (magnitude) of the waveforms is normalized in createBandLimitedTables().
244        switch (shape) {
245        case OscillatorNode::SINE:
246            // Standard sine wave function.
247            a = 0;
248            b = (n == 1) ? 1 : 0;
249            break;
250        case OscillatorNode::SQUARE:
251            // Square-shaped waveform with the first half its maximum value and the second half its minimum value.
252            a = 0;
253            b = invOmega * ((n & 1) ? 2 : 0);
254            break;
255        case OscillatorNode::SAWTOOTH:
256            // Sawtooth-shaped waveform with the first half ramping from zero to maximum and the second half from minimum to zero.
257            a = 0;
258            b = -invOmega * cos(0.5 * omega);
259            break;
260        case OscillatorNode::TRIANGLE:
261            // Triangle-shaped waveform going from its maximum value to its minimum value then back to the maximum value.
262            a = (4 - 4 * cos(0.5 * omega)) / (n * n * piFloat * piFloat);
263            b = 0;
264            break;
265        default:
266            ASSERT_NOT_REACHED();
267            a = 0;
268            b = 0;
269            break;
270        }
271
272        realP[n] = a;
273        imagP[n] = b;
274    }
275
276    createBandLimitedTables(realP, imagP, halfSize);
277}
278
279} // namespace WebCore
280
281#endif // ENABLE(WEB_AUDIO)
282