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 | |
39 | const unsigned PeriodicWaveSize = 4096; // This must be a power of two. |
40 | const unsigned NumberOfRanges = 36; // There should be 3 * log2(PeriodicWaveSize) 1/3 octave ranges. |
41 | const float CentsPerRange = 1200 / 3; // 1/3 Octave. |
42 | |
43 | namespace WebCore { |
44 | |
45 | using namespace VectorMath; |
46 | |
47 | Ref<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 | |
56 | Ref<PeriodicWave> PeriodicWave::createSine(float sampleRate) |
57 | { |
58 | Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate)); |
59 | waveTable->generateBasicWaveform(Type::Sine); |
60 | return waveTable; |
61 | } |
62 | |
63 | Ref<PeriodicWave> PeriodicWave::createSquare(float sampleRate) |
64 | { |
65 | Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate)); |
66 | waveTable->generateBasicWaveform(Type::Square); |
67 | return waveTable; |
68 | } |
69 | |
70 | Ref<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate) |
71 | { |
72 | Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate)); |
73 | waveTable->generateBasicWaveform(Type::Sawtooth); |
74 | return waveTable; |
75 | } |
76 | |
77 | Ref<PeriodicWave> PeriodicWave::createTriangle(float sampleRate) |
78 | { |
79 | Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate)); |
80 | waveTable->generateBasicWaveform(Type::Triangle); |
81 | return waveTable; |
82 | } |
83 | |
84 | PeriodicWave::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 | |
95 | void 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 | |
123 | unsigned PeriodicWave::maxNumberOfPartials() const |
124 | { |
125 | return m_periodicWaveSize / 2; |
126 | } |
127 | |
128 | unsigned 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. |
145 | void 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 | |
215 | void 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 | |