Obtenga el valor máximo de una señal si su frecuencia se encuentra entre dos centros de depósito

12

Suponga lo siguiente:

  • La frecuencia fundamental de una señal se ha estimado utilizando FFT y algunos métodos de estimación de frecuencia y se encuentra entre dos centros de depósito.
  • La frecuencia de muestreo es fija
  • El esfuerzo computacional no es un problema

Conociendo la frecuencia, ¿cuál es la forma más precisa de estimar el valor pico correspondiente de las señales fundamentales?

Una forma podría ser poner a cero la señal de tiempo para aumentar la resolución FFT de modo que el centro del depósito esté más cerca de la frecuencia estimada. En este escenario, un punto del que no estoy seguro es si puedo hacer cero-pad tanto como quiera o si hay algunos inconvenientes al hacerlo. Otro es el centro del contenedor que debería seleccionar después del relleno cero como el que estoy obteniendo el valor máximo (porque uno no puede alcanzar la frecuencia de interés exactamente, incluso después del relleno cero).

Sin embargo, también me pregunto si existe otro método que pueda ofrecer mejores resultados, por ejemplo, un estimador que utiliza los valores máximos de los dos centros de contenedores circundantes para estimar el valor máximo a la frecuencia de interés.

lR8n6i
fuente
2
cero relleno antes de FFT es unidireccional. Otra es aplicar una función de ventana que sea adecuada para sus neads. La ventana superior plana fue diseñada exactamente para este propósito. Por supuesto, si ya conoce la frecuencia exactamente y está interesado en una sola amplutida, probablemente haya formas más baratas de hacerlo que una FFT.
sellibitze
1
no se requiere relleno cero: la interpolación parabólica simple (con 3 puntos: imax-1, imax, imax + 1, donde imaxestá el pico FFT) le dará resultados precisos
Basj
Asegúrese de que la función de interpolación coincida con la función de ventana. Flat-top es trivial, de lo contrario, desea un par coincidente (por ejemplo, ventana rectangular + interpolación sinc, ventana gaussiana + interpolación gaussiana, etc.)
finnw
@CedronDawg esta pregunta y sus respuestas están relacionadas (pero no son las mismas) con su fórmula de frecuencia exacta. Puede ser que pueda encontrarlo interesante.
Fat32

Respuestas:

5

El primer algoritmo que me viene a la mente es el algoritmo de Goertzel . Ese algoritmo generalmente supone que la frecuencia de interés es un múltiplo entero de la frecuencia fundamental. Sin embargo, este documento aplica el algoritmo (generalizado) al caso que le interesa.


Otro problema es que el modelo de señal es incorrecto. Lo usa 2*%pi*(1:siglen)*(Fc/siglen). Debería usarse 2*%pi*(0:siglen-1)*(Fc/siglen)para que la fase salga correctamente.

También creo que hay un problema con la frecuencia Fc=21.3es muy baja. Las señales de baja frecuencia con valores reales tienden a mostrar sesgo cuando se trata de problemas de estimación de fase / frecuencia.

También probé una búsqueda de grilla gruesa para la estimación de fase, y da la misma respuesta que el algoritmo de Goertzel.

A continuación se muestra un gráfico que muestra el sesgo en ambas estimaciones (Goertzel: azul, grueso: rojo) para dos frecuencias diferentes: Fc=21.3(sólido) y Fc=210.3(discontinua). Como puede ver, el sesgo de la frecuencia más alta es mucho menor.

x2π

ingrese la descripción de la imagen aquí

Peter K.
fuente
Acabo de probar el código para el algoritmo Goerzel basado en el documento. Usando el valor de salida DTFT, el pico se puede obtener con mucha precisión. Sin embargo, hay un factor de escala de exactamente 1000. Entonces, si el pico original es 1,234, después de Goerzel será 1234. ¿Alguien sabe de dónde podría venir esto?
lR8n6i
Investigé un poco mientras tanto. Probablemente tiene que ver con la escala de amplitud: amplitud del dominio del tiempo de escala = coeficiente del dominio de frecuencia * 2 / N, donde N es la longitud de la señal. ¿Es correcta esta suposición?
lR8n6i
¡Hola! Acabo de descubrir que usando el algoritmo de Goertzel, la amplitud en el coeficiente complejo resultante es muy precisa, pero la fase es completamente incorrecta. ¿Alguien tiene una idea de dónde podría venir esto? Por "fase" me refiero al retraso de fase especificado en el fundamental de la señal original.
lR8n6i
1
sin(ω0t+ϕ)j2[ejϕδ~(ω+ω0+2πk)e+jϕδ~(ωω0+2πk)]π/2
4

Si está dispuesto a usar varios contenedores FFT vecinos, no solo 2, la interpolación Sinc en ventana entre los resultados del contenedor complejo puede producir una estimación muy precisa, dependiendo del ancho de la ventana.

La interpolación Sinc con ventana se encuentra comúnmente en muestreadores de audio de alta calidad, por lo que los documentos sobre ese tema tendrán fórmulas de interpolación adecuadas con análisis de errores.

hotpaw2
fuente
Gracias por el comentario. También intentaré este enfoque.
lR8n6i
4

sin(πx)(πx)

[1] JL Flanagan y RM Golden, "Phase vocoder", Bell Systems Technical Journal, vol. 45, págs. 1493–1509, 1966.

[2] K. Dressler, "Extracción sinusoidal utilizando una implementación eficiente de una FFT de resolución múltiple", en Proc. 9th Int. Conf. on Digital Audio Effects (DAFx-06), Montreal, Canadá, septiembre de 2006, págs. 247–252.

ederwander
fuente
¡Hola! Muchas gracias por todos sus comentarios. Extendí mi código (ver más abajo) para combinar el filtro Goertzel con la interpolación parabólica de picos para obtener la fase. Sin embargo, los resultados aún no son precisos (+ - 3-4deg). ¿Es esto lo más cercano posible o hay errores en la comprensión o la codificación?
lR8n6i
3

Tuve muchas dificultades con este problema exacto hace un par de años.

Publiqué esta pregunta:

/programming/4633203/extracting-precise-frequencies-from-fft-bins-using-phase-change-between-frames

Terminé haciendo los cálculos desde cero y publiqué una respuesta a mi propia pregunta.

Me sorprende que no haya podido encontrar ninguna exposición similar en Internet.

Publicaré la respuesta nuevamente aquí; tenga en cuenta que el código está diseñado para un escenario en el que estoy superponiendo mi ventana FFT por 4x.

π


Este rompecabezas tiene dos llaves para desbloquearlo.

Gráfico 3.3:

ingrese la descripción de la imagen aquí

Gráfico 3.4:

ingrese la descripción de la imagen aquí

Código:

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 bin1Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin1Freq;

    // 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;
}
Pi
fuente
Usted está interpolando la frecuencia, mientras que el OP conoce la frecuencia y quiere interpolar la amplitud.
finnw
2

Este código de Python le dará un resultado muy preciso (lo usé para muchas notas musicales y obtuve errores de menos del 0,01% de semitono) con interpolación parabólica (método utilizado con éxito por McAulay Quatieri, Serra, etc. en armónico + residual técnicas de separación)

import matplotlib.pyplot as plt
import numpy as np
from scipy.io.wavfile import read
from scipy.fftpack import fft, ifft
import math

(fs, x) = read('test.wav')
if (len(x.shape) == 2):    # if stereo we keep left channel only
 x = x[:,1]

n=x.size
freq = np.arange(n)*1.0/n*fs 
xfft = abs(fft(x))

imax=np.argmax(xfft)  
p=1.0/2*(xfft[imax-1]/xfft[imax]-xfft[imax+1]/xfft[imax])/(xfft[imax-1]/xfft[imax]-2+xfft[imax+1]/xfft[imax])   # parabolic interpolation 
print 'Frequence detectee avec interpolation parabolique :',(imax+p)*1.0/n*fs, 'Hz'
Basj
fuente
1
clear all
clc

for phase_orig = 0:pi/18:pi,

%% Specify and generate signal
Amp = 1;                     % Amplitude of signal
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
Fc = 21.3;                   % Hz
StopTime = 0.25;             % seconds
t = (0:dt:StopTime-dt)';     % seconds

siglen = length(t);
sig = Amp * 1.5 * sin(2*pi*(0:siglen-1)*(Fc/siglen) + phase_orig) + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 3) ...
  + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 5)+ 0.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 7) ...
  + 1.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 9)+ 1.4 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 11);

%% Estimate the peak value of the signals fundamental using Goertzel algorithm
peak = 0;
indvec = [Fc-1 Fc Fc+1];

% Check the input data
if ~isvector(sig) || isempty(sig)
  error('X must be a nonempty vector')
end

if ~isvector(indvec) || isempty(indvec)
  error('INDVEC must be a nonempty vector')
end
if ~isreal(indvec)
  error('INDVEC must contain real numbers')
end

% forcing x to be column
sig = reshape(sig,siglen,1);

% initialization
no_freq = length(indvec); %number of frequencies to compute
y = zeros(no_freq,1); %memory allocation for the output coefficients

% Computation via second-order system
% loop over the particular frequencies
for cnt_freq = 1:no_freq
  %for a single frequency:
  %a/ precompute the constants
  pik_term = 2*pi*(indvec(cnt_freq))/(siglen);
  cos_pik_term2 = cos(pik_term) * 2;
  cc = exp(-1i*pik_term); % complex constant
  %b/ state variables
  s0 = 0;
  s1 = 0;
  s2 = 0;
  %c/ 'main' loop
  for ind = 1:siglen-1 %number of iterations is (by one) less than the length of signal
    %new state
    s0 = sig(ind) + cos_pik_term2 * s1 - s2;  % (*)
    %shifting the state variables
    s2 = s1;
    s1 = s0;
  end
  %d/ final computations
  s0 = sig(siglen) + cos_pik_term2 * s1 - s2; %correspond to one extra performing of (*)
  y(cnt_freq) = s0 - s1*cc; %resultant complex coefficient

  %complex multiplication substituting the last iterationA
  %and correcting the phase for (potentially) non-integer valued
  %frequencies at the same time
  y(cnt_freq) = y(cnt_freq) * exp(-1i*pik_term*(siglen-1));
end

  % perfom amplitude scaling
  peak = abs(y(2)) * 2 / siglen

% perform parabolic interpolation to get the phase estimate
phase_orig=phase_orig*180/pi
ym1 = angle(unwrap(y(1)));
y0 = angle(unwrap(y(2)));
yp1 = angle(unwrap(y(3)));

p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1)); 
phase = y0 - 0.25*(ym1-yp1)*p;
phase_est = phase * 180/pi + 90;
phase_est = mod(phase_est+180,360)-180
end

Las frecuencias con las que está tratando (21.3Hz muestreado a 8kHz) son muy bajas. Debido a que estas son señales de valor real, exhibirán un sesgo en la estimación de fase para ** cualquier ** frecuencia.

Esta imagen muestra una gráfica del sesgo ( phase_est - phase_orig) para Fc = 210.3;(en rojo) versus el sesgo para Fc = 21.3;. Como puede ver, el desplazamiento es mucho más significativo para el 21.3caso.

Otra opción es reducir su frecuencia de muestreo. La curva verde muestra el sesgo para en Fs = 800lugar de 8000.

ingrese la descripción de la imagen aquí

lR8n6i
fuente
1
¡Gracias por la actualización! Mira mi trama; Todavía creo que cualquier estimador de fase va a tener un sesgo para esta baja frecuencia. Una forma de evitarlo es usar la frecuencia conocida (si se conoce) para corregir el sesgo de la estimación de fase a través de una tabla de búsqueda. Pero deberá tener cuidado: el sesgo cambiará con la frecuencia. Otra forma de hacerlo será reducir su frecuencia de muestreo.
Peter K.
1
¡Gracias a ti también! Sin embargo, si está usando Fs = 8000 Hz y Fc = 210 en lugar de 210.3, el sesgo se ve aún peor. ¿Alguna idea de dónde podría venir esto?
lR8n6i
1
Erk! Ni idea. Fwiw, el estimador Geortzel no tiene problemas: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Excavará un poco más. Mira este espacio.
Peter K.
1
La interpolación parabólica no está haciendo lo que crees que está haciendo. En particular, si reemplaza su cálculo de pcon p2 = (abs(y(3)) - abs(y(1)))/(2*(2*abs(y(2)) - abs(y(3)) - abs(y(1)))); phase2 = y0 - 0.25*(ym1-yp1)*p2;, obtendrá MUCHO mejores respuestas, incluso para Fc=210. No estoy del todo seguro de que la versión actual de ple dé algo sensato. La fórmula de interpolación es para la interpolación de la AMPLITUD de una parábola, pero pestá interpolando la fase que es simplemente ... extraña.
Peter K.
1
Todo eso está bien, EXCEPTO que la ubicación del pico ( p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))) será incorrecta algunas veces si está usando FASES en lugar de amplitudes. Esto se debe a que las fases pueden saltar alrededor del límite de +/- 180 grados. Todo lo que se necesita para arreglarlo para la fase es cambiar esa línea a mi p2cálculo anterior.
Peter K.