¿Existe un algoritmo para calcular la fase de una sola frecuencia?

Si tiene una función $f(t)=A \cdot \sin(\omega t+\phi)$ , y una onda referencia sin ( ω x $\sin(\omega x)$ ¿cuál sería un algoritmo rápido para calcular $\phi$ ?

Estaba mirando el algoritmo de Goertzel , pero no parece tratar con la fase.

Use un DFT a la frecuencia específica. Luego calcule la amplitud y la fase de las partes real / imag. Le proporciona la fase referenciada al inicio del tiempo de muestreo.

En una FFT 'normal' (o una DFT calculada para todos los N armónicos), normalmente se calcula la frecuencia con f = k * (sample_rate) / N, donde k es un número entero. Aunque pueda parecer sacrílego (especialmente para los miembros de la Iglesia del Entero Entero), en realidad puede usar valores no enteros de k cuando hace un solo DFT.

Por ejemplo, suponga que ha generado (u obtenido) N = 256 puntos de una onda sinusoidal de 27 Hz. (digamos, sample_rate = 200). Sus frecuencias 'normales' para una FFT de 256 puntos (o DFT de punto N) corresponderían a: f = k * (frecuencia de muestreo) / N = k * (200) / 256, donde k es un número entero. Pero una 'k' no entera de 34.56 correspondería a una frecuencia de 27 Hz., Utilizando los parámetros enumerados anteriormente. Es como crear un 'contenedor' DFT que esté exactamente centrado en la frecuencia de interés (27 Hz.). Algunos códigos C ++ (compilador DevC ++) pueden tener el siguiente aspecto:

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer
double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {
    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;
// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

(PD: Espero que lo anterior se traduzca bien en stackoverflow, algo de lo que podría envolver)

El resultado de lo anterior es una fase de -twopi / 4, como se muestra en los puntos reales generados (y el amplificador se duplica para reflejar la frecuencia pos / neg).

Algunas cosas a tener en cuenta: utilizo el coseno para generar la forma de onda de prueba e interpretar los resultados, debe tener cuidado con eso: la fase se refiere al tiempo = 0, que es cuando comenzó a muestrear (es decir, cuando recolectó r [0] ), y el coseno es la interpretación correcta).

El código anterior no es elegante ni eficiente (por ejemplo: use tablas de búsqueda para los valores sin / cos, etc.).

Sus resultados serán más precisos a medida que use N más grande, y hay un pequeño error debido al hecho de que la frecuencia de muestreo y N anteriores no son múltiplos entre sí.

Por supuesto, si desea cambiar su frecuencia de muestreo, N o f, tendría que cambiar el código y el valor de k. Puede colocar un contenedor DFT en cualquier lugar de la línea de frecuencia continua, solo asegúrese de estar usando un valor de k que corresponda a la frecuencia de interés.

El problema se puede formular como problema de mínimos cuadrados (no lineal):

$F(\phi)$$\phi$

La derivada es muy simple:

$F''(\phi)$

Obviamente, la función objetivo anterior tiene mínimos mínimos debido a la periodicidad, por lo tanto, se puede agregar algún término de penalización para discriminar otros mínimos (por ejemplo, agregar a la ecuación del modelo). Pero creo que la optimización solo convergerá a los mínimos más cercana y se puede actualizar el resultado de restar . 2 π k , k ∈ $\phi^{2}$$2\pi k, k\in N$

Hay varias formulaciones diferentes del algoritmo de Goertzel. Las que proporcionan 2 variables de estado (ortogonales o cercanas a), o una variable de estado compleja, como posibles salidas, a menudo se pueden usar para calcular o estimar la fase con referencia a algún punto de la ventana de Goertzel, como el centro. Los que proporcionan una sola salida escalar por lo general no pueden.

También necesitará saber dónde está su ventana de Goertzel en relación con su eje de tiempo.

Si su señal no es exactamente un número entero periódico en su ventana de Goertzel, la estimación de fase alrededor de un punto de referencia en el medio de la ventana puede ser más precisa que hacer referencia a la fase al principio o al final.

Una FFT completa es exagerada si conoce la frecuencia de su señal. Además, un Goertzel puede sintonizarse a una frecuencia no periódica en la longitud de la FFT, mientras que una FFT necesitará interpolación adicional o relleno de cero para frecuencias no periódicas en la ventana.

Un Goertzel complejo es equivalente a 1 contenedor de un DFT que utiliza una recurrencia para los vectores de coseno y seno seno o factores twiddle FFT.

Si sus señales no tienen ruido, puede identificar cruces por cero en ambos y determinar la frecuencia y la fase relativa.

Eso depende de cuál sea su definición de "rápido", qué tan precisa desea su estimación, si desea o la fase relativa a sus muestreos, y cuánto ruido hay en su función y onda sinusoidal de referencia.$\phi$

Una forma de hacerlo es tomar la FFT de y mirar el contenedor más cercano a . $f(t)$$\omega$ Sin embargo, esto dependerá de que esté cerca de la frecuencia central del contenedor.$\omega$

Entonces:

• ¿Qué quieres decir con "rápido"?
• ¿Qué tan precisa necesita la estimación?
• ¿Desea (fase relativa a la referencia) o fase relativa al inicio del muestreo? ¿Importa?$\phi$
• ¿Cuál es el nivel de ruido en cada señal?

PD: Supongo que quisiste decir , en lugar de .f ( t ) = A sin ( ω x + ϕ$f(t)=A\sin(ωt+ϕ)$$f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$

Punto de inicio:
1) multiplique su señal y la onda sin de referencia: = A⋅sin (ωt + ϕ) ⋅sin (ωt) = 0.5⋅A⋅ (cos (ϕ) - cos (2⋅ωt + ϕ) ) 2) encuentre integral en el período : 3) puede calcular :
T = π / ω I ( ϕ ) = ∫ T 0 F ( t ) d t = 0.5 ⋅ A ⋅ c o s ( ϕ ) ⋅ T ϕ c o s ( ϕ ) = I ( t ) / ( 0.5 ⋅ A ⋅ T $F(t)$
$T= \pi /\omega$
$I(\phi) = \int_0^T F(t)dt\ = 0.5⋅A⋅cos(ϕ) \cdot T$
$\phi$
$cos(\phi) = I(t)/(0.5 \cdot A \cdot T)$

Piensa en:
¿cómo medir A?
¿Cómo determinar en intervalo? (piense en " onda cos de referencia ")0 .. ( 2 ⋅ π $\phi$$0..(2 \cdot \pi)$

Para una señal discreta, cambie la integral a suma y elija cuidadosamente T!

También podría hacer esto (en notación numpy):

np.arctan( (signal*cos).sum() / (signal*sin).sum() ))

donde la señal es su señal de fase desplazada, cos y sin son las señales de referencia, y usted genera una aproximación de una integral durante un cierto tiempo mediante la suma de los dos productos.

Esta es una mejora en la sugerencia de @Kevin McGee de usar un DFT de frecuencia única con un índice de bin fraccionario. El algoritmo de Kevin no arroja grandes resultados: mientras que en medias bandejas y bandejas enteras es muy preciso, también cerca de las mitades y mitades también es bastante bueno, pero de lo contrario el error puede estar dentro del 5%, lo que probablemente no sea aceptable para la mayoría de las tareas .

Sugiero mejorar el algoritmo de Kevin ajustando , es decir, la longitud de la ventana DFT para que se acerque lo más posible a un todo. Esto funciona ya que a diferencia de FFT, DFT no requiere que sea ​​una potencia de 2.$N$$k$$N$

El siguiente código está en Swift, pero debe ser intuitivamente claro:

let f = 27.0 // frequency of the sinusoid we are going to generate
let S = 200.0 // sampling rate
let Nmax = 512 // max DFT window length
let twopi = 2 * Double.pi

// First, calculate k for Nmax, and then round it
var k = round(f * Double(Nmax) / S)

// The magic part: recalculate N to make k as close to whole as possible
// We also need to recalculate k once again due to rounding of N. This is important.
let N = Int(k * S / f)
k = f * Double(N) / S

// Generate the sinusoid
var r: [Double] = []
for i in 0..<N {
    let t = Double(i) / S
    r.append(sin(twopi * f * t))
}

// Compute single-frequency DFT
var R = 0.0, I = 0.0
let twopikn = twopi * k / Double(N)
for i in 0..<N {
    let x = Double(i) * twopikn
    R += r[i] * cos(x)
    I += r[i] * sin(x)
}
R /= Double(N)
I /= Double(N)

let amp = 2 * sqrt(R * R + I * I)
let phase = atan2(I, R) / twopi

print(String(format: "k = %.2f    R = %.8f    I = %.8f    A = %.8f    φ/2π = %.8f", k, R, I, amp, phase))