Home | History | Annotate | Download | only in webaudio
      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