src/third_party/WebKit/Source/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 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 blink {
  45
  46 using namespace VectorMath;
  47
  48 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         PeriodicWave* periodicWave = 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 PeriodicWave* PeriodicWave::createSine(float sampleRate)
  62 {
  63     PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
  64     periodicWave->generateBasicWaveform(OscillatorNode::SINE);
  65     return periodicWave;
  66 }
  67
  68 PeriodicWave* PeriodicWave::createSquare(float sampleRate)
  69 {
  70     PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
  71     periodicWave->generateBasicWaveform(OscillatorNode::SQUARE);
  72     return periodicWave;
  73 }
  74
  75 PeriodicWave* PeriodicWave::createSawtooth(float sampleRate)
  76 {
  77     PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
  78     periodicWave->generateBasicWaveform(OscillatorNode::SAWTOOTH);
  79     return periodicWave;
  80 }
  81
  82 PeriodicWave* PeriodicWave::createTriangle(float sampleRate)
  83 {
  84     PeriodicWave* periodicWave = 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     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 buffers.
 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         OwnPtr<AudioFloatArray> table = adoptPtr(new AudioFloatArray(m_periodicWaveSize));
 201         m_bandLimitedTables.append(table.release());
 202
 203         // Apply an inverse FFT to generate the time-domain table data.
 204         float* data = m_bandLimitedTables[rangeIndex]->data();
 205         frame.doInverseFFT(data);
 206
 207         // For the first range (which has the highest power), calculate its peak value then compute normalization scale.
 208         if (!rangeIndex) {
 209             float maxValue;
 210             vmaxmgv(data, 1, &maxValue, m_periodicWaveSize);
 211
 212             if (maxValue)
 213                 normalizationScale = 1.0f / maxValue;
 214         }
 215
 216         // Apply normalization scale.
 217         vsmul(data, 1, &normalizationScale, data, 1, m_periodicWaveSize);
 218     }
 219 }
 220
 221 void PeriodicWave::generateBasicWaveform(int shape)
 222 {
 223     unsigned fftSize = periodicWaveSize();
 224     unsigned halfSize = fftSize / 2;
 225
 226     AudioFloatArray real(halfSize);
 227     AudioFloatArray imag(halfSize);
 228     float* realP = real.data();
 229     float* imagP = imag.data();
 230
 231     // Clear DC and Nyquist.
 232     realP[0] = 0;
 233     imagP[0] = 0;
 234
 235     for (unsigned n = 1; n < halfSize; ++n) {
 236         float piFactor = 2 / (n * piFloat);
 237
 238         // All waveforms are odd functions with a positive slope at time 0. Hence the coefficients
 239         // for cos() are always 0.
 240
 241         // Fourier coefficients according to standard definition:
 242         // b = 1/pi*integrate(f(x)*sin(n*x), x, -pi, pi)
 243         //   = 2/pi*integrate(f(x)*sin(n*x), x, 0, pi)
 244         // since f(x) is an odd function.
 245
 246         float b; // Coefficient for sin().
 247
 248         // Calculate Fourier coefficients depending on the shape. Note that the overall scaling
 249         // (magnitude) of the waveforms is normalized in createBandLimitedTables().
 250         switch (shape) {
 251         case OscillatorNode::SINE:
 252             // Standard sine wave function.
 253             b = (n == 1) ? 1 : 0;
 254             break;
 255         case OscillatorNode::SQUARE:
 256             // Square-shaped waveform with the first half its maximum value and the second half its
 257             // minimum value.
 258             //
 259             // See http://mathworld.wolfram.com/FourierSeriesSquareWave.html
 260             //
 261             // b[n] = 2/n/pi*(1-(-1)^n)
 262             //      = 4/n/pi for n odd and 0 otherwise.
 263             //      = 2*(2/(n*pi)) for n odd
 264             b = (n & 1) ? 2 * piFactor : 0;
 265             break;
 266         case OscillatorNode::SAWTOOTH:
 267             // Sawtooth-shaped waveform with the first half ramping from zero to maximum and the
 268             // second half from minimum to zero.
 269             //
 270             // b[n] = -2*(-1)^n/pi/n
 271             //      = (2/(n*pi))*(-1)^(n+1)
 272             b = piFactor * ((n & 1) ? 1 : -1);
 273             break;
 274         case OscillatorNode::TRIANGLE:
 275             // Triangle-shaped waveform going from 0 at time 0 to 1 at time pi/2 and back to 0 at
 276             // time pi.
 277             //
 278             // See http://mathworld.wolfram.com/FourierSeriesTriangleWave.html
 279             //
 280             // b[n] = 8*sin(pi*k/2)/(pi*k)^2
 281             //      = 8/pi^2/n^2*(-1)^((n-1)/2) for n odd and 0 otherwise
 282             //      = 2*(2/(n*pi))^2 * (-1)^((n-1)/2)
 283             if (n & 1) {
 284                 b = 2 * (piFactor * piFactor) * ((((n - 1) >> 1) & 1) ? -1 : 1);
 285             } else {
 286                 b = 0;
 287             }
 288             break;
 289         default:
 290             ASSERT_NOT_REACHED();
 291             b = 0;
 292             break;
 293         }
 294
 295         realP[n] = 0;
 296         imagP[n] = b;
 297     }
 298
 299     createBandLimitedTables(realP, imagP, halfSize);
 300 }
 301
 302 } // namespace blink
 303
 304 #endif // ENABLE(WEB_AUDIO)