383af60dbe69742ac20a164b0fbcf0783fe87f97
[WebKit-https.git] / Source / WebCore / Modules / webaudio / PeriodicWave.cpp
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
40 const unsigned PeriodicWaveSize = 4096; // This must be a power of two.
41 const unsigned NumberOfRanges = 36; // There should be 3 * log2(PeriodicWaveSize) 1/3 octave ranges.
42 const float CentsPerRange = 1200 / 3; // 1/3 Octave.
43
44 namespace WebCore {
45     
46 using namespace VectorMath;
47
48 RefPtr<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
61 Ref<PeriodicWave> PeriodicWave::createSine(float sampleRate)
62 {
63     Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
64     waveTable->generateBasicWaveform(OscillatorNode::SINE);
65     return waveTable;
66 }
67
68 Ref<PeriodicWave> PeriodicWave::createSquare(float sampleRate)
69 {
70     Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
71     waveTable->generateBasicWaveform(OscillatorNode::SQUARE);
72     return waveTable;
73 }
74
75 Ref<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate)
76 {
77     Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
78     waveTable->generateBasicWaveform(OscillatorNode::SAWTOOTH);
79     return waveTable;
80 }
81
82 Ref<PeriodicWave> PeriodicWave::createTriangle(float sampleRate)
83 {
84     Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
85     waveTable->generateBasicWaveform(OscillatorNode::TRIANGLE);
86     return waveTable;
87 }
88
89 PeriodicWave::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
100 void 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
128 unsigned PeriodicWave::maxNumberOfPartials() const
129 {
130     return m_periodicWaveSize / 2;
131 }
132
133 unsigned 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.
150 void 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
220 void 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)