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