como - Extracción de frecuencias precisas de los contenedores FFT utilizando el cambio de fase entre cuadros

como usar fft en matlab

Question

He estado mirando este fantástico artículo: http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

Si bien es fantástico, es extremadamente difícil y pesado. Este material realmente me está estirando.

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

Antes de profundizar en el código, déjame establecer la escena:

• Digamos que configuramos fftFrameSize = 1024, por lo que estamos tratando con 512 + 1 contenedores

• Como ejemplo, la frecuencia ideal de Bin [1] se ajusta a una sola onda en el marco. A una frecuencia de muestreo de 40KHz, tOneFrame = 1024 / 40K segundos = 1 / 40s, por lo que Bin [1] idealmente estaría recolectando una señal de 40Hz.

• Si configuramos osamp (overSample) = 4, avanzamos a lo largo de nuestra señal de entrada en los pasos de 256. Por lo tanto, el primer análisis examina los bytes de cero a 1023, luego de 256 a 1279, etc. Tenga en cuenta que 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;

(EDITAR :) Ahora lo poco que no consigo:

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

Simplemente no puedo verlo claramente, a pesar de que parece estar mirando a la cara. ¿Podría alguien explicar este proceso desde cero, paso a paso?

Answer 1

El principio básico es muy simple. Si un componente dado coincide exactamente con una frecuencia bin, entonces su fase no cambiará de un FT al siguiente. Sin embargo, si la frecuencia no se corresponde exactamente con la frecuencia del contenedor, habrá un cambio de fase entre los FT sucesivos. El delta de frecuencia es justo:

delta_freq = delta_phase / delta_time

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

freq_est = bin_freq + delta_freq

Answer 2

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

Si observa un solo punto en una onda sinusoidal (frecuencia fija y amplitud fija) en el tiempo, la fase avanzará con el tiempo en una cantidad proporcional a la frecuencia. O puede hacer lo contrario: si mide cuánto cambia la fase de una sinusoide en cualquier unidad de tiempo, puede calcular la frecuencia de esa sinusoide.

Un vocoder de fase usa dos FFT para estimar la 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. Desde allí, tiene su estimación de frecuencia para ese contenedor de FFT (un contenedor de FFT es aproximadamente un filtro para aislar un componente sinusoidal u otra señal de banda estrecha que se ajuste dentro de ese contenedor).

Para que este método funcione, el espectro cerca del contenedor de FFT en uso debe ser bastante estacionario, por ejemplo, no cambiar la frecuencia, etc. Ese es el supuesto que requiere un vocoder de fase.

Answer 3

Finalmente he descubierto esto; Realmente tuve que derivarlo desde cero. Sabía que habría una forma sencilla de derivarlo, mi error (habitual) era intentar seguir la lógica de otras personas en lugar de usar mi propio sentido común.

Este puzzle tiene dos claves para desbloquearlo.

• La primera clave es entender cómo el muestreo excesivo introduce una rotación en la fase del contenedor.

• La segunda clave viene de Graph 3.3 & 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;
}

Answer 4

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

Answer 5

He implementado este algoritmo para Performous yo mismo. Cuando toma otra FFT en un desplazamiento de tiempo, espera que la fase cambie de acuerdo con el desplazamiento, es decir, dos FFT tomadas con 256 muestras separadas deben tener una diferencia de fase de 256 muestras para todas las frecuencias presentes en la señal (esto supone que las señales en sí mismas son estables, lo que es una buena suposición para períodos cortos como 256 muestras).

Ahora, los valores de fase reales que obtiene de FFT no están en las muestras sino en el ángulo de fase, por lo que serán diferentes dependiendo de la frecuencia. En el siguiente código, el valor de PhaseStep es el factor de conversión necesario por cada bin, es decir, para la frecuencia correspondiente al bin x, el cambio de fase será x * phaseStep. Para las frecuencias centrales de los contenedores, x sería un número entero (el número de bin) pero para las frecuencias reales detectadas 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 contenedor tiene la frecuencia central del contenedor y luego calcula el cambio de fase esperado para eso. Este cambio esperado se resta del cambio real, dejando el error. Se toma un resto (módulo 2 pi) (-pi al rango pi) y la frecuencia final se calcula con el centro de bin + 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

Observe que muchas bandejas adyacentes a menudo terminan corregidas a la misma frecuencia porque la corrección delta puede ser de hasta 0,5 * bandejas FFT_N / FFT_STEP de cualquier manera, por lo que la FFT_STEP más pequeña que use, más lejos serán posibles las correcciones (pero esto aumenta la potencia de procesamiento necesario, así como imprecisión debido a inexactitudes).

Espero que esto ayude :)

Answer 6

Tal vez esto ayude. Piense en las bandejas FFT como si se tratara de pequeños relojes o rotores, cada uno girando a la frecuencia de la bandeja. Para una señal estable, la siguiente posición (teórica) del rotor se puede predecir usando las matemáticas en el bit que no se obtiene. Contra esta posición "debería ser" (ideal), puede calcular varias cosas útiles: (1) la diferencia con la fase en un intervalo de una trama adyacente, que es utilizada por un vocodificador de fase para estimar mejor la frecuencia de intervalo, o (2) más generalmente la desviación de fase , que es un indicador positivo del inicio de una nota o algún otro evento en el audio.