Extraer el cambio de fase y la ganancia de una información de series de tiempo

Esto parece ser una gran oportunidad para usar la discreta serie de Fourier. Es posible evitar todas las transformaciones desagradables y solo extraer un término de la serie, y correlacionar sus datos con una sola frecuencia deseada para la comparación. Este método utiliza esencialmente solo un término de la serie de Fourier discreta para correlacionar los datos de series temporales con una sola frecuencia discretizada. Las matemáticas se vuelven bastante desagradables, y todavía no estoy muy bien con el script de matemáticas en este sitio, así que les proporcionaré un enlace externo a una explicación elocuente de la teoría que la respalda. Te lo advierto, es muy complicado y francamente innecesario para aplicaciones simples.

enlace

Yo mismo proporcionaré algunos ejemplos de aplicaciones con código matlab, que se pueden trasladar fácilmente a algo como python si no tienes acceso.

Lo que debemos hacer es tomar los datos de la serie de tiempo, y multiplicar cada punto de datos por el valor de una onda de coseno en ese punto de tiempo, y sumarlos todos. Luego repite con una onda de pecado. Sí, todas esas matemáticas desagradables, y se reduce a algo tan simple. Eso nos dará los coeficientes para la correlación entre nuestras series de tiempo y las ondas de pecado y coseno de la frecuencia deseada. Una vez que tengamos estos dos números, use las ecuaciones estándar de magnitud y fase:

mag = sqrt (a ^ 2 + b ^ 2)

fase = atan2 (a, b)

Aquí hay un código que genera una onda sinusoidal, simula una salida de un sistema lineal arbitrario y compara su magnitud y fase. Configuré mis datos para tener 1000 puntos de datos y un paso de tiempo de .001, necesitarás cambiar estos parámetros para que coincidan con tus frecuencias de muestreo:

N=1000; %simulation parameters, adjust to your measurement frequency
T = 0.001;

freq=5*2*pi; %5 Hz

n = 0:1:(N-1);
t = n*T;


y=sin(freq*t);

G = tf(100,[1 100]); %arbitrary system for example

y1=lsim(G,y,t)'; %get a linear simulation output from the example system

plot(t,y1,t,y)



a_sum=0;
b_sum=0;
w = 5*2*pi; %declare relevant frequency for comparison
for (k=0:(N-1)) %iterate over all samples
    time = k*T;
    a_sum = a_sum + y(k+1)*cos(w*time);%perform fourier correlation with ONLY THE RELEVANT FREQUENCY
    b_sum = b_sum + y(k+1)*sin(w*time);
end
a = a_sum*2/N;
b = b_sum*2/N;
mag = sqrt(a^2+b^2)  %magnitude formula
phase = atan2(a,b)  %phase is arctangent of fourier coefficients


a_sum=0;  %now do it all again for the output wave
b_sum=0;
w = 5*2*pi;
for (k=0:(N-1))
    time = k*T;
    a_sum = a_sum + y1(k+1)*cos(w*time);
    b_sum = b_sum + y1(k+1)*sin(w*time);
end
a = a_sum*2/N;
b = b_sum*2/N;
mag1 = sqrt(a^2+b^2)
phase1 = atan2(a,b)

salidas de código:

mag =

    1.0000


phase =

  -6.5239e-17


mag1 =

    0.9539


phase1 =

   -0.2984

Como puede ver, generé una onda sinusoidal de 5Hz (o 5 * 2 * pi radian), simulé una salida no muy diferente a la que se ve en su serie de tiempo y realicé las correlaciones necesarias para que los bucles calculen la magnitud y la fase. Para obtener la magnitud y el cambio de fase, solo necesitas tomar la diferencia entre las dos salidas de magnitud y fase, y el auge, lo tienes.

Avíseme si puedo editar mi respuesta para aclarar cualquier cosa, he metido una gran cantidad de material de análisis de sistemas de tiempo discreto en un espacio muy pequeño aquí, así que me complace expandirlo si es necesario.

Hay varios métodos para extraer la ganancia & diferencia de fase, cada uno con su propio conjunto de pros y contras

Método 1: cruce por cero y amp; pico, rms

El método que casi siempre me viene a la mente cuando primero se enfrentó con este desafío es medir la diferencia de tiempo entre el cruce por cero y amp; luego un simple \ $ \ Theta = 360 f \ Delta T \ $. En la práctica, determinar el cruce por cero, especialmente con los datos adquiridos, es extremadamente propenso a errores principalmente debido al ruido. Dependiendo del período de tiempo de la adquisición, el cruce por cero real puede ser mayor

La detección de picos para luego determinar la ganancia es igualmente propensa a errores debido al mismo problema mencionado anteriormente.

Sin embargo, la comparación de los valores RMS de ambas formas de onda es muy robusta para determinar la magnitud

Método 2: Transformada de Fourier

Por lo general, es tentador determinar el contenido espectral de las dos señales para luego calcular la magnitud y la diferencia de fase. Este es un método válido, pero se debe tener cuidado con respecto a los datos adquiridos, la frecuencia de muestreo & ventanas Si adquirió varios ciclos completos, estos métodos son muy precisos. Sin embargo,

Método 3 Bloqueo

Para extraer claramente la fase y la magnitud, se requieren dos osciladores, en cuadratura.

\ $ V_ {sig} = Asin (\ omega t + \ phi) \ $

\ $ V_ {osc0} = Xsin (\ omega t) \ $

\ $ V_ {osc90} = Xcos (\ omega t) \ $

\ $ V_0 = Xsin (\ omega t) Asin (\ omega t + \ phi) = \ frac {XA} {2} (cos (\ phi) - cos (2 \ omega t + \ phi)) \ PS \ $ V_ {90} = Xcos (\ omega t) Asin (\ omega t + \ phi) = \ frac {XA} {2} (sin (\ phi) + sin (2 \ omega t + \ phi)) \ $

Filtre estas señales para eliminar el componente de portadora dos veces

\ $ V_ {0f} = \ frac {XA} {2} (cos (\ phi)) \ $

\ $ V_ {90f} = \ frac {XA} {2} (sin (\ phi)) \ $

a través de trig:

\ $ \ phi = atan (\ frac {V_ {90f}} {V_ {0f}}) \ $ \ $ A = \ frac {2} {X} \ sqrt {V_ {0f} ^ 2 + V_ {90f} ^ 2} \ $

EJEMPLO:

%% Generate the Stimulus and Signal 
F = 50;


t = linspace( 0, 5/50,10000);
stim = sin(2*pi*F*t + pi/10);  % stimulus. PHase shift for indication purposes and post-processing 
sig = lsim( tf([50*2*pi],[1 50*2*pi]),stim,t)'; % signal is the stimulus post 50Hz 1st order filter, gain should be -3dB & phase should be 45deg 

%%%%%%%%%%%%%%%%%%%%%%%%%
%% Post Processing 
%%%%%%%%%%%%%%%%%%%%%%%%%

osc1 = sin(2*pi*F*t );
osc2 = cos(2*pi*F*t );

% Calculate components w.r.t. stimulus
V0 = osc1.*stim;
V90 = osc2.*stim;
tmp = ceil(1/(F*t(2)));
B = (tmp)*ones(1,tmp); % moving average filter,over the fundemental period. Suppress the 2nd harmonic that is created

V0f = filter( B,1,V0);
V90f = filter(B,1, V90);

angle1 = atan2(V90f,V0f);
Amp1 = 2*sqrt(V0f.^2 + V90f.^2);

% Calculate components w.r.t. signal
V0 = osc1.*sig;
V90 = osc2.*sig;
tmp = ceil(1/(F*t(2)));
B = (tmp)*ones(1,tmp); % moving average filter,over the fundemental period. Suppress the 2nd harmonic that is created

V0f = filter( B,1,V0);
V90f = filter(B,1, V90);

angle2 = atan2(V90f,V0f);
Amp2 = 2*sqrt(V0f.^2 + V90f.^2);


% calculate gain & phase:  
gain = 20*log10(Amp2./Amp1);
phase = angle2 - angle1;

subplot(3,1,1);
plot(t,sig,t,stim);
grid on

subplot(3,1,2);
plot(t,gain);
grid on;

subplot(3,1,3);
plot(t,rad2deg(phase));
grid on;

Como puede ver, requiere un ciclo completo para establecerse y esto se debe a que el filtro de promedio móvil requiere un ciclo completo para acumular suficientes muestras. Se necesita un filtro para suprimir el término \ $ 2 \ omega \ $ término $ este podría ser un filtro de primer orden con una frecuencia de corte muy baja, PERO entonces la cantidad de ciclos necesarios para permitir que el filtro se asiente sería más larga. Asimismo, suprimiría cualquier componente de CA en la fase (algo que podría ocurrir en otros sistemas).