Estrazione frequenze precise da FFT scomparti con cambiamento di fase tra fotogrammi

https://stackoverflow.com/questions/4633203

08-10-2019
|

Domanda

Ho cercato attraverso questo fantastico articolo: http: / /blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

Pur essendo fantastico, è estremamente duro e pesante in corso. Questo materiale è veramente mi stretching.

Ho estratto i calcoli del modulo di codice di Stefan che calcola la frequenza per un assegnato bin. Ma io non capisco l'ultimo calcolo. Qualcuno può spiegare a me la costruzione matematica alla fine?

Prima di scavare nel codice, mi permetta di impostare la scena:

Diciamo che abbiamo impostato fftFrameSize = 1024, quindi abbiamo a che fare con 512 + 1 scomparti
Come esempio, Bin [1] 's frequenza ideale inserisce una singola onda nel telaio. Ad una frequenza di campionamento di 40 KHz, tOneFrame = 1024 / 40K secondi = 1 / 40s, così Bin [1] andrebbe idealmente raccogliendo un segnale 40Hz.
Impostazione osamp (sovracampionamento) = 4, progrediamo lungo il segnale di ingresso a passi di 256. Così i primi esamina analisi byte zero a 1023, quindi da 256 a 1279, ecc Si noti ciascun galleggiante viene processato 4 volte.

...

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 :) Ora il bit non ottengo:

        // 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;
    }
}

non riesco proprio a vederlo chiaramente, anche se sembra essere guardando in faccia. Qualcuno potrebbe spiegare questo processo da zero, passo dopo passo?

Soluzione 4

Finalmente ho capito questo; veramente ho dovuto ricavare da zero. Sapevo che ci sarebbe stato un modo semplice per ricavare, il mio (al solito) errore è stato quello di tentare di seguire la logica di altre persone, piuttosto che utilizzare il mio buon senso.

Questo puzzle ha due tasti per sbloccarlo.

La prima chiave è di capire come over-sampling introduce una rotazione su fase di bin.
La seconda chiave viene da Graph 3.3 e 3.4 qui: 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;
}

Altri suggerimenti

Il principio di base è molto semplice. Se un dato componente corrisponde esattamente una frequenza bin allora la sua fase non cambierà da un FT a quella successiva. Tuttavia, se la frequenza non corrisponde esattamente alla frequenza bin allora ci sarà un cambiamento di fase tra FTS successive. Il delta di frequenza è solo:

delta_freq = delta_phase / delta_time

e la stima raffinato della frequenza della componente saranno quindi:

freq_est = bin_freq + delta_freq

Ho implementato questo algoritmo per Performous me stesso. Quando si prende un'altra FFT con un offset di tempo, si prevede che la fase di cambiare secondo l'offset, ossia due FFT prelevati 256 campioni a parte dovrebbe avere una differenza di fase di 256 campioni per tutte le frequenze presenti nel segnale (ciò presuppone che i segnali stessi sono costante, che è un buon presupposto per periodi brevi come 256 campioni).

Ora, i valori di fase effettivi si ottiene da FFT non sono in campioni ma in angolo di fase, in modo che saranno diverse a seconda della frequenza. Nel codice seguente il valore phaseStep è il fattore di conversione necessaria per intervallo, cioè per frequenza corrispondente x bin, lo sfasamento sarà x * phaseStep. Per frequenze centrali bin x sarebbe un numero intero (il numero bin) ma per frequenze reali rilevate può essere qualsiasi numero reale.

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

La correzione funziona ipotizzando che il segnale in un bidone ha la frequenza centrale bin e poi calcolando lo sfasamento previsto per questo. Questo cambiamento atteso è sottratto dal trasferimento effettivo, lasciando l'errore. Un resto (modulo 2 pi) è presa (-pi a gamma pi) e la frequenza finale viene calcolato con cuore bin + correzione.

// 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

Si noti che molti bidoni adiacenti spesso finiscono corretto alla stessa frequenza, poiché la correzione delta può essere fino a 0,5 bidoni * FFT_N / FFT_STEP entrambe le direzioni in modo che il FFT_STEP minore si utilizza, più lontano sarà correzione possibile (ma questo aumenta la potenza di elaborazione necessaria nonché imprecisioni dovute a imprecisioni).

Spero che questo aiuta:)

Questa è la tecnica di stima di frequenza utilizzata dai metodi fase vocoder.

Se si guarda a un singolo punto su uno (frequenza fissa e ampiezza fissa) dell'onda sinusoidale nel tempo, la fase avanzerà nel tempo di una quantità proporzionale alla frequenza. Oppure si può fare il contrario:. Se si misura quanto la fase di un cambiamento sinusoide oltre qualsiasi unità di tempo, è possibile calcolare la frequenza di quel sinusoide

Un vocoder fase utilizza due FFT per stimare fase con riferimento a due finestre FFT, e l'offset dei due FFT è la distanza tra le misurazioni 2 di fase nel tempo. Di là, avete la vostra stima di frequenza per che bin FFT (uno scomparto di FFT di essere più o meno un filtro per isolare un componente sinusoidale o altro segnale abbastanza a banda stretta che si inserisce all'interno che bin).

Per questo metodo di lavoro, lo spettro vicino al bidone FFT in uso deve essere abbastanza stazionario, ad esempio non cambiando in frequenza, ecc Questo è il presupposto di un vocoder fase richiede.

Forse questo aiuterà. Pensate ai bidoni FFT come specificando piccoli orologi o rotori, ciascuno di filatura alla frequenza del bidone. Per un segnale stabile, il (teorico) successiva posizione del rotore può essere previsto utilizzando la matematica nel bit non si ottiene. Contro questa posizione "dovrebbe essere" (ideale), è possibile calcolare parecchie cose utili: (1) la differenza con la fase in un bidone di un telaio adiacente, che viene utilizzato da un fase di vocoder per meglio stima della frequenza bin, o (2) di più generalmente deviazione di fase , che è un indicatore positivo di un esordio nota o qualche altro evento audio.

frequenze di segnale che cadono esattamente su una fase avanzata bidone frequenza bin di multipli interi di 2p. Poiché le fasi bin che corrispondono alle frequenze bin sono multipli di 2p dovute alla natura periodica della FFT non c'è variazione di fase in questo caso. L'articolo si parla anche lo spiega.

Autorizzato sotto: CC-BY-SA insieme a attribuzione

Non affiliato a StackOverflow