/*
 * Copyright (C) 2012 Google Inc. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1.  Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 * 2.  Redistributions in binary form must reproduce the above copyright
 *     notice, this list of conditions and the following disclaimer in the
 *     documentation and/or other materials provided with the distribution.
 * 3.  Neither the name of Apple Inc. ("Apple") nor the names of
 *     its contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "config.h"

#if ENABLE(WEB_AUDIO)

#include "PeriodicWave.h"

#include "BaseAudioContext.h"
#include "FFTFrame.h"
#include "VectorMath.h"
#include <algorithm>

// The number of bands per octave. Each octave will have this many entries in the wave tables.
constexpr unsigned NumberOfOctaveBands = 3;

// The max length of a periodic wave. This must be a power of two greater than
// or equal to 2048 and must be supported by the FFT routines.
constexpr unsigned MaxPeriodicWaveSize = 16384;

constexpr float CentsPerRange = 1200 / NumberOfOctaveBands;

namespace WebCore {

Ref<PeriodicWave> PeriodicWave::create(float sampleRate, Float32Array& real, Float32Array& imaginary)
{
    ASSERT(real.length() == imaginary.length());

    auto waveTable = adoptRef(*new PeriodicWave(sampleRate));
    waveTable->createBandLimitedTables(real.data(), imaginary.data(), real.length());
    return waveTable;
}

ExceptionOr<Ref<PeriodicWave>> PeriodicWave::create(BaseAudioContext& context, PeriodicWaveOptions&& options)
{
    Vector<float> real;
    Vector<float> imag;
    
    if (options.real && options.imag) {
        if (options.real->size() != options.imag->size())
            return Exception { IndexSizeError, "real and imag have different lengths"_s };
        if (options.real->size() < 2)
            return Exception { IndexSizeError, "real's length cannot be less than 2"_s };
        if (options.imag->size() < 2)
            return Exception { IndexSizeError, "imag's length cannot be less than 2"_s };
        real = WTFMove(*options.real);
        imag = WTFMove(*options.imag);
    } else if (options.real) {
        if (options.real->size() < 2)
            return Exception { IndexSizeError, "real's length cannot be less than 2"_s };
        real = WTFMove(*options.real);
        imag.fill(0, real.size());
    } else if (options.imag) {
        if (options.imag->size() < 2)
            return Exception { IndexSizeError, "imag's length cannot be less than 2"_s };
        imag = WTFMove(*options.imag);
        real.fill(0, imag.size());
    } else {
        real.fill(0, 2);
        imag.fill(0, 2);
        imag[1] = 1;
    }

    real[0] = 0;
    imag[0] = 0;
    
    auto waveTable = adoptRef(*new PeriodicWave(context.sampleRate()));
    waveTable->createBandLimitedTables(real.data(), imag.data(), real.size(), options.disableNormalization ? ShouldDisableNormalization::Yes : ShouldDisableNormalization::No);
    return waveTable;
}

Ref<PeriodicWave> PeriodicWave::createSine(float sampleRate)
{
    Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
    waveTable->generateBasicWaveform(Type::Sine);
    return waveTable;
}

Ref<PeriodicWave> PeriodicWave::createSquare(float sampleRate)
{
    Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
    waveTable->generateBasicWaveform(Type::Square);
    return waveTable;
}

Ref<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate)
{
    Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
    waveTable->generateBasicWaveform(Type::Sawtooth);
    return waveTable;
}

Ref<PeriodicWave> PeriodicWave::createTriangle(float sampleRate)
{
    Ref<PeriodicWave> waveTable = adoptRef(*new PeriodicWave(sampleRate));
    waveTable->generateBasicWaveform(Type::Triangle);
    return waveTable;
}

PeriodicWave::PeriodicWave(float sampleRate)
    : m_sampleRate(sampleRate)
    , m_numberOfRanges(0.5 + NumberOfOctaveBands * log2f(periodicWaveSize()))
{
    float nyquist = 0.5 * m_sampleRate;
    m_lowestFundamentalFrequency = nyquist / maxNumberOfPartials();
    m_rateScale = periodicWaveSize() / m_sampleRate;
}

void PeriodicWave::waveDataForFundamentalFrequency(float fundamentalFrequency, float* &lowerWaveData, float* &higherWaveData, float& tableInterpolationFactor)
{
    // Negative frequencies are allowed, in which case we alias to the positive frequency.
    fundamentalFrequency = fabsf(fundamentalFrequency);

    // Calculate the pitch range.
    float ratio = fundamentalFrequency > 0 ? fundamentalFrequency / m_lowestFundamentalFrequency : 0.5;
    float centsAboveLowestFrequency = log2f(ratio) * 1200;

    // Add one to round-up to the next range just in time to truncate partials before aliasing occurs.
    float pitchRange = 1 + centsAboveLowestFrequency / CentsPerRange;

    pitchRange = std::max(pitchRange, 0.0f);
    pitchRange = std::min(pitchRange, static_cast<float>(m_numberOfRanges - 1));

    // The words "lower" and "higher" refer to the table data having the lower and higher numbers of partials.
    // It's a little confusing since the range index gets larger the more partials we cull out.
    // So the lower table data will have a larger range index.
    unsigned rangeIndex1 = static_cast<unsigned>(pitchRange);
    unsigned rangeIndex2 = rangeIndex1 < m_numberOfRanges - 1 ? rangeIndex1 + 1 : rangeIndex1;

    lowerWaveData = m_bandLimitedTables[rangeIndex2]->data();
    higherWaveData = m_bandLimitedTables[rangeIndex1]->data();
    
    // Ranges from 0 -> 1 to interpolate between lower -> higher.
    tableInterpolationFactor = pitchRange - rangeIndex1;
}

unsigned PeriodicWave::maxNumberOfPartials() const
{
    return periodicWaveSize() / 2;
}

unsigned PeriodicWave::numberOfPartialsForRange(unsigned rangeIndex) const
{
    // Number of cents below nyquist where we cull partials.
    float centsToCull = rangeIndex * CentsPerRange;

    // A value from 0 -> 1 representing what fraction of the partials to keep.
    float cullingScale = pow(2, -centsToCull / 1200);

    // The very top range will have all the partials culled.
    unsigned numberOfPartials = cullingScale * maxNumberOfPartials();

    return numberOfPartials;
}

// Convert into time-domain wave tables.
// One table is created for each range for non-aliasing playback at different playback rates.
// Thus, higher ranges have more high-frequency partials culled out.
void PeriodicWave::createBandLimitedTables(const float* realData, const float* imagData, unsigned numberOfComponents, ShouldDisableNormalization disableNormalization)
{
    float normalizationScale = 0.5;

    unsigned fftSize = periodicWaveSize();
    unsigned halfSize = fftSize / 2;
    
    numberOfComponents = std::min(numberOfComponents, halfSize);

    m_bandLimitedTables.reserveCapacity(m_numberOfRanges);

    for (unsigned rangeIndex = 0; rangeIndex < m_numberOfRanges; ++rangeIndex) {
        // This FFTFrame is used to cull partials (represented by frequency bins).
        FFTFrame frame(fftSize);
        auto& realP = frame.realData();
        auto& imagP = frame.imagData();

        RELEASE_ASSERT(realP.size() >= numberOfComponents);
        RELEASE_ASSERT(imagP.size() >= numberOfComponents);

        // Copy from loaded frequency data and scale.
        VectorMath::multiplyByScalar(realData, fftSize, realP.data(), numberOfComponents);
        VectorMath::multiplyByScalar(imagData, -static_cast<float>(fftSize), imagP.data(), numberOfComponents);

        // Find the starting bin where we should start culling.
        // We need to clear out the highest frequencies to band-limit the waveform.
        unsigned numberOfPartials = numberOfPartialsForRange(rangeIndex);

        // If fewer components were provided than 1/2 FFT size, then clear the
        // remaining bins. We also need to cull the aliasing partials for this
        // pitch range.
        unsigned clampedNumberOfComponents = std::min(numberOfComponents, numberOfPartials + 1);
        if (clampedNumberOfComponents < halfSize) {
            size_t numBytes = (halfSize - clampedNumberOfComponents) * sizeof(float);
            memset(&realP[clampedNumberOfComponents], 0, numBytes);
            memset(&imagP[clampedNumberOfComponents], 0, numBytes);
        }

        // Clear packed-nyquist and any DC-offset.
        realP[0] = 0;
        imagP[0] = 0;

        // Create the band-limited table.
        unsigned waveSize = periodicWaveSize();
        m_bandLimitedTables.append(makeUnique<AudioFloatArray>(waveSize));

        // Apply an inverse FFT to generate the time-domain table data.
        float* data = m_bandLimitedTables[rangeIndex]->data();
        frame.doInverseFFT(data);

        // For the first range (which has the highest power), calculate its peak value then compute normalization scale.
        if (disableNormalization == ShouldDisableNormalization::No) {
            if (!rangeIndex) {
                float maxValue = VectorMath::maximumMagnitude(data, fftSize);

                if (maxValue)
                    normalizationScale = 1.0f / maxValue;
            }
        }

        // Apply normalization scale.
        VectorMath::multiplyByScalar(data, normalizationScale, data, fftSize);
    }
}

void PeriodicWave::generateBasicWaveform(Type shape)
{
    unsigned fftSize = periodicWaveSize();
    unsigned halfSize = fftSize / 2;

    AudioFloatArray real(halfSize);
    AudioFloatArray imag(halfSize);
    float* realP = real.data();
    float* imagP = imag.data();

    // Clear DC and Nyquist.
    realP[0] = 0;
    imagP[0] = 0;

    for (unsigned n = 1; n < halfSize; ++n) {
        float piFactor = 2 / (n * piFloat);

        // All waveforms are odd functions with a positive slope at time 0. Hence
        // the coefficients for cos() are always 0.

        // Fourier coefficients according to standard definition:
        // b = 1/pi*integrate(f(x)*sin(n*x), x, -pi, pi)
        //   = 2/pi*integrate(f(x)*sin(n*x), x, 0, pi)
        // since f(x) is an odd function.

        float b; // Coefficient for sin().

        // Calculate Fourier coefficients depending on the shape.
        // Note that the overall scaling (magnitude) of the waveforms is normalized in createBandLimitedTables().
        switch (shape) {
        case Type::Sine:
            // Standard sine wave function.
            b = (n == 1) ? 1 : 0;
            break;
        case Type::Square:
            // Square-shaped waveform with the first half its maximum value and the
            // second half its minimum value.
            //
            // See http://mathworld.wolfram.com/FourierSeriesSquareWave.html
            //
            // b[n] = 2/n/pi*(1-(-1)^n)
            //      = 4/n/pi for n odd and 0 otherwise.
            //      = 2*(2/(n*pi)) for n odd
            b = (n & 1) ? 2 * piFactor : 0;
            break;
        case Type::Sawtooth:
            // Sawtooth-shaped waveform with the first half ramping from zero to
            // maximum and the second half from minimum to zero.
            //
            // b[n] = -2*(-1)^n/pi/n
            //      = (2/(n*pi))*(-1)^(n+1)
            b = piFactor * ((n & 1) ? 1 : -1);
            break;
        case Type::Triangle:
            // Triangle-shaped waveform going from 0 at time 0 to 1 at time pi/2 and
            // back to 0 at time pi.
            //
            // See http://mathworld.wolfram.com/FourierSeriesTriangleWave.html
            //
            // b[n] = 8*sin(pi*k/2)/(pi*k)^2
            //      = 8/pi^2/n^2*(-1)^((n-1)/2) for n odd and 0 otherwise
            //      = 2*(2/(n*pi))^2 * (-1)^((n-1)/2)
            if (n & 1)
                b = 2 * (piFactor * piFactor) * ((((n - 1) >> 1) & 1) ? -1 : 1);
            else
                b = 0;
            break;
        }

        realP[n] = 0;
        imagP[n] = b;
    }

    createBandLimitedTables(realP, imagP, halfSize);
}

unsigned PeriodicWave::periodicWaveSize() const
{
    // Choose an appropriate wave size for the given sample rate. This allows us
    // to use shorter FFTs when possible to limit the complexity. The breakpoints
    // here are somewhat arbitrary, but we want sample rates around 44.1 kHz or so
    // to have a size of 4096 to preserve backward compatibility.
    if (m_sampleRate <= 24000)
        return 2048;

    if (m_sampleRate <= 88200)
        return 4096;

    return MaxPeriodicWaveSize;
}

} // namespace WebCore

#endif // ENABLE(WEB_AUDIO)