La extracción de frecuencias precisas de FFT Bins utilizando el cambio de fase entre tramas

https://stackoverflow.com/questions/4633203

08-10-2019
|

Pregunta

He estado buscando a través de este fantástico artículo: http: / /blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

Si bien es fantástica, es extremadamente duro y pesado ir. Este material me está estirando.

Me han extraído las matemáticas de módulo de código de Stefan que calcula la frecuencia exacta para un contenedor dado. Pero no entiendo el último cálculo. ¿Puede alguien explicarme la construcción matemática al final?

Antes de la excavación en el código, vamos a poner la escena:

Supongamos que establecemos fftFrameSize = 1,024, por lo que estamos tratando con 512 contenedores + 1
Como ejemplo, Bin [1] 's frecuencia ideal se ajusta a una sola onda en el marco. A una frecuencia de muestreo de 40 KHz, tOneFrame = 1024 / 40K segundos = 1/40, de modo Bin [1] idealmente ser recogida de una señal de 40 Hz.
Configuración de osamp (sobremuestreo) = 4, que el progreso a lo largo de nuestra señal de entrada en pasos de 256. Así que la primera examina de análisis cero bytes a través de 1023, a continuación, 256 a través de 1279, etc. Nota cada flotador se procesa 4 veces.

...

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 :) Ahora el bit que no entiendo:

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

Yo simplemente no puede ver con claridad, a pesar de que parece estar mirando a la cara. Podría alguien explicar este proceso desde cero, paso a paso?

Solución 4

Finalmente he dado cuenta de esto; Realmente tuve que derivar a partir de cero. Yo sabía que habría alguna forma sencilla de una derivación del mismo, mi (habitual) error fue tratar de seguir la lógica de otras personas en lugar de usar mi propio sentido común.

Este puzzle tiene dos teclas para desbloquearlo.

La primera clave es entender cómo introduce sobre-muestreo de una rotación en fase de bin.
La segunda clave proviene de Gráfico 3.3 y 3.4 aquí: 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;
}

Otros consejos

El principio básico es muy simple. Si un determinado componente coincide exactamente con una frecuencia bin entonces su fase no cambiará de un FT a la siguiente. Sin embargo, si la frecuencia no se corresponde exactamente con la frecuencia bin entonces habrá un cambio de fase entre FTS sucesivas. El delta de frecuencia es simplemente:

delta_freq = delta_phase / delta_time

y la estimación refinada de la frecuencia de la componente serán entonces:

freq_est = bin_freq + delta_freq

Me han implementado este algoritmo para Performous mí mismo. Cuando se toma otra FFT a una compensación de tiempo, se espera que la fase de cambio de acuerdo a la compensación, es decir, dos FFT tomado 256 muestras, aparte debe tener una diferencia de fase de 256 muestras para todas las frecuencias presentes en la señal (esto supone que las señales propias son constantes, que es una buena suposición para períodos cortos como 256 muestras).

Ahora, los valores de fase reales que obtiene de la FFT no están en muestras pero en ángulo de fase, por lo que será diferente en función de la frecuencia. En el siguiente código el valor phaseStep es el factor de conversión necesaria por bin, es decir, para la frecuencia correspondiente a x bin, el desplazamiento de fase será x * phaseStep. Para frecuencias centrales bin x sería un número entero (el número bin) pero para frecuencias detectadas reales puede ser cualquier número real.

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

La corrección funciona suponiendo que la señal en un cubo tiene la frecuencia central bin y luego calculando el cambio de fase esperada para eso. Este cambio esperado se sustrajo el cambio real, dejando el error. Un resto (módulo 2 pi) se toma (-pi a la gama pi) y la frecuencia final se calcula con bin centro + corrección.

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

Tenga en cuenta que muchos contenedores adyacentes a menudo terminan corregido a la misma frecuencia porque la corrección delta puede ser de hasta 0,5 papeleras * FFT_N / FFT_STEP de cualquier manera por lo que el FFT_STEP más pequeña que utiliza, cuanto más lejos será posible correcciones (pero esto aumenta la potencia de procesamiento necesaria, así como la imprecisión debido a imprecisiones).

Espero que esto ayude:)

Esta es la técnica de estimación de frecuencia utilizada por los métodos de vocoder de fase.

Si nos fijamos en un solo punto en un (frecuencia y amplitud fijas fijo) de onda sinusoidal en el tiempo, la fase avanzará con el tiempo en una cantidad proporcional a la frecuencia. O puede hacerlo a la inversa:. Si se mide la cantidad de la fase de una sinusoide cambios a lo largo de cualquier unidad de tiempo, se puede calcular la frecuencia de esa sinusoide

Un vocoder de fase utiliza dos FFT para estimar fase con referencia a dos ventanas de FFT, y el desplazamiento de las dos FFT es la distancia entre las 2 mediciones de fase en el tiempo. A partir de ahí, usted tiene su estimación de la frecuencia para que bin FFT (FFT siendo una bandeja de más o menos un filtro para aislar un componente sinusoidal u otra señal suficientemente estrecha que se ajusta dentro de que Bin).

Para que este método de trabajo, el espectro cerca del bin FFT en uso tiene que ser bastante estacionaria, por ejemplo no cambiar de frecuencia, etc. Esa es la suposición de un vocoder fase requiere.

Tal vez esto ayude. Piense en las casillas FFT como la especificación de pequeños relojes o rotores, cada uno de hilado a la frecuencia de la bin. Para una señal estable, la siguiente posición (teórica) del rotor puede predecirse utilizando las matemáticas en el bit no obtiene. Contra esta posición "debe ser" (ideal), puede calcular varias cosas útiles: (1) la diferencia con la fase en un compartimiento de un marco adyacente, que es utilizado por un fase vocoder a una mejor estimación de la frecuencia bin, o (2) de manera más general desviación de fase , que es un indicador positivo de un inicio nota o algún otro evento en el audio.

frecuencias de señal que caen exactamente sobre una fase bin avance frecuencia bin por múltiplos enteros de 2p. Puesto que las fases de basura que corresponden a las frecuencias de las barras son múltiplos de 2p debido a la naturaleza periódica de la FFT no hay ningún cambio de fase en este caso. El artículo que mencionas también explica esto.

Licenciado bajo: CC-BY-SA con atribución

No afiliado a StackOverflow