Hoorsimulatie.php

hoorsimulatie.php

<?php

define('TWO_PI', 2 * M_PI);
define('FOUR_PI', 4 * M_PI);
define('FSAM', 44100);
define('N_OLA', 2);
define('N_FFT', 2048);
define('N_FFT2', N_FFT / 2);
define('N_TAP', N_FFT2);
define('N_ZERO', N_FFT2);
define('N_SAM', N_FFT - N_ZERO);
define('N_SHIFT', N_SAM / N_OLA);

define('MISDAT', -11);
define('MIN_FREQ', 1);
define('MAX_FREQ', 13);

$FREQ = [125, 187, 250, 375, 500, 750, 1000, 1500, 2000, 3000, 4000, 6000, 8000, 12000, 16000];
$SMEARING = [0.4, 0.5, 0.9, 1.2, 1.5, 1.8, 2.0, 2.3, 3.1, 3.9];

$Sum = array_fill(0, N_FFT, 0.0);
$FFT = new FFT(N_FFT);
$Int = array_fill(0, N_FFT, 0);
$Sam = array_fill(0, N_FFT, 0);
$Sam2 = array_fill(0, 2 * N_FFT, 0);
$yFir = array_fill(0, N_FFT, 0.0);
$Fir = array_fill(0, N_FFT, 0.0);
$fFFT = array_fill(0, N_FFT, 0.0);
$xFFT = array_fill(0, N_FFT, 0.0);
$yFFT = array_fill(0, N_FFT, 0.0);
$aFFT = array_fill(0, N_FFT, 0.0);
$zFFT = array_fill(0, N_FFT, 0.0);
$eFFT = array_fill(0, N_FFT, 0.0);
$lFFT = array_fill(0, N_FFT, 0.0);
$W = array_fill(0, N_FFT, 0.0);
$wFFT = array_fill(0, N_FFT, 0.0);
$rPow2 = array_fill(0, N_FFT, 0.0);
$iPow2 = array_fill(0, N_FFT, 0);

function LowPassTimeSeries($Fc, &$Fir) {
    $M = N_TAP;
    $Sum = 0.0;

    for ($i = 0; $i < $M; $i++) {
        if ($i - $M / 2 == 0) {
            $Fir[$i] = $Fc * TWO_PI;
        } else {
            $Fir[$i] = sin($Fc * TWO_PI * ($i - $M / 2)) / ($i - $M / 2);  // sinc(x)
        }
        $Fir[$i] *= 0.42 - 0.5 * cos(TWO_PI * $i / $M) + 0.08 * cos(FOUR_PI * $i / $M);  // windowing
    }

    for ($i = 0; $i < $M; $i++) {
        $Sum += $Fir[$i];
    }

    for ($i = 0; $i < $M; $i++) {
        $Fir[$i] /= $Sum;
    }
}

function HearingLossFilter(&$yFir, $T, $Ear) {
    global $FREQ, $FFT;

    $rGain = -1;
    $Gain = array_fill(0, MAX_FREQ + 3, 0.0);
    $Fir = array_fill(0, MAX_FREQ + 3, array_fill(0, N_FFT, 0.0));
    $xFir = array_fill(0, N_FFT, 0.0);

    if ($T[3] == MISDAT) $T[3] = $T[5];
    if ($T[1] == MISDAT) $T[1] = $T[3];
    if ($T[13] == MISDAT) $T[13] = $T[11];
    $T[14] = $T[13];
    $T[15] = $T[13];

    for ($iFreq = MIN_FREQ + 1; $iFreq <= MAX_FREQ + 1; $iFreq++) {
        if ($iFreq % 2 == 1) continue;
        if ($T[$iFreq] == MISDAT) {
            $T[$iFreq] = round(0.5 * ($T[$iFreq - 1] + $T[$iFreq + 1]));
        }
    }

    for ($iFreq = MIN_FREQ; $iFreq <= MAX_FREQ + 2; $iFreq++) {
        $Gain[$iFreq] = pow(10, -0.05 * $T[$iFreq]);
        LowPassTimeSeries($FREQ[$iFreq] * pow(2, 0.25) / FSAM, $Fir[$iFreq]);
    }

    $Gain[MAX_FREQ + 3] = $Gain[MAX_FREQ + 2];
    LowPassTimeSeries(0.5, $Fir[MAX_FREQ + 3]);

    for ($i = 0; $i < N_FFT; $i++) {
        $xFir[$i] = 0;
        $yFir[$i] = 0;
    }

    for ($iFreq = MAX_FREQ + 3; $iFreq >= MIN_FREQ + 1; $iFreq--) {
        for ($i = 0; $i < N_TAP; $i++) {
            $xFir[$i] += ($Fir[$iFreq][$i] - $Fir[$iFreq - 1][$i]) * $Gain[$iFreq];
        }
    }

    for ($i = 0; $i < N_TAP; $i++) {
        $xFir[$i] += $Fir[1][$i] * $Gain[1];
    }

    $FFT->DoFFT($xFir, $yFir);
}

function SmearingFilter(&$wFFT, $B) {
    global $FFT, $rPow2, $iPow2, $W;

    for ($i = 0; $i < N_FFT2; $i++) {
        $rPow2[$i] = pow(2, $i * log(N_FFT2) / log(2) / N_FFT2);
        $iPow2[$i] = round($rPow2[$i]);
    }

    if ($B < 3.0) $B = 3.0;

    $Sum = 0;
    for ($i = 0; $i < N_TAP; $i++) {
        $W[$i] = exp(-M_PI * pow(($i - N_TAP / 2) / $B, 2));
        $Sum += $W[$i];
    }

    for ($i = 0; $i < N_TAP; $i++) {
        $W[$i] /= $Sum;
    }

    for ($i = N_TAP; $i < N_FFT; $i++) {
        $W[$i] = 0;
    }

    $FFT->DoFFT($W, $wFFT);
}

class FFT {
    private $nFFT;

    public function __construct($nFFT) {
        $this->nFFT = $nFFT;
    }

    public function DoFFT($input, &$output) {
        // Placeholder for actual FFT implementation
    }

    public function DoIFFT($input, &$output) {
        // Placeholder for actual IFFT implementation
    }

    public function DoFilter($filter, &$data) {
        // Placeholder for filtering logic
    }

    public function DoHamming(&$data, $nSam) {
        // Apply Hamming window to the data
        for ($i = 0; $i < $nSam; $i++) {
            $data[$i] *= 0.54 - 0.46 * cos(TWO_PI * $i / ($nSam - 1));
        }
    }
}

// Initialize FFT instance
$FFT = new FFT(N_FFT);

// Initialization
$FFT = new FFT(N_FFT);

// Finalization
unset($FFT);