/*
 * 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 Computer, 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 "PeriodicWave.h"

#include <algorithm>
#include <cmath>
#include <numbers>

#include "mozilla/FFTBlock.h"

const unsigned MinPeriodicWaveSize = 4096;  // This must be a power of two.
const unsigned MaxPeriodicWaveSize = 8192;  // This must be a power of two.
const float CentsPerRange = 1200 / 3;       // 1/3 Octave.

using namespace mozilla;
using mozilla::dom::OscillatorType;

namespace WebCore {

already_AddRefed<PeriodicWave> PeriodicWave::create(float sampleRate,
                                                    const float* real,
                                                    const float* imag,
                                                    size_t numberOfComponents,
                                                    bool disableNormalization) {
  bool isGood = real && imag && numberOfComponents > 0;
  MOZ_ASSERT(isGood);
  if (isGood) {
    RefPtr<PeriodicWave> periodicWave =
        new PeriodicWave(sampleRate, numberOfComponents, disableNormalization);

    // Limit the number of components used to those for frequencies below the
    // Nyquist of the fixed length inverse FFT.
    size_t halfSize = periodicWave->m_periodicWaveSize / 2;
    numberOfComponents = std::min(numberOfComponents, halfSize);
    periodicWave->m_numberOfComponents = numberOfComponents;
    periodicWave->m_realComponents =
        MakeUnique<AudioFloatArray>(numberOfComponents);
    periodicWave->m_imagComponents =
        MakeUnique<AudioFloatArray>(numberOfComponents);
    memcpy(periodicWave->m_realComponents->Elements(), real,
           numberOfComponents * sizeof(float));
    memcpy(periodicWave->m_imagComponents->Elements(), imag,
           numberOfComponents * sizeof(float));

    return periodicWave.forget();
  }
  return nullptr;
}

already_AddRefed<PeriodicWave> PeriodicWave::createSine(float sampleRate) {
  RefPtr<PeriodicWave> periodicWave =
      new PeriodicWave(sampleRate, MinPeriodicWaveSize, false);
  periodicWave->generateBasicWaveform(OscillatorType::Sine);
  return periodicWave.forget();
}

already_AddRefed<PeriodicWave> PeriodicWave::createSquare(float sampleRate) {
  RefPtr<PeriodicWave> periodicWave =
      new PeriodicWave(sampleRate, MinPeriodicWaveSize, false);
  periodicWave->generateBasicWaveform(OscillatorType::Square);
  return periodicWave.forget();
}

already_AddRefed<PeriodicWave> PeriodicWave::createSawtooth(float sampleRate) {
  RefPtr<PeriodicWave> periodicWave =
      new PeriodicWave(sampleRate, MinPeriodicWaveSize, false);
  periodicWave->generateBasicWaveform(OscillatorType::Sawtooth);
  return periodicWave.forget();
}

already_AddRefed<PeriodicWave> PeriodicWave::createTriangle(float sampleRate) {
  RefPtr<PeriodicWave> periodicWave =
      new PeriodicWave(sampleRate, MinPeriodicWaveSize, false);
  periodicWave->generateBasicWaveform(OscillatorType::Triangle);
  return periodicWave.forget();
}

PeriodicWave::PeriodicWave(float sampleRate, size_t numberOfComponents,
                           bool disableNormalization)
    : m_sampleRate(sampleRate),
      m_centsPerRange(CentsPerRange),
      m_maxPartialsInBandLimitedTable(0),
      m_normalizationScale(1.0f),
      m_disableNormalization(disableNormalization) {
  float nyquist = 0.5 * m_sampleRate;

  if (numberOfComponents <= MinPeriodicWaveSize) {
    m_periodicWaveSize = MinPeriodicWaveSize;
  } else {
    unsigned npow2 = fdlibm_exp2f(floorf(
        fdlibm_logf(numberOfComponents - 1.0) / fdlibm_logf(2.0f) + 1.0f));
    m_periodicWaveSize = std::min(MaxPeriodicWaveSize, npow2);
  }

  m_numberOfRanges =
      (unsigned)(3.0f * fdlibm_logf(m_periodicWaveSize) / fdlibm_logf(2.0f));
  m_bandLimitedTables.SetLength(m_numberOfRanges);
  m_lowestFundamentalFrequency = nyquist / maxNumberOfPartials();
  m_rateScale = m_periodicWaveSize / m_sampleRate;
}

size_t PeriodicWave::sizeOfIncludingThis(
    mozilla::MallocSizeOf aMallocSizeOf) const {
  size_t amount = aMallocSizeOf(this);

  amount += m_bandLimitedTables.ShallowSizeOfExcludingThis(aMallocSizeOf);
  for (size_t i = 0; i < m_bandLimitedTables.Length(); i++) {
    if (m_bandLimitedTables[i]) {
      amount +=
          m_bandLimitedTables[i]->ShallowSizeOfIncludingThis(aMallocSizeOf);
    }
  }

  return amount;
}

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);

  // We only need to rebuild to the tables if the new fundamental
  // frequency is low enough to allow for more partials below the
  // Nyquist frequency.
  unsigned numberOfPartials = numberOfPartialsForRange(0);
  float nyquist = 0.5 * m_sampleRate;
  if (fundamentalFrequency != 0.0) {
    numberOfPartials =
        std::min(numberOfPartials, (unsigned)(nyquist / fundamentalFrequency));
  }
  if (numberOfPartials > m_maxPartialsInBandLimitedTable) {
    for (unsigned rangeIndex = 0; rangeIndex < m_numberOfRanges; ++rangeIndex) {
      m_bandLimitedTables[rangeIndex] = nullptr;
    }

    // We need to create the first table to determine the normalization
    // constant.
    createBandLimitedTables(fundamentalFrequency, 0);
    m_maxPartialsInBandLimitedTable = numberOfPartials;
  }

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

  // Add one to round-up to the next range just in time to truncate
  // partials before aliasing occurs.
  float pitchRange = 1 + centsAboveLowestFrequency / m_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;

  if (!m_bandLimitedTables[rangeIndex1].get())
    createBandLimitedTables(fundamentalFrequency, rangeIndex1);

  if (!m_bandLimitedTables[rangeIndex2].get())
    createBandLimitedTables(fundamentalFrequency, rangeIndex2);

  lowerWaveData = m_bandLimitedTables[rangeIndex2]->Elements();
  higherWaveData = m_bandLimitedTables[rangeIndex1]->Elements();

  // Ranges from 0 -> 1 to interpolate between lower -> higher.
  tableInterpolationFactor = rangeIndex2 - pitchRange;
}

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

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

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

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

  return numberOfPartials;
}

// Convert into time-domain wave buffers.
// 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(float fundamentalFrequency,
                                           unsigned rangeIndex) {
  unsigned fftSize = m_periodicWaveSize;
  unsigned i;

  const float* realData = m_realComponents->Elements();
  const float* imagData = m_imagComponents->Elements();

  // This FFTBlock is used to cull partials (represented by frequency bins).
  FFTBlock frame(fftSize);

  // Find the starting bin where we should start culling the aliasing
  // partials for this pitch range.  We need to clear out the highest
  // frequencies to band-limit the waveform.
  unsigned numberOfPartials = numberOfPartialsForRange(rangeIndex);
  // Also limit to the number of components that are provided.
  numberOfPartials = std::min(numberOfPartials, m_numberOfComponents - 1);

  // Limit number of partials to those below Nyquist frequency
  float nyquist = 0.5 * m_sampleRate;
  if (fundamentalFrequency != 0.0) {
    numberOfPartials =
        std::min(numberOfPartials, (unsigned)(nyquist / fundamentalFrequency));
  }

  // Copy from loaded frequency data and generate complex conjugate
  // because of the way the inverse FFT is defined.
  // The coefficients of higher partials remain zero, as initialized in
  // the FFTBlock constructor.
  for (i = 0; i < numberOfPartials + 1; ++i) {
    frame.RealData(i) = realData[i];
    frame.ImagData(i) = -imagData[i];
  }

  // Clear any DC-offset.
  frame.RealData(0) = 0;
  // Clear value which has no effect.
  frame.ImagData(0) = 0;

  // Create the band-limited table.
  m_bandLimitedTables[rangeIndex] =
      MakeUnique<AlignedAudioFloatArray>(m_periodicWaveSize);

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

  // For the first range (which has the highest power), calculate
  // its peak value then compute normalization scale.
  if (m_disableNormalization) {
    // See Bug 1424906, results need to be scaled by 0.5 even
    // when normalization is disabled.
    m_normalizationScale = 0.5;
  } else if (!rangeIndex) {
    float maxValue;
    maxValue = AudioBufferPeakValue(data, m_periodicWaveSize);

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

  // Apply normalization scale.
  AudioBufferInPlaceScale(data, m_normalizationScale, m_periodicWaveSize);
}

void PeriodicWave::generateBasicWaveform(OscillatorType shape) {
  const float piFloat = std::numbers::pi_v<float>;
  unsigned fftSize = periodicWaveSize();
  unsigned halfSize = fftSize / 2;

  m_numberOfComponents = halfSize;
  m_realComponents = MakeUnique<AudioFloatArray>(halfSize);
  m_imagComponents = MakeUnique<AudioFloatArray>(halfSize);
  float* realP = m_realComponents->Elements();
  float* imagP = m_imagComponents->Elements();

  // Clear DC and imag value which is ignored.
  realP[0] = 0;
  imagP[0] = 0;

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

    // Fourier coefficients according to standard definition.
    float a;  // Coefficient for cos().
    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 OscillatorType::Sine:
        // Standard sine wave function.
        a = 0;
        b = (n == 1) ? 1 : 0;
        break;
      case OscillatorType::Square:
        // Square-shaped waveform with the first half its maximum value
        // and the second half its minimum value.
        a = 0;
        b = invOmega * ((n & 1) ? 2 : 0);
        break;
      case OscillatorType::Sawtooth:
        // Sawtooth-shaped waveform with the first half ramping from
        // zero to maximum and the second half from minimum to zero.
        a = 0;
        b = -invOmega * fdlibm_cos(0.5 * omega);
        break;
      case OscillatorType::Triangle:
        // Triangle-shaped waveform going from its maximum value to
        // its minimum value then back to the maximum value.
        a = 0;
        if (n & 1) {
          b = 2 * (2 / (n * piFloat) * 2 / (n * piFloat)) *
              ((((n - 1) >> 1) & 1) ? -1 : 1);
        } else {
          b = 0;
        }
        break;
      default:
        MOZ_ASSERT_UNREACHABLE("invalid oscillator type");
        a = 0;
        b = 0;
        break;
    }

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

}  // namespace WebCore