Extraction des fréquences précises de FFT Bins en utilisant un changement de phase entre trames

https://stackoverflow.com/questions/4633203

08-10-2019
|

Question

Je suis à la recherche à travers cet article fantastique: http: / /blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

Tout en étant fantastique, il est extrêmement difficile et lourd va. Ce matériel est vraiment me stretching.

J'ai extrait les mathématiques à partir du module de code de Stefan qui calcule la fréquence exacte pour un bac donné. Mais je ne comprends pas le dernier calcul. Quelqu'un peut-il me expliquer la construction mathématique à la fin?

Avant de creuser dans le code, permettez-moi de planter le décor:

Disons que nous avons mis en fftFrameSize = 1024, nous avons donc affaire à 512 + 1 bacs
A titre d'exemple, Bin [1] de 'la fréquence idéale correspond une onde unique dans la trame. A une fréquence d'échantillonnage de 40KHz, tOneFrame = 1024 / 40K seconde = 1/40, de sorte Bin [1] serait idéalement recueille un signal de 40Hz.
Réglage osamp (SURÉCHANTILLON) = 4, nous progressons le long de notre signal d'entrée dans les étapes de 256. La première Étudie analyse octets zéro à 1023, puis 256 à 1279, etc. Notez chaque flotteur obtient traité 4 fois.

...

void calcBins( 
              long fftFrameSize, 
              long osamp, 
              float sampleRate, 
              float * floats, 
              BIN * bins
              )
{
    /* initialize our static arrays */
    static float gFFTworksp[2*MAX_FRAME_LENGTH];
    static float gLastPhase[MAX_FRAME_LENGTH/2+1];

    static long gInit = 0;
    if (! gInit) 
    {
        memset(gFFTworksp, 0, 2*MAX_FRAME_LENGTH*sizeof(float));
        memset(gLastPhase, 0, (MAX_FRAME_LENGTH/2+1)*sizeof(float));
        gInit = 1;
    }

    /* do windowing and re,im interleave */
    for (long k = 0; k < fftFrameSize; k++) 
    {
        double window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
        gFFTworksp[2*k] = floats[k] * window;
        printf("sinValue: %f", gFFTworksp[2*k]);
        gFFTworksp[2*k+1] = 0.;
    }

    /* do transform */
    smbFft(gFFTworksp, fftFrameSize, -1);

    printf("\n");

    /* this is the analysis step */
    for (long k = 0; k <= fftFrameSize/2; k++) 
    {
        /* de-interlace FFT buffer */
        double real = gFFTworksp[2*k];
        double imag = gFFTworksp[2*k+1];

        /* compute magnitude and phase */
        double magn = 2.*sqrt(real*real + imag*imag);
        double phase = atan2(imag,real);

        /* compute phase difference */
        double phaseDiff = phase - gLastPhase[k];
        gLastPhase[k] = phase;

        /* subtract expected phase difference */
        double binPhaseOffset = M_TWOPI * (double)k / (double)osamp;
        double deltaPhase = phaseDiff - binPhaseOffset;

        /* map delta phase into [-Pi, Pi) interval */
        // better, but obfuscatory...
        //    deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);

        while (deltaPhase >= M_PI)
            deltaPhase -= M_TWOPI;
        while (deltaPhase < -M_PI)
            deltaPhase += M_TWOPI;

(EDIT :) Maintenant le bit Je ne comprends pas:

        // Get deviation from bin frequency from the +/- Pi interval 
        // Compute the k-th partials' true frequency    

        // Start with bin's ideal frequency
        double bin0Freq = (double)sampleRate / (double)fftFrameSize;
        bins[k].idealFreq = (double)k * bin0Freq;

        // Add deltaFreq
        double sampleTime = 1. / (double)sampleRate;
        double samplesInStep = (double)fftFrameSize / (double)osamp;
        double stepTime = sampleTime * samplesInStep;
        double deltaTime = stepTime;        

        // Definition of frequency is rate of change of phase, i.e. f = dϕ/dt
        // double deltaPhaseUnit = deltaPhase / M_TWOPI; // range [-.5, .5)
        double freqAdjust = (1. / M_TWOPI) * deltaPhase / deltaTime; 

        // Actual freq <-- WHY ???
        bins[k].freq = bins[k].idealFreq + freqAdjust;
    }
}

Je ne peux pas le voir clairement, même si elle semble regarder en face. Quelqu'un pourrait-il expliquer s'il vous plaît ce processus à partir de zéro, étape par étape?

La solution 4

Enfin j'ai compris cela; vraiment je devais tirer à partir de zéro. Je savais qu'il y aurait une façon simple de dériver, mon (habituelle) erreur a été de tenter de suivre la logique des autres plutôt que d'utiliser mon bon sens.

Ce casse-tête prend deux touches pour le déverrouiller.

La première clé est de comprendre comment introduit suréchantillonnage une rotation de la phase bin.
La deuxième clé provient du graphique 3.3 et 3.4 ici : http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/

...

for (int k = 0; k <= fftFrameSize/2; k++) 
{
    // compute magnitude and phase 
    bins[k].mag = 2.*sqrt(fftBins[k].real*fftBins[k].real + fftBins[k].imag*fftBins[k].imag);
    bins[k].phase = atan2(fftBins[k].imag, fftBins[k].real);

    // Compute phase difference Δϕ fo bin[k]
    double deltaPhase;
    {
        double measuredPhaseDiff = bins[k].phase - gLastPhase[k];
        gLastPhase[k] = bins[k].phase;

        // Subtract expected phase difference <-- FIRST KEY
        // Think of a single wave in a 1024 float frame, with osamp = 4
        //   if the first sample catches it at phase = 0, the next will 
        //   catch it at pi/2 ie 1/4 * 2pi
        double binPhaseExpectedDiscrepancy = M_TWOPI * (double)k / (double)osamp;
        deltaPhase = measuredPhaseDiff - binPhaseExpectedDiscrepancy;

        // Wrap delta phase into [-Pi, Pi) interval 
        deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
    }

    // say sampleRate = 40K samps/sec, fftFrameSize = 1024 samps in FFT giving bin[0] thru bin[512]
    // then bin[1] holds one whole wave in the frame, ie 44 waves in 1s ie 44Hz ie sampleRate / fftFrameSize
    double bin0Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin0Freq;

    // Consider Δϕ for bin[k] between hops.
    // write as 2π / m.
    // so after m hops, Δϕ = 2π, ie 1 extra cycle has occurred   <-- SECOND KEY
    double m = M_TWOPI / deltaPhase;

    // so, m hops should have bin[k].idealFreq * t_mHops cycles.  plus this extra 1.
    // 
    // bin[k].idealFreq * t_mHops + 1 cycles in t_mHops seconds 
    //   => bins[k].actualFreq = bin[k].idealFreq + 1 / t_mHops
    double tFrame = fftFrameSize / sampleRate;
    double tHop = tFrame / osamp;
    double t_mHops = m * tHop;

    bins[k].freq = bins[k].idealFreq + 1. / t_mHops;
}

Autres conseils

Le principe de base est très simple. Si un composant donné correspond exactement à une fréquence bin alors sa phase ne changera pas d'un FT à l'autre. Toutefois, si la fréquence ne correspond pas exactement à la fréquence bin alors il y aura un changement de phase entre les tables de tir successifs. Le delta de fréquence est juste:

delta_freq = delta_phase / delta_time

et l'estimation affinée de la fréquence de la composante sera alors:

freq_est = bin_freq + delta_freq

J'ai mis en œuvre cet algorithme pour Performous moi-même . Quand on prend un autre FFT à un décalage de temps, on attend la phase de changer en fonction du décalage, à savoir deux FFT prises 256 échantillons devraient part avoir une différence de phase de 256 échantillons pour toutes les fréquences présentes dans le signal (ce qui suppose que les signaux eux-mêmes sont stables, ce qui est une bonne supposition pour de courtes périodes comme 256 échantillons).

Maintenant, les valeurs de phase réelles que vous obtenez de FFT ne sont pas dans des échantillons, mais l'angle de phase, donc ils seront différents en fonction de la fréquence. Dans le code suivant la valeur de phaseStep est le facteur de conversion nécessaire par bac, à savoir la fréquence correspondant au bac x, le déphasage sera x * phaseStep. Pour les fréquences centrales bin x serait un entier (le numéro bin), mais pour les fréquences détectées réelles, il peut être un nombre réel.

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

La correction agit en supposant que le signal dans un bac a la fréquence centrale du casier, puis en calculant le décalage de phase prévu pour cela. Ce changement attendu est soustrait du changement réel, en laissant l'erreur. Un reste (modulo 2 pi) est prise (-pi à la plage de pi) et la fréquence finale est calculée avec le centre bin + correction.

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

Notez que plusieurs bacs adjacents finissent souvent corrigées à la même fréquence que la correction delta peut être jusqu'à 0,5 * FFT_N / bacs FFT_STEP ou l'autre manière si l'utilisation FFT_STEP plus petit, le plus loin sera possible corrections (mais cela augmente la puissance de traitement nécessaire ainsi que l'imprécision due à des erreurs).

J'espère que cette aide:)

Ceci est la technique d'estimation de fréquence utilisée par des méthodes de vocodeur de phase.

Si vous regardez un seul point sur une (fréquence fixe et d'amplitude fixe) onde sinusoïdale dans le temps, la phase avancera avec le temps d'un montant proportionnel à la fréquence. Ou vous pouvez faire l'inverse. Si vous mesurez combien la phase d'un changement sinusoïde sur une unité de temps, vous pouvez calculer la fréquence de cette sinusoïde

A vocodeur de phase utilise deux FFT pour estimer la phase de référence à deux fenêtres FFT, et le décalage des deux FFT est la distance entre les 2 mesures de phase dans le temps. De là, vous avez votre estimation de fréquence pour cette case FFT (un bac FFT étant à peu près un filtre pour isoler une composante sinusoïdale ou un autre signal suffisamment bas débit qui correspond à l'intérieur de ce bac).

Pour que cette méthode de travail, le spectre près de la case FFT en cours d'utilisation doit être assez stationnaire, par exemple ne change pas en fréquence, etc. C'est l'hypothèse d'un vocodeur de phase nécessite.

Peut-être que cela vous aidera. Pensez aux bacs de FFT en spécifiant peu horloges ou rotors, chaque rotation à la fréquence du bac. Pour un signal stable, le (théorique) position suivante du rotor peut être prédite en utilisant le calcul dans le bit, vous ne recevez pas. Dans ce « devrait être » la position (idéal), vous pouvez calculer plusieurs choses utiles: (1) la différence avec la phase dans un bac d'un cadre adjacent, qui est utilisé par vocodeur de phase pour mieux estimation de la fréquence de casier, ou (2) plus généralement écart de phase , qui est un indicateur positif d'un début de note ou d'un autre événement dans le signal audio.

fréquences de signaux qui se situent exactement sur une fréquence de bin phase de bac à l'avance par des multiples entiers de 2p. Depuis les phases de bin qui correspondent aux fréquences bin sont des multiples de 2p en raison de la nature périodique de la FFT il n'y a pas de changement de phase dans ce cas. L'article que vous mentionnez explique aussi cela.

Licencié sous: CC-BY-SA avec attribution

Non affilié à StackOverflow