Electronic – Extracting phase-shift and gain from a time series information

seriestimevoltage

I have many files consisting input and output voltage signals for a linear filter circuit. For each file the input and output signal is at different constant frequency. All sampled for 2 sec duration with a sampling rate Fs. So basically by using a function generation each time I applied sinusoidal input signal Vin at a constant freq. to the input and I obtain a sinusoidal output. I want to plot the phase-shift versus frequency. Below is an example for 5Hz input in time-series.

enter image description here

By using the time domain data above how can I algorithmically extract the phase shift? What comes to my mind is to find the first maximum point for each sinusoid and the difference will reveal the phase shift. But I dont know how to implement this. And for the freq. vs amplitude I need to find the pk-pk values(because sinusoids have some offsets) and extract the ratio.

Is there a way to implement this?

Best Answer

This seems like a great chance to use the discrete fourier series. It is possible to avoid all the nasty transformations and only pull one term from the series, and correlate your data with only one desired frequency for comparison. This method essentially uses only one term of the discrete fourier series to correlate the timeseries data with only one discretized frequency. The math gets pretty nasty, and I'm not great with the mathscript on this site yet, so I'll provide you an external link to an eloquent explanation of the theory behind it. I'll warn you, it's very involved, and frankly unnecessary for simple applications.

https://www.math.ubc.ca/~feldman/m267/dft.pdf

I myself will provide some application examples with matlab code, easily ported to something like python if you don't have access.

What we need to do is take the time series data, and multiply each data point by the value of a cosine wave at that time point, and sum them all up. Then repeat with a sin wave. Yep, all that nasty math, and it boils down to something that simple. That will give us coefficients for the correlation between our time series and the sin and cosine waves of the desired frequency. Once we have these two numbers, use the standard magnitude and phase equations:

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

phase = atan2(a,b)

Here is some code that generates a sine wave, simulates an output from an arbitrary linear system, and compares their magnitude and phase. I set up my data to have 1000 data points and a timestep of .001, you will need to change these parameters to match your sample rates:

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)

code outputs:

Input and Output Sine Waves (arbitrary)

mag =

    1.0000


phase =

  -6.5239e-17


mag1 =

    0.9539


phase1 =

   -0.2984

As you can see, I generated a 5Hz (or 5*2*pi radian) sine wave, simulated an output not unlike the one you see in your time series, and performed the necessary correlations in for loops to compute the magnitude and phase. To get the magnitude and phase shift you only need to take the difference between the two magnitude and phase outputs, and boom, you've got it.

Let me know if I can edit my answer to clear anything up, I've crammed a lot of discrete-time systems analysis material into a very small space here so I'm happy to expand if necessary.

Related Topic