#include "PhiScientificModel.h"

/*
 ===========================================================
    PhiScientificModel.cpp — VERSION SCIENTIFIQUEMENT COHÉRENTE
    
    MODULE: Analyse spectrale IMU + Audio/EMI
    
    OBJECTIF SCIENTIFIQUE:
    1. Quantifier énergie vibratoire mécanique (IMU)
    2. Classifier signaux acoustiques (Audio vs EMI)
    3. Détecter couplages vibro-acoustiques
 ===========================================================
*/

PhiScientificModel::PhiScientificModel()
: _fftSize(128),
  _fsIMU(90.0f),      // ✅ FIX: 90 Hz cohérent avec PhiIMU
  _fsAUD(16000.0f),
  _eIMU(0), _eAUD(0),
  _fIMU(0), _fAUD(0),
  _corr(0),
  _audioType(PhiAudio::TYPE_SILENCE),
  _spectralFlatness(0.0f),
  _previousAudioType(PhiAudio::TYPE_SILENCE),
  _fftIMU(nullptr),
  _fftAUD(nullptr)
{}

void PhiScientificModel::begin(size_t fftSize, float fsIMU, float fsAUD) {
    _fftSize = fftSize;
    _fsIMU = fsIMU;
    _fsAUD = fsAUD;

    _im_real = new float[_fftSize];
    _im_imag = new float[_fftSize];
    _au_real = new float[_fftSize];
    _au_imag = new float[_fftSize];
    _hann = new float[_fftSize];

    _fftIMU = new ArduinoFFT<float>(_im_real, _im_imag, _fftSize, _fsIMU);
    _fftAUD = new ArduinoFFT<float>(_au_real, _au_imag, _fftSize, _fsAUD);

    for (size_t i = 0; i < _fftSize; i++)
        _hann[i] = 0.5f * (1.0f - cosf(2.0f * PI * i / (_fftSize - 1)));

    Serial.println("✓ PhiScientificModel (Audio/EMI spectral classification)");
}

// ═══════════════════════════════════════════════════════════════════
// EXTRACTION FRÉQUENCE DOMINANTE
// ═══════════════════════════════════════════════════════════════════
// Retourne fréquence du bin FFT de magnitude maximale
// Usage: identification harmoniques fondamentales
// ═══════════════════════════════════════════════════════════════════

float PhiScientificModel::_dominantFreq(const float* mag, float fs) {
    float maxAmp = 0;
    int maxBin = 0;

    for (size_t i = 1; i < _fftSize/2; i++) {
        if (mag[i] > maxAmp) {
            maxAmp = mag[i];
            maxBin = i;
        }
    }

    return (fs * maxBin) / float(_fftSize);
}

// ═══════════════════════════════════════════════════════════════════
// EXTRACTION AMPLITUDE SPECTRALE DOMINANTE
// ═══════════════════════════════════════════════════════════════════
// Retourne magnitude maximale du spectre FFT
// 
// CLARIFICATION VOCABULAIRE:
// "Amplitude dominante" (pas "énergie") car il s'agit du maximum
// ponctuel, pas d'une intégration sur bande.
// 
// Énergie spectrale totale serait: E = Σ|X(f)|²
// Amplitude dominante est: A_max = max|X(f)|
// 
// Usage: proxy rapide d'intensité spectrale (pic principal)
// ═══════════════════════════════════════════════════════════════════

float PhiScientificModel::_dominantAmplitude(const float* mag) {
    float maxAmp = 0;
    for (size_t i = 1; i < _fftSize/2; i++) {
        if (mag[i] > maxAmp)
            maxAmp = mag[i];
    }
    return maxAmp;
}

// ═══════════════════════════════════════════════════════════════════
// COUPLAGE VIBRO-ACOUSTIQUE (normalisé explicitement)
// ═══════════════════════════════════════════════════════════════════
// DÉFINITION PHYSIQUE:
// Le couplage mesure la cohérence entre vibrations mécaniques (IMU)
// et variations de pression acoustique (micro).
// 
// MÉTHODE:
// 1. Corrélation fréquentielle normalisée [0,1]
//    → Proximité fréquences dominantes
//    → Décroissance exponentielle: exp(-Δf/σ)
//    → Constante σ = 100 Hz (largeur spectrale typique phonème)
// 
// 2. Corrélation énergétique normalisée [0,1]
//    → Moyenne géométrique des amplitudes relatives
//    → e1, e2 normalisés par référence fixe AVANT produit
//    → Référence = 0.1 (amplitude FFT typique signal moyen)
// 
// 3. Score final [0,1]: produit des 2 corrélations
//    → 0 = découplé (fréquences différentes OU énergies faibles)
//    → 1 = couplage parfait (même fréquence + amplitudes fortes)
// 
// NORMALISATION EXPLICITE:
// • Δf normalisé par σ=100 Hz → exp(-Δf/100) ∈ [0,1]
// • e1, e2 normalisés par REF_AMPLITUDE → (e/0.1) ∈ [0,∞[
// • Saturation à 1.0 après normalisation → garantie [0,1]
// • Produit [0,1] × [0,1] ∈ [0,1] garanti
// ═══════════════════════════════════════════════════════════════════

float PhiScientificModel::_computeCorrelation(float f1, float f2, float e1, float e2) {
    // ───────────────────────────────────────────────────────────────
    // 1. Corrélation fréquentielle
    // ───────────────────────────────────────────────────────────────
    const float FREQ_SIGMA = 100.0f;  // Hz - largeur spectrale phonème
    float df = fabsf(f1 - f2);
    float freqCorr = expf(-df / FREQ_SIGMA);
    freqCorr = constrain(freqCorr, 0.0f, 1.0f);
    
    // ───────────────────────────────────────────────────────────────
    // 2. Corrélation énergétique (normalisée EXPLICITEMENT)
    // ───────────────────────────────────────────────────────────────
    // Référence: amplitude FFT typique d'un signal "moyen"
    // Calibrée expérimentalement: voix normale → ~0.05-0.15
    const float REF_AMPLITUDE = 0.1f;
    
    // Normalisation des amplitudes par référence
    float e1_norm = e1 / REF_AMPLITUDE;
    float e2_norm = e2 / REF_AMPLITUDE;
    
    // Moyenne géométrique (insensible échelle mais normalisée)
    float energyCorr = sqrtf(e1_norm * e2_norm);
    
    // Saturation à 1.0 (signaux très forts donnent corrélation max)
    energyCorr = constrain(energyCorr, 0.0f, 1.0f);
    
    // ───────────────────────────────────────────────────────────────
    // 3. Score final [0,1]
    // ───────────────────────────────────────────────────────────────
    return freqCorr * energyCorr;
}

// ═══════════════════════════════════════════════════════════════════
// PLATITUDE SPECTRALE (Spectral Flatness Measure)
// ═══════════════════════════════════════════════════════════════════
// DÉFINITION PHYSIQUE:
// Ratio moyenne géométrique / moyenne arithmétique du spectre de puissance
// 
// INTERPRÉTATION:
// • SF → 1 : spectre plat (bruit blanc) → EMI typique
// • SF → 0 : spectre concentré (tonales) → son structuré typique
// 
// FORMULE:
// SF = exp(⟨log|X(f)|⟩) / ⟨|X(f)|⟩
// 
// Robustesse: floor 1e-10 évite log(0), bins nuls ignorés
// ═══════════════════════════════════════════════════════════════════

float PhiScientificModel::_computeSpectralFlatness(const float* spectrum, size_t n) {
    float logSum = 0.0f;
    float arithmeticSum = 0.0f;
    uint16_t validBins = 0;
    
    const float FLOOR = 1e-10f;
    
    for (size_t i = 1; i < n; i++) {
        float mag = spectrum[i];
        
        if (mag > FLOOR) {
            logSum += log(mag);
            arithmeticSum += mag;
            validBins++;
        }
    }
    
    if (validBins < 2)
        return 0.0f;
    
    float logMean = logSum / validBins;
    float geometricMean = exp(logMean);
    float arithmeticMean = arithmeticSum / validBins;
    
    if (arithmeticMean < FLOOR)
        return 0.0f;
    
    float flatness = geometricMean / arithmeticMean;
    return constrain(flatness, 0.0f, 1.0f);
}

// ═══════════════════════════════════════════════════════════════════
// CLASSIFICATION AUDIO vs EMI (décision purement spectrale)
// ═══════════════════════════════════════════════════════════════════
// PRINCIPE PHYSIQUE:
// Un microphone capte indifféremment:
// • Variations acoustiques structurées (voix, musique) → SF faible
// • Interférences électromagnétiques acoustiques (secteur, switching) → SF élevé
// 
// MÉTHODE:
// Classification basée UNIQUEMENT sur platitude spectrale (SF)
// Le seuil temporel (peak dans PhiAudio) sert uniquement à détecter silence numérique
// 
// SEUILS CALIBRÉS:
// • SF > 0.45 → EMI (bruit blanc/rose dominant)
// • SF < 0.40 → AUDIO (harmoniques structurées)
// • 0.40 < SF < 0.45 → hystérèsis (conserve état précédent)
// 
// ÉNERGIE THRESHOLD:
// • Amplitude < 0.002 → bruit thermique micro → classé EMI (bruit ambiant)
// ═══════════════════════════════════════════════════════════════════

void PhiScientificModel::_classifyAudio() {
    _spectralFlatness = _computeSpectralFlatness(_au_real, _fftSize / 2);
    
    const float FLATNESS_EMI_HIGH = 0.45f;    // Seuil EMI (bruit blanc)
    const float FLATNESS_AUDIO_LOW = 0.40f;   // Seuil Audio (harmoniques)
    const float AMPLITUDE_THRESHOLD = 0.002f; // Plancher bruit thermique
    
    // ───────────────────────────────────────────────────────────────
    // CAS 1: Amplitude sous plancher bruit → bruit thermique ambiant
    // ───────────────────────────────────────────────────────────────
    if (_eAUD < AMPLITUDE_THRESHOLD) {
        _audioType = PhiAudio::TYPE_EMI;
        _previousAudioType = PhiAudio::TYPE_EMI;
        return;
    }
    
    // ───────────────────────────────────────────────────────────────
    // CAS 2: Amplitude détectable → classification spectrale pure
    // ───────────────────────────────────────────────────────────────
    // Hystérèsis à 2 seuils pour stabilité (évite oscillations)
    
    if (_previousAudioType == PhiAudio::TYPE_EMI) {
        // État EMI → basculer AUDIO si SF nettement structuré
        if (_spectralFlatness < FLATNESS_AUDIO_LOW) {
            _audioType = PhiAudio::TYPE_AUDIO;
        } else {
            _audioType = PhiAudio::TYPE_EMI;
        }
    } else {
        // État AUDIO → basculer EMI si SF nettement plat
        if (_spectralFlatness > FLATNESS_EMI_HIGH) {
            _audioType = PhiAudio::TYPE_EMI;
        } else {
            _audioType = PhiAudio::TYPE_AUDIO;
        }
    }
    
    _previousAudioType = _audioType;
}

// ═══════════════════════════════════════════════════════════════════
// SPECTRE 8 BANDES HARMONIQUES (0-1000 Hz)
// ═══════════════════════════════════════════════════════════════════
// Découpage aligné résolution FFT (125 Hz/bin @ 16 kHz, N=128)
// Bandes centrées sur harmoniques naturelles: 50, 125, 250...875 Hz
// 
// Normalisation FIXE (pas adaptive) pour cohérence scientifique:
// → Même niveau bruit ambiant toujours visible
// → Pas de "gonflage" artificiel signaux faibles
// ═══════════════════════════════════════════════════════════════════

void PhiScientificModel::getSpectrum8Bands(uint8_t* bands) const {
    const uint16_t bandRanges[9] = {
        0, 125, 250, 375, 500, 625, 750, 875, 1000
    };
    
    const float hzPerBin = _fsAUD / (float)_fftSize;  // 125 Hz/bin
    const float FIXED_SCALE = 50.0f;  // Échelle calibrée expérimentalement
    
    for (uint8_t b = 0; b < 8; b++) {
        uint16_t binStart = (uint16_t)(bandRanges[b] / hzPerBin);
        uint16_t binEnd = (uint16_t)(bandRanges[b + 1] / hzPerBin);
        
        binStart = constrain(binStart, 1, _fftSize / 2 - 1);
        binEnd = constrain(binEnd, binStart + 1, _fftSize / 2);
        
        // Énergie moyenne bande
        float bandEnergy = 0.0f;
        for (uint16_t i = binStart; i < binEnd; i++) {
            bandEnergy += _au_real[i];
        }
        bandEnergy /= (binEnd - binStart);
        
        // Normalisation fixe + compression logarithmique douce
        float normalized = bandEnergy * FIXED_SCALE;
        if (normalized > 1.0f) {
            normalized = 1.0f + log10f(normalized) / 2.0f;
        }
        
        bands[b] = (uint8_t)constrain((int)(normalized * 255.0f), 0, 255);
    }
}

// ═══════════════════════════════════════════════════════════════════
// TRAITEMENT PRINCIPAL
// ═══════════════════════════════════════════════════════════════════

void PhiScientificModel::process(
    const float* imuBuf,
    const float* audioBuf
) {
    // Application fenêtre Hann + chargement buffers FFT
    for (size_t i = 0; i < _fftSize; i++) {
        _im_real[i] = imuBuf[i] * _hann[i];
        _im_imag[i] = 0.0f;
        _au_real[i] = audioBuf[i] * _hann[i];
        _au_imag[i] = 0.0f;
    }

    // Transformées de Fourier
    _fftIMU->compute(FFT_FORWARD);
    _fftIMU->complexToMagnitude();
    _fftAUD->compute(FFT_FORWARD);
    _fftAUD->complexToMagnitude();

    // Extraction métriques spectrales
    _fIMU = _dominantFreq(_im_real, _fsIMU);
    _eIMU = _dominantAmplitude(_im_real);  
    _fAUD = _dominantFreq(_au_real, _fsAUD);
    _eAUD = _dominantAmplitude(_au_real); 

    // Classification Audio/EMI (décision purement spectrale)
    _classifyAudio();

    // Couplage vibro-acoustique (normalisé explicitement [0,1])
    _corr = _computeCorrelation(_fIMU, _fAUD, _eIMU, _eAUD);
}