<?php
/**
 * SeekQuarry/Yioop --
 * Open Source Pure PHP Search Engine, Crawler, and Indexer
 *
 * Copyright (C) 2009 - 2026  Chris Pollett chris@pollett.org
 *
 * LICENSE:
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 *
 * END LICENSE
 *
 * @author Chris Pollett chris@pollett.org
 * @license https://www.gnu.org/licenses/ GPL3
 * @link https://www.seekquarry.com/
 * @copyright 2009 - 2026
 * @filesource
 */
namespace seekquarry\yioop\library\av_processing;

/**
 * Mdct turns overlapping stretches of sound into tones and back, in a way that
 * loses nothing despite the stretches overlapping. Sound is cut into stretches
 * before being turned into tones, because what a sound is made of changes as it
 * goes on. Cutting it leaves joins, and a join heard on its own is a click, so
 * the stretches are faded in and out and made to overlap. Fading and
 * overlapping would normally mean storing twice as many numbers as there were
 * samples. This transform avoids that: it gives back half as many tones as it
 * was given samples, and the halves that appear to be missing come back when
 * two neighboring stretches are added together. What one stretch gets wrong at
 * its edge, the next gets wrong in the opposite direction, and the two cancel
 * exactly. That cancelling only works if the fade satisfies one condition:
 * squared, a fade and its own reverse must add to one everywhere. Both fades
 * offered here satisfy it. Every part of this work needs the transform. Opus
 * uses it to turn stored tones back into sound, and both of the formats being
 * written out use it again in the other direction.
 */
class Mdct
{
    /**
     * plans stores how the work for each length is split, worked out once per
     * length and kept for reuse
     *
     * @var array
     */
    public static $plans = [];
    /**
     * size stores how many tones this gives back, which is half the samples it
     * is given.
     * @var int
     */
    public $size;
    /**
     * $fourier stores the transform that splits a run of numbers into
     * the tones it is made of, which this longer transform is built
     * on top of.
     * @var object
     */
    public $fourier;
    /**
     * turn_across stores the sideways part of the turning factors used before
     * and after the splitting.
     * @var array
     */
    public $turn_across;
    /**
     * turn_up stores the upward part of those turning factors.
     * @var array
     */
    public $turn_up;
    /**
     * forSize works out how to handle stretches giving a given number of tones,
     * reusing the answer where one has been worked out before
     *
     * @param int $size how many tones to give back
     * @return object something that can handle stretches of that size
     */
    public static function forSize($size)
    {
        if (isset(self::$plans[$size])) {
            return self::$plans[$size];
        }
        $plan = new self($size);
        self::$plans[$size] = $plan;
        return $plan;
    }
    /**
     * __construct works out the turning factors a given size needs
     *
     * @param int $size how many tones to give back
     */
    public function __construct($size)
    {
        if ($size < 2 || $size % 2 != 0) {
            throw new \Exception("The number of tones must be even");
        }
        $this->size = $size;
        $half = intdiv($size, 2);
        $this->fourier = Fft::forSize($half);
        $this->turn_across = [];
        $this->turn_up = [];
        /* The turning is split evenly between before and after the
           splitting, so the two together come to a quarter step. */
        for ($i = 0; $i < $half; $i++) {
            $angle = -M_PI * ($i + 0.125) / $size;
            $this->turn_across[$i] = cos($angle);
            $this->turn_up[$i] = sin($angle);
        }
    }
    /**
     * forward turns a stretch of sound into the tones that make it up. The
     * stretch is twice as long as the run of tones that comes back.
     *
     * @param array $sound the stretch of sound, already faded
     * @return array the tones it is made of
     */
    public function forward($sound)
    {
        $size = $this->size;
        $quarter = intdiv($size, 2);
        $folded = array_fill(0, $size, 0.0);
        $three_halves = intdiv(3 * $size, 2);
        /* The stretch is folded onto itself first. The transform that
           does the real work takes as many numbers as it gives back,
           and this is what turns twice as many into that many without
           losing anything. */
        for ($i = 0; $i < $quarter; $i++) {
            $folded[$i] = -$sound[$three_halves - 1 - $i] -
                $sound[$three_halves + $i];
            $folded[$quarter + $i] = $sound[$i] - $sound[$size - 1 - $i];
        }
        return $this->tones($folded);
    }
    /**
     * inverse turns tones back into a stretch of sound, twice as long as the
     * run of tones given to its neighbor
     *
     * @param array $tones the tones to turn back
     * @return array the stretch of sound, ready to be faded and added
     */
    public function inverse($tones)
    {
        $size = $this->size;
        $quarter = intdiv($size, 2);
        $spread = $this->tones($tones);
        $scale = 2.0 / $size;
        $sound = array_fill(0, 2 * $size, 0.0);
        $three_halves = intdiv(3 * $size, 2);
        for ($i = 0; $i < $quarter; $i++) {
            $front = $spread[$quarter + $i] * $scale;
            $back = $spread[$i] * $scale;
            $sound[$i] = $front;
            $sound[$size - 1 - $i] = -$front;
            $sound[$three_halves - 1 - $i] = -$back;
            $sound[$three_halves + $i] = -$back;
        }
        return $sound;
    }
    /**
     * tones runs the transform both directions are built on. It turns a
     * run of numbers into the strengths of the tones that make it up,
     * and running it again gives the numbers back apart from a scaling,
     * which is why one method serves both ways.
     *
     * @param array $given The run of numbers to transform.
     * @return array The strength of each tone.
     */
    public function tones($given)
    {
        $size = $this->size;
        $half = intdiv($size, 2);
        $across = array_fill(0, $half, 0.0);
        $up = array_fill(0, $half, 0.0);
        /* The run is read as pairs, one from the front and one from
           the back, so that a run of a given length can be handled by
           splitting a run of half that length. */
        for ($i = 0; $i < $half; $i++) {
            $first = $given[2 * $i];
            $second = $given[$size - 1 - 2 * $i];
            $across[$i] = $first * $this->turn_across[$i] -
                $second * $this->turn_up[$i];
            $up[$i] = $first * $this->turn_up[$i] +
                $second * $this->turn_across[$i];
        }
        $this->fourier->run($across, $up);
        $out = array_fill(0, $size, 0.0);
        for ($i = 0; $i < $half; $i++) {
            $turned_across = $across[$i] * $this->turn_across[$i] -
                $up[$i] * $this->turn_up[$i];
            $turned_up = $across[$i] * $this->turn_up[$i] +
                $up[$i] * $this->turn_across[$i];
            $out[2 * $i] = $turned_across;
            $out[$size - 1 - 2 * $i] = -$turned_up;
        }
        return $out;
    }
    /**
     * shapedFade a fade that falls away faster at its ends, which lets less of
     * one stretch leak into the next at the cost of a wider middle
     *
     * @param int $length how long the fade is
     * @param float $sharpness how fast it falls away at the ends
     * @return array the fade, one value for each sample
     */
    public static function shapedFade($length, $sharpness = 4.0)
    {
        $half = intdiv($length, 2);
        $running = [];
        $total = 0.0;
        for ($i = 0; $i <= $half; $i++) {
            $offset = 2.0 * $i / $half - 1.0;
            $total += self::bessel(M_PI * $sharpness *
                sqrt(max(0.0, 1.0 - $offset * $offset)));
            $running[$i] = $total;
        }
        $fade = [];
        for ($i = 0; $i < $half; $i++) {
            $value = sqrt($running[$i] / $total);
            $fade[$i] = $value;
            $fade[$length - 1 - $i] = $value;
        }
        return $fade;
    }
    /**
     * bessel the shape the fade above is built from, worked out by adding up
     * terms until they stop making a difference
     *
     * @param float $at where to work the shape out
     * @return float the shape's value there
     */
    public static function bessel($at)
    {
        $total = 1.0;
        $term = 1.0;
        for ($i = 1; $i < 50; $i++) {
            $term *= ($at / (2.0 * $i)) * ($at / (2.0 * $i));
            $total += $term;
            if ($term < $total * 1e-16) {
                break;
            }
        }
        return $total;
    }
}
