# Extracting phase-shift and gain from a time series information

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.

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?

• because it's a packaged function in MATLAB, scipy, or almost any other language's signal processing library, I'd tend to take the FFT of each, and compare the phases of the peak. Just the DFT at the single frequency would be faster to run, but probably slower to write. – Neil_UK Nov 22 '17 at 16:26
• Wouldn't it be easier to find the zero crossing point and use that information to calculate the phase shift? – Harry Svensson Nov 22 '17 at 17:43
• @HarrySvensson, that method is easy to understand, but tends to give noisier results than using all the waveform data. – The Photon Nov 22 '17 at 17:46
• This problem has been heavily researched for applications in remote sensing (radar). There are whole books written on it. – The Photon Nov 22 '17 at 17:52
• @Neil_UK, that can be tricky to get right if there are not an integer number of signal periods in the sampling window. – The Photon Nov 22 '17 at 17:59

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:

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.

• Thanks I checked your code works fine but at some angles I get peculiar results for phase shifts. Your code phase1 = atan2(a,b); and phase = atan2(a,b); SO when I want to plot phase shift isnt it (phase - phase1)? At some frequencies I get jumps. But the plots follows the logical shift, Do you have an idea ? – user1999 Nov 23 '17 at 12:36
• Hmm, I recall once seeing a problem similar to what you've described caused by some offsets of 2*pi. Im on mobile so I can't give you a detailed explanation right now, but you can look up the matlab "unwrap" function in the meantime, it's a really sophisticated function for automatically correcting erroneous phase offsets that result from these types of algorithms sometimes. – phillipd94 Nov 23 '17 at 15:51
• diffinrad = phase1-phase; d= rad2deg(diffinrad ); if d>0 d= d-360; end – user1999 Nov 24 '17 at 11:09
• this solved the problem i dont know why – user1999 Nov 24 '17 at 11:09
• Oh, I see what the problem was then. At steady state, a phase shift of -60 degrees is the same as a phase shift of +300 degrees, so it's possible for your numbers to be off by 360 degrees (or any multiple thereof) and still represent the same sin wave. That's because a sin wave has a period of 360 degrees, and if I shift a sin wave left or right by 360 degrees it is still a sin wave. But, these extra 360 offsets can cause jumps like what you saw in your numbers. The code you added removed the offset and made all of your phase values in the range of [-360,0]. – phillipd94 Nov 24 '17 at 15:14

There are a number of methods to extract the gain & phase difference each with their own set of pro's and con's

Method1: Zero crossing & peak,rms

The method that almost always comes to mind when 1st met with this challenge is measure the time difference between zero crossing & then a simple $\Theta = 360 f \Delta T$. In practice determining the zero crossing, especially with acquired data, is extremely error prone primarily due to noise. Depending on the time-step of the acquisition, the actual zero crossing can be further

Peak detection to then determine the gain is equally prone to errors due to the same issue mentioned above.

Comparing the RMS values of both waveforms however is very robust for determining the magnitude

Method2: Fourier transform

It is usually tempting to determine the spectral content of the two signals to then work out the magnitude and phase difference. This is a valid method but care must be taken with respect to the data acquired, the sample rate & windowing. If you acquired multiple complete cycles then such methods are very accurate. However,

Method3 Lock-in

To clearly extract the phase and magnitude then two oscillators, in quadrature, are required.

$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))$ $V_{90} = Xcos(\omega t)Asin(\omega t +\phi) = \frac{XA}{2}(sin(\phi) + sin(2\omega t + \phi))$

Filter these signals to remove the twice carrier component

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

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

via trig:

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

## EXAMPLE:

%% 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);

As you can see, it requires a complete cycle to settle and this is due to the moving average filter requiring a fully cycle to accumulate enough samples. A filter is needed to suppress the $2\omega$ term \$ this could be a 1st order filter with a very low cutoff frequency BUT then the number of cycles needed to allow the filter to settle would be longer. Likewise it would suppress any AC component to the phase (something that could occur in other systems).