Formas de estimar el cambio de fase

La forma fácil de seguir pero con uso intensivo de la CPU

Si quieres algo que sea obvio, funciona bien y es fácil de entender, y no te importa lanzarle muchos ciclos de CPU (parece que no, si estás usando python), entonces podrías simplemente ajuste las ondas sinusoidales a las dos señales y lea la eliminación gradual de la función de ajuste.

Probablemente querrá usar scipy.optimise.curve_fit() para ajustar algo como \ $ V = A \ times \ sin \ left (2 \ pi f t + \ phi \ right) \ $ con \ $ \ phi \ $ como Solo parámetro libre. Esto podría converger en varios valores diferentes de \ $ \ phi \ $, así que tómelo modulo \ $ 2 \ pi \ $, luego haga lo mismo para la corriente y reste los dos valores de \ $ \ phi \ $, teniendo cuidado de agregar o reste \ $ 2 \ pi \ $ para obtenerlo en el rango significativo.

Esto requerirá mucho, mucho más uso de la CPU que el filtrado, o la detección homodina en tiempo real, pero es fácil de seguir, solo tiene unas pocas líneas y los ciclos de la CPU son baratos.

Si no puedes pagar esos ciclos de CPU

A continuación, puede intentar generar dos ondas sinusoidales de cuateratura en el software y luego usarlas para realizar la demodulación de IQ de las dos señales. Luego calcula las fases y resta. Esto no será tan fácil de seguir, pero podría hacerse muy rápido (incluso en python) y será muy, muy preciso.

Tome una parte de los datos que es un número entero de ciclos (asumiendo que su cuadrícula de 50Hz es agradable y estable, solo tome 1s de datos), y multiplíquelo punto por punto con cada una de las ondas sinusoidales generadas. Luego, integre cada uno usando numpy.trapz() para obtener dos escalares. Tome estos como los argumentos de numpy.arctan2() , y voila, tiene la fase de esa señal en relación con la onda sinusoidal generada. Haga lo mismo con la segunda señal y reste las dos para obtener la diferencia. Como se indicó anteriormente, agregue o reste \ $ 2 \ pi \ $ para obtenerlo en el rango significativo.

No lo menciona, pero según sus variables parece que está intentando medir el cambio de fase entre la corriente y el voltaje en una línea de 50 Hz.

Para realizar las mediciones de fase adecuadas, es esencial que la señal esté limpia y usted mencionó específicamente que no lo está. Para obtener una señal limpia, puede hacer un filtro de paso bajo de forma analógica antes del ADC o publicar digitalmente el ADC. Naturalmente, cualquier cambio de fase desde el filtrado debe ser igual para ambas señales o el error de fase introducido debe ser predecible.

Dos medios comunes para implementar un filtro de paso bajo digital son los algoritmos IIR y FIR. Afortunadamente, la biblioteca scipy los tiene incorporados. Por ejemplo, mire las funciones scipy ifilter. Sin embargo, tenga en cuenta que sus datos se consideran datos de tiempo no uniformes, por lo que primero deberá volver a muestrear sus datos utilizando una spline u otro algoritmo de interpolación apropiado.

Una vez que las señales han pasado por el filtro LP, la técnica más común para encontrar la diferencia de fase es encontrar el pico de cada señal, calcular la diferencia de tiempo entre picos y convertirla en grados o radianes.

Voy a dar dos opciones. No sé cuál dará mejores resultados.

1. Interpola tus datos en una cuadrícula de tiempo regular. Luego use todos los diversos métodos conocidos para la estimación de fase (cruces por cero, correlación cruzada, ...) que usted y otras respuestas mencionaron.

2. Construya una nueva matriz de seno y coseno en su cuadrícula de tiempo. Estimar

$$ A = \ int_0 ^ T f (t) \ sin \ left (2 \ pi f t \ right) {\ rm d} t $$ y $$ B = \ int_0 ^ T f (t) \ cos \ left (2 \ pi f t \ right) {\ rm d} t $$

para cada una de tus dos señales, usa trapz como lo hiciste para estimar la potencia promedio.

Ahora la fase de cada señal en relación con el coseno de referencia se puede estimar como \ $ \ tan ^ {- 1} \ left (\ frac {A} {B} \ right) \ $, o en Python numpy.arctan2(B, A) . Y la diferencia de fase entre las dos señales se puede obtener mediante subracción.

Advertencia: existe la posibilidad de que haya cometido un error en lo anterior que dé un error de \ $ \ pi / 2 \ $ en las estimaciones de fase de las señales individuales. También tenga cuidado de la posible \ $ 2 \ pi \ $ envolvente en la resta final. Verifique que los resultados sean razonables antes de usarlos.

Una forma es utilizar la detección sensible a la fase. A 50 Hz con 1000 muestras, está obteniendo 20 muestras por ciclo. Use bloques de 20 muestras y reste la suma de las segundas 10 muestras de la suma de las primeras 10 muestras. Esto le dará un valor relacionado con la relación de fase entre el límite del bloque y la señal mediante la función seno (tendrá que escalar la amplitud; puede usar un promedio del voltaje P-P). Estará en un pico positivo cuando esté exactamente en fase y un pico negativo cuando esté a 180. El enfoque elimina cualquier componente DC. Luego puede promediar este valor a lo largo de muchos ciclos para mitigar su ruido. Si mide ambas señales con el mismo tiempo de bloque, obtiene la relación de fase entre cada señal y su sistema de muestreo, lo que le da una respuesta. Alternativamente, puede crear un bucle que mueva el límite del bloque (cambiando las muestras) y mantenga la salida en cero, pero si toma este enfoque, su resolución es de 1/20 de un ciclo completo, que es de 18 grados (asumiendo que los tiempos de muestra son fijado a 1 kHz).

En el mundo del procesamiento de señales, la diferencia de fase se suele estimar con la correlación cruzada de dos señales; llamémoslos a (t) y b (t). Dado que a (t) y b (t) tienen formas de onda con formas similares, la correlación cruzada debería ser un enfoque confiable.

Aquí hay una secuencia de comandos de Python (parte de la cual se inspiró en esta SO respuesta ) que usa correlate que se encuentra en la biblioteca numpy para calcular el cambio de fase entre dos seno , a (t) y b (t ). Debería poder extrapolar el contenido de este script para su propia aplicación:

import numpy as np

# Create the time axis (seconds)
num_samples = 10001
samples_per_second = 1000
freq_Hz = 50.0
t = np.linspace(0.0, ((num_samples - 1) / samples_per_second), num_samples)
# Create a sine wave, a(t), with a frequency of 1 Hz
a = np.sin((2.0 * np.pi) * freq_Hz * t)
# Create b(t), a (pi / 2.0) phase-shifted replica of a(t)
b_shift = (np.pi / 2.0)
b = np.sin((2.0 * np.pi) * freq_Hz * t + b_shift)

# Cross-correlate the signals, a(t) & b(t)
ab_corr = np.correlate(a, b, "full")
dt = np.linspace(-t[-1], t[-1], (2 * num_samples) - 1)
# Calculate time & phase shifts
t_shift_alt = (1.0 / samples_per_second) * ab_corr.argmax() - t[-1]
t_shift = dt[ab_corr.argmax()]
# Limit phase_shift to [-pi, pi]
phase_shift = ((2.0 * np.pi) * ((t_shift / (1.0 / freq_Hz)) % 1.0)) - np.pi

manual_t_shift = (b_shift / (2.0 * np.pi)) / freq_Hz

# Print out applied & calculated shifts
print ("Manual time shift: {}".format(manual_t_shift))            --> 0.005
print ("Alternate calculated time shift: {}".format(t_shift_alt)) --> 0.005000000000000782
print ("Calculated time shift: {}".format(t_shift))               --> 0.005000000000000782                         
print ("Manual phase shift: {}".format(b_shift))                  --> 1.5707963267948966, b(t) relative to a(t)                            
print ("Calculated phase shift: {}".format(phase_shift))          --> -1.570796326794651, a(t) relative to b(t)

Cuando phase_shift es negativo, entonces sabe que la primera entrada de señal en correlate() se está retrasando con respecto a la segunda entrada de señal. En nuestro caso, esto significa que a (t) se está quedando atrás de b (t) por (pi / 2). Cuando phase_shift es positivo, entonces usted sabe que la primera entrada de señal está adelantando a la segunda entrada de señal.

Como puede ver, la correlación cruzada de dos señales se puede usar simplemente para detectar la diferencia de fase entre dos señales. El vector de error entre el real y amp; los cambios de tiempo / fase calculados permanecerán bajos siempre que sus señales estén relacionadas linealmente, muestreadas a una tasa por encima de su límite de Nyquist y tengan una SNR decente.

Dato curioso, aún puede utilizar la correlación cruzada para esta aplicación incluso si las señales son periódicas. Tenga cuidado al tomar la correlación cruzada entre dos señales que están relacionadas de manera no lineal.

Si su único problema con la medición es la ambigüedad, puede resolverlo volviendo a los datos, usando algún método para armonizar los cruces por cero demasiado cerca uno del otro (tome su tiempo de llegada promedio, tal vez), y luego calcule El ángulo utilizando el método de cruce por cero. El resultado que desee será cualquiera de los dos cambios de fase posibles que esté más cerca del resultado del cruce por cero.

Mientras estoy en el tema, un método para encontrar el ángulo de fase que nunca he intentado, pero debería funcionar, lo que también da resultados ambiguos, es medir el voltaje de cada onda sinusoidal y luego medir la diferencia de voltaje entre ellos. Eso te daría los 3 lados de un triángulo y puedes encontrar los ángulos con la ley de los senos / ley de los cosenos.