Analyse spectrale

[XT95]

Bonjour a tous,
j'aurais besoin d'un petit coup de main pour calculer le spectre audio à la shadertoy.
Apres moult recherches ici et la sur le net, j'ai fait ceci :

- Récupération du flux audio (mise à l'échelle [-1,1]), application d'un fenetrage de Hann.
- Calcule d'une transformée de fourier (rapide) avec en entrée une suite de nombre complexe (réel = flux audio, imaginaire = 0)
- Pour chaque fréquence (sortie de la FFT)
  -- Calcul de la magnitude en utilisant le module du nombre complexe
  -- Calcule de l'amplitude en prenant le log10 de la magnitude et bound de l'amplitude entre [0, 1]

Et enfin j'utilise mon tableau d'amplitude pour afficher des bars de hauteur plus ou moins haute.

C'est ceux qu'ils font ici, et çà fonctionne très bien pour eux.
En revanche je n'obtient pas du tout le même résultat..
La seule différence que j'ai trouvé, est qu'on utilise pas la même bibliothèque pour calculer la FFT, a par ca j'ai bien l'impression de faire la même chose qu'eux  .

Voila, je viens de passer ma soirée dessus, je craque !!
Peut être qu'il manque une étape que je ne vois pas, quelqu'un aurait il une idée (Ponce ? Xtrium ?  )
Je poste mon code au cas ou vous trouvez quelque chose :

Code: [Sélectionner]

void SoundManager::run()
{
    fftw_complex *in1, *out1;
    fftw_plan p1;
    in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * SOUND_DATA_SIZE);
    out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * SOUND_DATA_SIZE);
    p1 = fftw_plan_dft_1d(SOUND_DATA_SIZE, in1, out1, FFTW_FORWARD, FFTW_ESTIMATE);


    ALshort samples[SOUND_DATA_SIZE];
    while(true)
    {
        if(m_captureDevice)
        {
            ALCint sampleSize;
            alcGetIntegerv(m_captureDevice, ALC_CAPTURE_SAMPLES, 1, &sampleSize);
            if(sampleSize > SOUND_DATA_SIZE)
            {
                alcCaptureSamples(m_captureDevice, samples, SOUND_DATA_SIZE);

                for(int i=0; i<SOUND_DATA_SIZE; i++)
                {
                    in1[i][0] = double(samples[i])/32768.0 * 0.5 *(1-cos(2*M_PI*i/(SOUND_DATA_SIZE-1))); //Hann Window
                    in1[i][1] = 0.;
                }
                fftw_execute(p1);
                emit newSpectrumDataChannel1(out1);
            }
        }
        this->msleep(5);
    }
    fftw_destroy_plan(p1);
    fftw_free(in1); fftw_free(out1);
}

Code: [Sélectionner]

void RenderWidget::updateTextureSpectrumChannel1(fftw_complex *out)
{
    for(int i=2; i<=SOUND_DATA_SIZE/2; i++)
    {
        float magnitude = sqrt(out[i][0]*out[i][0]+out[i][1]*out[i][1]);
        float amplitude = 0.15 * log(magnitude);
        amplitude = fmax(amplitude,0.0);
        amplitude = fmin(amplitude,1.0);
        m_dataSpectrum[i] = amplitude;
    }
}

Mercii !

[Patapom]

J'y connais rien mais juste une remarque : la magnitude d'un nombre complexe avec partie imaginaire nulle c'est juste la valeur absolue de de la partie réelle hein...

[XT95]

Non je me suis mal exprimé!
Je calcul la magnitude a partir du tableau de complexes "out" qui est en sortie de la FFT, et donc la partie imaginaire n'est plus nulle .

[ponce]

En revanche je n'obtient pas du tout le même résultat..

Tu obtiens quoi justement ? Tu peux nous montrer une image de \( \lvert FFT(x) \rvert \) ?

En première approximation je dirais que tu dois obtenir le spectre de puissance mais ça va être un spectre "miroir" seule la première moitié t'intéresse (la deuxième moitié est redondante parce que ton input est réel) (EDiT: ah ok tu l'as vu).
Ensuite il y l'histoire du scaling, des fois tu as besoin de diviser par N après une FFT ou IFFT et je sais jamais trop où 

[ponce]

Sinon c'est bien de fenétrer mais si tu ne fais pas d'overlap des segments fenetrés tu va perdre de l'information ! Je posterai ce soir du code qui fait ça correctement.

[ponce]

Code: [Sélectionner]

#include <complex>
#include <vector>
#include <algorithm>
#include <assert.h>

#define _USE_MATH_DEFINES
#include <math.h>

// Return true if power of 2.
inline bool isPowerOf2(size_t i)
{
    assert(i >= 0);
    return (i != 0) && ((i & (i - 1)) == 0);
}

// return x so that (1 << x) >= i
int iFloorLog2(size_t i)
{
    assert(i > 0);
    
    #if defined(_WIN32)
        unsigned long y;
        _BitScanReverse(&y, i);
        return y;
    #else            
        int result = 0;
        while (i > 1)
        {
            i = i / 2;
            result = result + 1;
        }
        return result;
    #endif       
}

class Window
{
public:
    enum Type
    {
        TYPE_RECT,
        TYPE_HANN
    };

    template<typename T>
    static void generate(Type type, size_t N, T out[])
    {
        for (size_t i = 0; i < N; ++i)
        {
            out[i] = static_cast<T>(eval(type, i, N));
        }
    }

    static double eval(Type type, size_t n, size_t N)
    {
        switch(type)
        {
        case TYPE_RECT:
            return 1.0;

        case TYPE_HANN:
            return 0.5 - 0.5 * cos((2 * M_PI * n) / (N - 1));

        default:
            assert(false);
        }
    }
};


enum FFTDirection
{
    FFT_FORWARD = 0,
    FFT_REVERSE = 1
};

template<typename T>
void FFT(std::complex<T>* inout, size_t size, FFTDirection direction)
{
    assert(isPowerOf2(size));
    int m = iFloorLog2(size);

    // do the bit reversal
    int i2 = size / 2;
    int j = 0;
    for (int i = 0; i < (int)size - 1; ++i)
    {
        if (i < j)
            std::swap(inout[i], inout[j]);

        int k = i2;
        while(k <= j)
        {
            j = j - k;
            k = k / 2;
        }
        j += k;
    }

    // compute the FFT
    std::complex<T> c = -1;
    int l2 = 1;
    for (int l = 0; l < m; ++l)
    {
        int l1 = l2;
        l2 = l2 * 2;
        std::complex<T> u = 1;
        for (int j2 = 0; j2 < l1; ++j2)
        {
            int i = j2;
            while (i < (int)size)
            {
                int i1 = i + l1;
                std::complex<T> t1 = u * inout[i1];
                inout[i1] = inout[i] - t1;
                inout[i] += t1;
                i += l2;
            }
            u = u * c;
        }

        T newImag = sqrt((1 - real(c)) / 2);
        if (direction == FFT_FORWARD)
            newImag = -newImag;
        T newReal = sqrt((1 + real(c)) / 2);
        c = std::complex<T>(newReal, newImag);
    }

    // scaling for forward transformation
    if (direction == FFT_FORWARD)
    {
        for (int i = 0; i < (int)size; ++i)
            inout[i] = inout[i] / std::complex<T>((T)size, 0);
    }
}

// Send Short term FFT data to a listener.
// Variable overlap.
// Introduced approx. windowSize/2 samples delay.
class FFTAnalyzer
{
public:

    // if zeroPhaseWindowing = true, "zero phase" windowing is used
    // (center of window is at first sample, zero-padding happen at center)
    FFTAnalyzer(Window::Type windowType, bool zeroPhaseWindowing, bool correctWindowLoss)
    {
        _windowType = windowType;
        _zeroPhaseWindowing = zeroPhaseWindowing;
        _setupWasCalled = false;
        _correctWindowLoss = correctWindowLoss;
    }

    size_t windowSize() const
    {
        return _windowSize;
    }

    size_t analysisPeriod() const
    {
        return _analysisPeriod;
    }

    // To call at initialization and whenever samplerate changes
    // windowSize = size of analysis window, expressed in samples
    // fftSize = size of FFT. Must be power-of-two and >= windowSize. Missing samples are zero-padded in time domain.
    // analysisPeriod = period of analysis results, allow to be more precise frequentially, expressed in samples
    // Basic overlap is achieved with windowSize = 2 * analysisPeriod
    void setup(size_t windowSize, size_t fftSize, size_t analysisPeriod)
    {
        assert(isPowerOf2(fftSize));
        assert(fftSize >= windowSize);
        _setupWasCalled = true;

        assert(windowSize != 1);
        assert(analysisPeriod <= windowSize); // no support for zero overlap

        // 1-sized FFT support
        if (analysisPeriod == 0)
            analysisPeriod = 1;

        _windowSize = windowSize;
        _fftSize = fftSize;
        _analysisPeriod = analysisPeriod;

        // clear input delay
        _audioBuffer.resize(_windowSize);
        _index = 0;

        _windowBuffer.resize(_windowSize);
        Window::generate(_windowType, _windowSize, _windowBuffer.data());

        _windowGainCorrFactor = 0;
        for (size_t i = 0; i < _windowSize; ++i)
            _windowGainCorrFactor += _windowBuffer[i];
        _windowGainCorrFactor = _windowSize / _windowGainCorrFactor;

        if (_correctWindowLoss)
        {
            for (size_t i = 0; i < _windowSize; ++i)
                _windowBuffer[i] *= _windowGainCorrFactor;
        }
            
    }

    // Process one sample, eventually return the result of short-term FFT
    // in a given Buffer
    bool feed(float x, std::vector<std::complex<float>>* fftData)
    {
        assert(_setupWasCalled);
        _audioBuffer[_index] = x;
        _index = _index + 1;
        if (_index >= _windowSize)
        {
            fftData->resize(_fftSize);

            if (_zeroPhaseWindowing)
            {
                // "Zero Phase" windowing
                // Through clever reordering, phase of ouput coefficients will relate to the
                // center of the window
                //_
                // \_                   _/
                //   \                 /
                //    \               /
                //     \_____________/____
                size_t center = (_windowSize - 1) / 2; // position of center bin
                size_t nLeft = _windowSize - center;
                for (size_t i = 0; i < nLeft; ++i)
                    (*fftData)[i] = _audioBuffer[center + i] * _windowBuffer[center + i];

                size_t nPadding = _fftSize - _windowSize;
                for (size_t i = 0; i < nPadding; ++i)
                    (*fftData)[nLeft + i] = 0.0f;

                for (size_t i = 0; i < center; ++i)
                    (*fftData)[nLeft + nPadding + i] = _audioBuffer[i] * _windowBuffer[i];
            }
            else
            {
                // "Normal" windowing
                // Phase of ouput coefficient will relate to the start of the buffer
                //      _
                //    _/ \_
                //   /     \
                //  /       \
                //_/         \____________

                // fill FFT buffer and multiply by window
                for (size_t i = 0; i < _windowSize; ++i)
                    (*fftData)[i] = _audioBuffer[i] * _windowBuffer[i];

                // zero-padding
                for (size_t i = _windowSize; i < _fftSize; ++i)
                    (*fftData)[i] = 0.0f;
            }

            // perform forward FFT on this slice
            FFT(fftData->data(), _fftSize, FFT_FORWARD);

            // rotate buffer
            {
                size_t const samplesToDrop = _analysisPeriod;
                assert(0 < samplesToDrop && samplesToDrop <= _windowSize);
                size_t remainingSamples = _windowSize - samplesToDrop;

                // TODO: use ring buffer instead of copy
                memmove(_audioBuffer.data(), _audioBuffer.data() + samplesToDrop, sizeof(float) * remainingSamples);
                _index = remainingSamples;

            }
            return true;
        }
        else
        {
            return false;
        }
    }

private:
    std::vector<float> _audioBuffer;
    std::vector<float> _windowBuffer;

    size_t _fftSize;        // in samples
    size_t _windowSize;     // in samples
    size_t _analysisPeriod; // in samples

    Window::Type _windowType;
    bool _zeroPhaseWindowing;

    size_t _index;
    bool _setupWasCalled;

    // should we multiply by _windowGainCorrFactor?
    bool _correctWindowLoss;
    // the factor by which to multiply transformed data to get in range results
    float _windowGainCorrFactor; 
};

Avec du windowing zero-phase optionnel pour que la phase que tu obtienne corresponde au milieu du boût de signal signal fenêtré. Je pense que la correction du gain pour pallier au fenêtrage ne marche pas par contre.


EDIT: et comme XT me dit que le FFTAnalyzer est chaud à configurer, voilà comment faire

Code: [Sélectionner]

double desiredAnalysisWindowInMs = 20; // pour une fenêtre de 20 ms
size_t windowSizeInSamples = (size_t)(0.5 + desiredAnalysisWindowInMs * sampleRate / 1000);                       
if (!impair(windowSizeInSamples))
  ++windowSizeInSamples;
int FFT_OVERSAMPLE = 1; // mettre plus pour avoir une FFT lissée mais sans info supplémentaire
int fftSize = nextPowerOf2(windowSizeInSamples) * FFT_OVERSAMPLE;
int analysisPeriodInSamples = (windowSizeInSamples + 1) / 2;  // 50% overlap mais vous pouvez diviser plus pour analyser plus
fftAnalyzer.setup(windowSizeInSamples, fftSize, analysisPeriodInSamples);

[Patapom]

 \( \lvert FFT(x) \rvert \) ?

Genre on a du \( L_aT^ex \) maintenant ?!!

[ponce]

Eh oui ! Depuis http://bbs.demoscene.fr/aide/modifications-du-forum/msg13178/#msg13178

Genre tu écris

Code: [Sélectionner]

 \(\sqrt{3x-1}+(1+x)^2\)

et ça donne

\(\sqrt{3x-1}+(1+x)^2\)

[Zavie]

Genre on a du \( L_aT^ex \) maintenant ?!!

Tu voulais sûrement dire \(\LaTeX\).
(en vrai ceci était un test ; ponce tu roxes !)

[ponce]

Ah ben c'est pas moi :p c'est MathJax (http://www.mathjax.org/) qui en plus prêtent gracieusement leur CDN. C'est très facile à intégrer dans un site en fait.

[Patapom]

Comment je fais pour l'utiliser pour générer des images depuis PHP ?? C'est possible ? Faut que je RTFM ?

[ponce]

Je sais pas trop mais apparemment tu peux self-hoster le truc, j'ai pas compris comment.

[Zavie]

C'est cool en plus ça génère un vrai texte sélectionnable et non pas une malheureuse image. Merci pour le tuyau !

[ulrick]

J'savais pas qu'on pouvait frimer avec des $$a^{a^a}.$$ C'est cool car j'ai fait ma thèse en \LaTeX $$q^\ast=\underset{q}{\textit argmax}\{U(q)-pq\},$$ s.t. $$pq\le r.$$ Mais maintenant, j'aime bien Word, je frime plus .