michael@0: /*
michael@0:  * Copyright (C) 2012 Google Inc. All rights reserved.
michael@0:  *
michael@0:  * Redistribution and use in source and binary forms, with or without
michael@0:  * modification, are permitted provided that the following conditions
michael@0:  * are met:
michael@0:  *
michael@0:  * 1.  Redistributions of source code must retain the above copyright
michael@0:  *     notice, this list of conditions and the following disclaimer.
michael@0:  * 2.  Redistributions in binary form must reproduce the above copyright
michael@0:  *     notice, this list of conditions and the following disclaimer in the
michael@0:  *     documentation and/or other materials provided with the distribution.
michael@0:  * 3.  Neither the name of Apple Computer, Inc. ("Apple") nor the names of
michael@0:  *     its contributors may be used to endorse or promote products derived
michael@0:  *     from this software without specific prior written permission.
michael@0:  *
michael@0:  * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
michael@0:  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
michael@0:  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
michael@0:  * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
michael@0:  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
michael@0:  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
michael@0:  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
michael@0:  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
michael@0:  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
michael@0:  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0:  */
michael@0: 
michael@0: #include "PeriodicWave.h"
michael@0: #include <algorithm>
michael@0: #include <cmath>
michael@0: #include "mozilla/FFTBlock.h"
michael@0: 
michael@0: const unsigned PeriodicWaveSize = 4096; // This must be a power of two.
michael@0: const unsigned NumberOfRanges = 36; // There should be 3 * log2(PeriodicWaveSize) 1/3 octave ranges.
michael@0: const float CentsPerRange = 1200 / 3; // 1/3 Octave.
michael@0: 
michael@0: using namespace mozilla;
michael@0: using mozilla::dom::OscillatorType;
michael@0: 
michael@0: namespace WebCore {
michael@0: 
michael@0: PeriodicWave* PeriodicWave::create(float sampleRate,
michael@0:                                    const float* real,
michael@0:                                    const float* imag,
michael@0:                                    size_t numberOfComponents)
michael@0: {
michael@0:     bool isGood = real && imag && numberOfComponents > 0 &&
michael@0:          numberOfComponents <= PeriodicWaveSize;
michael@0:     MOZ_ASSERT(isGood);
michael@0:     if (isGood) {
michael@0:         PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
michael@0:         periodicWave->createBandLimitedTables(real, imag, numberOfComponents);
michael@0:         return periodicWave;
michael@0:     }
michael@0:     return 0;
michael@0: }
michael@0: 
michael@0: PeriodicWave* PeriodicWave::createSine(float sampleRate)
michael@0: {
michael@0:       PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
michael@0:           periodicWave->generateBasicWaveform(OscillatorType::Sine);
michael@0:               return periodicWave;
michael@0: }
michael@0: 
michael@0: PeriodicWave* PeriodicWave::createSquare(float sampleRate)
michael@0: {
michael@0:       PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
michael@0:           periodicWave->generateBasicWaveform(OscillatorType::Square);
michael@0:               return periodicWave;
michael@0: }
michael@0: 
michael@0: PeriodicWave* PeriodicWave::createSawtooth(float sampleRate)
michael@0: {
michael@0:       PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
michael@0:           periodicWave->generateBasicWaveform(OscillatorType::Sawtooth);
michael@0:               return periodicWave;
michael@0: }
michael@0: 
michael@0: PeriodicWave* PeriodicWave::createTriangle(float sampleRate)
michael@0: {
michael@0:       PeriodicWave* periodicWave = new PeriodicWave(sampleRate);
michael@0:           periodicWave->generateBasicWaveform(OscillatorType::Triangle);
michael@0:               return periodicWave;
michael@0: }
michael@0: 
michael@0: PeriodicWave::PeriodicWave(float sampleRate)
michael@0:     : m_sampleRate(sampleRate)
michael@0:     , m_periodicWaveSize(PeriodicWaveSize)
michael@0:     , m_numberOfRanges(NumberOfRanges)
michael@0:     , m_centsPerRange(CentsPerRange)
michael@0: {
michael@0:     float nyquist = 0.5 * m_sampleRate;
michael@0:     m_lowestFundamentalFrequency = nyquist / maxNumberOfPartials();
michael@0:     m_rateScale = m_periodicWaveSize / m_sampleRate;
michael@0: }
michael@0: 
michael@0: size_t PeriodicWave::sizeOfIncludingThis(mozilla::MallocSizeOf aMallocSizeOf) const
michael@0: {
michael@0:     size_t amount = aMallocSizeOf(this);
michael@0: 
michael@0:     amount += m_bandLimitedTables.SizeOfExcludingThis(aMallocSizeOf);
michael@0:     for (size_t i = 0; i < m_bandLimitedTables.Length(); i++) {
michael@0:         if (m_bandLimitedTables[i]) {
michael@0:             amount += m_bandLimitedTables[i]->SizeOfIncludingThis(aMallocSizeOf);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     return amount;
michael@0: }
michael@0: 
michael@0: void PeriodicWave::waveDataForFundamentalFrequency(float fundamentalFrequency, float* &lowerWaveData, float* &higherWaveData, float& tableInterpolationFactor)
michael@0: {
michael@0:     // Negative frequencies are allowed, in which case we alias
michael@0:     // to the positive frequency.
michael@0:     fundamentalFrequency = fabsf(fundamentalFrequency);
michael@0: 
michael@0:     // Calculate the pitch range.
michael@0:     float ratio = fundamentalFrequency > 0 ? fundamentalFrequency / m_lowestFundamentalFrequency : 0.5;
michael@0:     float centsAboveLowestFrequency = logf(ratio)/logf(2.0f) * 1200;
michael@0: 
michael@0:     // Add one to round-up to the next range just in time to truncate
michael@0:     // partials before aliasing occurs.
michael@0:     float pitchRange = 1 + centsAboveLowestFrequency / m_centsPerRange;
michael@0: 
michael@0:     pitchRange = std::max(pitchRange, 0.0f);
michael@0:     pitchRange = std::min(pitchRange, static_cast<float>(m_numberOfRanges - 1));
michael@0: 
michael@0:     // The words "lower" and "higher" refer to the table data having
michael@0:     // the lower and higher numbers of partials. It's a little confusing
michael@0:     // since the range index gets larger the more partials we cull out.
michael@0:     // So the lower table data will have a larger range index.
michael@0:     unsigned rangeIndex1 = static_cast<unsigned>(pitchRange);
michael@0:     unsigned rangeIndex2 = rangeIndex1 < m_numberOfRanges - 1 ? rangeIndex1 + 1 : rangeIndex1;
michael@0: 
michael@0:     lowerWaveData = m_bandLimitedTables[rangeIndex2]->Elements();
michael@0:     higherWaveData = m_bandLimitedTables[rangeIndex1]->Elements();
michael@0: 
michael@0:     // Ranges from 0 -> 1 to interpolate between lower -> higher.
michael@0:     tableInterpolationFactor = pitchRange - rangeIndex1;
michael@0: }
michael@0: 
michael@0: unsigned PeriodicWave::maxNumberOfPartials() const
michael@0: {
michael@0:     return m_periodicWaveSize / 2;
michael@0: }
michael@0: 
michael@0: unsigned PeriodicWave::numberOfPartialsForRange(unsigned rangeIndex) const
michael@0: {
michael@0:     // Number of cents below nyquist where we cull partials.
michael@0:     float centsToCull = rangeIndex * m_centsPerRange;
michael@0: 
michael@0:     // A value from 0 -> 1 representing what fraction of the partials to keep.
michael@0:     float cullingScale = pow(2, -centsToCull / 1200);
michael@0: 
michael@0:     // The very top range will have all the partials culled.
michael@0:     unsigned numberOfPartials = cullingScale * maxNumberOfPartials();
michael@0: 
michael@0:     return numberOfPartials;
michael@0: }
michael@0: 
michael@0: // Convert into time-domain wave buffers.
michael@0: // One table is created for each range for non-aliasing playback
michael@0: // at different playback rates. Thus, higher ranges have more
michael@0: // high-frequency partials culled out.
michael@0: void PeriodicWave::createBandLimitedTables(const float* realData, const float* imagData, unsigned numberOfComponents)
michael@0: {
michael@0:     float normalizationScale = 1;
michael@0: 
michael@0:     unsigned fftSize = m_periodicWaveSize;
michael@0:     unsigned halfSize = fftSize / 2 + 1;
michael@0:     unsigned i;
michael@0: 
michael@0:     numberOfComponents = std::min(numberOfComponents, halfSize);
michael@0: 
michael@0:     m_bandLimitedTables.SetCapacity(m_numberOfRanges);
michael@0: 
michael@0:     for (unsigned rangeIndex = 0; rangeIndex < m_numberOfRanges; ++rangeIndex) {
michael@0:         // This FFTBlock is used to cull partials (represented by frequency bins).
michael@0:         FFTBlock frame(fftSize);
michael@0:         nsAutoArrayPtr<float> realP(new float[halfSize]);
michael@0:         nsAutoArrayPtr<float> imagP(new float[halfSize]);
michael@0: 
michael@0:         // Copy from loaded frequency data and scale.
michael@0:         float scale = fftSize;
michael@0:         AudioBufferCopyWithScale(realData, scale, realP, numberOfComponents);
michael@0:         AudioBufferCopyWithScale(imagData, scale, imagP, numberOfComponents);
michael@0: 
michael@0:         // If fewer components were provided than 1/2 FFT size,
michael@0:         // then clear the remaining bins.
michael@0:         for (i = numberOfComponents; i < halfSize; ++i) {
michael@0:             realP[i] = 0;
michael@0:             imagP[i] = 0;
michael@0:         }
michael@0: 
michael@0:         // Generate complex conjugate because of the way the
michael@0:         // inverse FFT is defined.
michael@0:         float minusOne = -1;
michael@0:         AudioBufferInPlaceScale(imagP, minusOne, halfSize);
michael@0: 
michael@0:         // Find the starting bin where we should start culling.
michael@0:         // We need to clear out the highest frequencies to band-limit
michael@0:         // the waveform.
michael@0:         unsigned numberOfPartials = numberOfPartialsForRange(rangeIndex);
michael@0: 
michael@0:         // Cull the aliasing partials for this pitch range.
michael@0:         for (i = numberOfPartials + 1; i < halfSize; ++i) {
michael@0:             realP[i] = 0;
michael@0:             imagP[i] = 0;
michael@0:         }
michael@0:         // Clear nyquist if necessary.
michael@0:         if (numberOfPartials < halfSize)
michael@0:             realP[halfSize-1] = 0;
michael@0: 
michael@0:         // Clear any DC-offset.
michael@0:         realP[0] = 0;
michael@0: 
michael@0:         // Clear values which have no effect.
michael@0:         imagP[0] = 0;
michael@0:         imagP[halfSize-1] = 0;
michael@0: 
michael@0:         // Create the band-limited table.
michael@0:         AudioFloatArray* table = new AudioFloatArray(m_periodicWaveSize);
michael@0:         m_bandLimitedTables.AppendElement(table);
michael@0: 
michael@0:         // Apply an inverse FFT to generate the time-domain table data.
michael@0:         float* data = m_bandLimitedTables[rangeIndex]->Elements();
michael@0:         frame.PerformInverseFFT(realP, imagP, data);
michael@0: 
michael@0:         // For the first range (which has the highest power), calculate
michael@0:         // its peak value then compute normalization scale.
michael@0:         if (!rangeIndex) {
michael@0:             float maxValue;
michael@0:             maxValue = AudioBufferPeakValue(data, m_periodicWaveSize);
michael@0: 
michael@0:             if (maxValue)
michael@0:                 normalizationScale = 1.0f / maxValue;
michael@0:         }
michael@0: 
michael@0:         // Apply normalization scale.
michael@0:         AudioBufferInPlaceScale(data, normalizationScale, m_periodicWaveSize);
michael@0:     }
michael@0: }
michael@0: 
michael@0: void PeriodicWave::generateBasicWaveform(OscillatorType shape)
michael@0: {
michael@0:     const float piFloat = M_PI;
michael@0:     unsigned fftSize = periodicWaveSize();
michael@0:     unsigned halfSize = fftSize / 2 + 1;
michael@0: 
michael@0:     AudioFloatArray real(halfSize);
michael@0:     AudioFloatArray imag(halfSize);
michael@0:     float* realP = real.Elements();
michael@0:     float* imagP = imag.Elements();
michael@0: 
michael@0:     // Clear DC and Nyquist.
michael@0:     realP[0] = 0;
michael@0:     imagP[0] = 0;
michael@0:     realP[halfSize-1] = 0;
michael@0:     imagP[halfSize-1] = 0;
michael@0: 
michael@0:     for (unsigned n = 1; n < halfSize; ++n) {
michael@0:         float omega = 2 * piFloat * n;
michael@0:         float invOmega = 1 / omega;
michael@0: 
michael@0:         // Fourier coefficients according to standard definition.
michael@0:         float a; // Coefficient for cos().
michael@0:         float b; // Coefficient for sin().
michael@0: 
michael@0:         // Calculate Fourier coefficients depending on the shape.
michael@0:         // Note that the overall scaling (magnitude) of the waveforms
michael@0:         // is normalized in createBandLimitedTables().
michael@0:         switch (shape) {
michael@0:         case OscillatorType::Sine:
michael@0:             // Standard sine wave function.
michael@0:             a = 0;
michael@0:             b = (n == 1) ? 1 : 0;
michael@0:             break;
michael@0:         case OscillatorType::Square:
michael@0:             // Square-shaped waveform with the first half its maximum value
michael@0:             // and the second half its minimum value.
michael@0:             a = 0;
michael@0:             b = invOmega * ((n & 1) ? 2 : 0);
michael@0:             break;
michael@0:         case OscillatorType::Sawtooth:
michael@0:             // Sawtooth-shaped waveform with the first half ramping from
michael@0:             // zero to maximum and the second half from minimum to zero.
michael@0:             a = 0;
michael@0:             b = -invOmega * cos(0.5 * omega);
michael@0:             break;
michael@0:         case OscillatorType::Triangle:
michael@0:             // Triangle-shaped waveform going from its maximum value to
michael@0:             // its minimum value then back to the maximum value.
michael@0:             a = (4 - 4 * cos(0.5 * omega)) / (n * n * piFloat * piFloat);
michael@0:             b = 0;
michael@0:             break;
michael@0:         default:
michael@0:             NS_NOTREACHED("invalid oscillator type");
michael@0:             a = 0;
michael@0:             b = 0;
michael@0:             break;
michael@0:         }
michael@0: 
michael@0:         realP[n] = a;
michael@0:         imagP[n] = b;
michael@0:     }
michael@0: 
michael@0:     createBandLimitedTables(realP, imagP, halfSize);
michael@0: }
michael@0: 
michael@0: } // namespace WebCore