Measure phase between two discrete signals

Good morning,

I have two sampled signals, corresponding to two sines. These signals are sampled with a low sampling rate, i.e. I get few samples within the measurement. Furthermore they are real signals and I have a certain amount of noise compared to the amplitude of the signal.

I would like to get the amplitude, which I think is done well and the phase difference between them, which is the part where I am getting problems.(Because I know that the value it returns is not correct)

Can you give me some advice?

Thank you!

clc
clear all
close all

%% Reading CSV from CC
%Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 3);

% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ";";

% Specify column names and types
opts.VariableNames = ["Time", "AmpPert", "Vout"];
opts.VariableTypes = ["double", "double", "double"];

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";


% Import the data
ADC = readtable("D:\Doctorado\10. Pruebas\2. RespuestaFrec\3. Digital\divisor\SegundoIntento\160kHz.csv", opts);
%% Clear temporary variables
clear opts

%% Entry of characteristic data
fs = 663e3;                     % Sampling frequency (samples per second)
dt = 1/fs;                      % seconds per sample
AdcTime = ADC.Time*dt;          % Samples convertsion to time

OffsetVout = mean(ADC.Vout);    %Measurement of fignal offset
ADC.Vout = ADC.Vout-OffsetVout; %Removing DC
AVout = smoothdata(ADC.Vout,"gaussian",15);  %smoothing the signal 
% AVout=ADC.Vout;           

OffsetPert = mean(ADC.AmpPert); %Measurement of fignal offset
ADC.AmpPert = ADC.AmpPert-OffsetPert; %Removing DC
APert = smoothdata(ADC.AmpPert,"gaussian",15);    %smoothing the signal
% APert=ADC.AmpPert;

figure(1)
plot(AdcTime,AVout,'r')
hold on
plot(AdcTime,50*APert,'b')
legend('\DeltaVout','\DeltaPert')

%% Extraction of disturbance information
Y = fft(APert);                         %Calculation of the FFT
L =length(Y);                           

%Magnitude acquisition
P2 = abs(Y/L);                          %Calculation of the magnitude of harmonics
P1Pert = P2(1:L/2+1);                   %Only the positive spectrum
P1Pert(1:end-1) = 2*P1Pert(1:end-1);    %Add the part of the negative part - which is the reflection - therefore by 2
f = fs*(0:(L/2))/L;

%Phase adquisition
tol = 3e-3;
Y(abs(Y) < tol) = 0;                    %Remove harmonics minor to tolerance
thetaPert = angle(Y);                   
thetaPert = thetaPert(1:L/2+1);
thetaPert = abs(thetaPert*180/pi);

%% Extraction of Vout information
YVout = fft(AVout);                     %Calculation of the FFT

%Magnitude acquisition
P2Vout = abs(YVout/L);                  %Calculation of the magnitude of harmonics
P1Vout = P2Vout(1:L/2+1);               %Only the positive spectrum
P1Vout(1:end-1) = 2*P1Vout(1:end-1);    %Add the part of the negative part - which is the reflection - therefore by 2 

%Phase adquisition
tol = 0.3;
YVout(abs(YVout) < tol) = 0;            %Remove harmonics minor to tolerance
thetaVout = angle(YVout);                   
thetaVout = thetaVout(1:L/2+1);
thetaVout = abs(thetaVout*180/pi);

%% Gain and phase calculations
% Disturbance signal
AmplitudePert = max(P1Pert);
Index = find(P1Pert==AmplitudePert);
Frec = f(Index)

%Vout signal
AmplitudeVout = P1Vout(Index);

%Gain calculation
Gain =20*log10(AmplitudeVout/AmplitudePert)

%Phase difference calculation
fase = (thetaVout(Index) - thetaPert(Index))

figure(2)
loglog(f,P1Vout,'r');
hold on
grid on; grid minor
loglog(f,P1Pert,'b')

%Another method to determine phase
[pksPert,locsPert] = findpeaks(APert);  %Find peaks of disturbance
[pksVout,locsVout] = findpeaks(AVout);  %Find peaks of disturbance

figure(1)
plot(AdcTime(locsVout), AVout(locsVout), '+r')
plot(AdcTime(locsPert), 50*APert(locsPert), '+r')

%Not working well
if length(locsVout)>length(locsPert)
    difference =  locsVout(2:end)-locsPert;
elseif length(locsVout)==length(locsPert)
    difference =  locsVout-locsPert;
else
    difference = NaN;
end

difference = mean(difference);
differenceTime =  difference*dt;
phasePeaks = differenceTime*Frec*360;

result = ['La frecuencia es ',num2str(Frec),' La ganancia es  ',num2str(Gain),' .La fase es ',num2str(fase), ' o mediate picos ', num2str(FasePeaks)];
disp(result)

Answer 1

[Mathieu NOE]

hello

I tried this :

• averaged fft spectrum, to find a dominant peak at 160570 Hz

• band pass filtered both signals
• zero crossing detection to measure period and time lag between signals , and found an average phase shift of 54 degrees (or time lag delay_AVout_vs_APert = 9.4084e-07 second)

• NB that the inverse of the signal period found above matches the fft peak at 160570 Hz (fortunately !)

Code :

clc
clear all
close all
%% Reading CSV from CC
%Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 3);
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ";";
% Specify column names and types
opts.VariableNames = ["Time", "AmpPert", "Vout"];
opts.VariableTypes = ["double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
ADC = readtable("160kHz.csv", opts);
%% Clear temporary variables
clear opts
%% Entry of characteristic data
fs = 663e3;                     % Sampling frequency (samples per second)
dt = 1/fs;                      % seconds per sample
AdcTime = ADC.Time*dt;          % Samples convertsion to time
OffsetVout = mean(ADC.Vout);    %Measurement of fignal offset
ADC.Vout = ADC.Vout-OffsetVout; %Removing DC
% AVout = smoothdata(ADC.Vout,"gaussian",15);  %smoothing the signal 
% AVout=ADC.Vout;           
OffsetPert = mean(ADC.AmpPert); %Measurement of fignal offset
ADC.AmpPert = ADC.AmpPert-OffsetPert; %Removing DC
% APert = smoothdata(ADC.AmpPert,"gaussian",15);    %smoothing the signal
% APert=ADC.AmpPert;

%% band pass  filter section %%%%%%
%bandpass filter
fc = 160570; % center frequency
a = 0.1; % relative bandwith
f_low = fc/(1+a);
f_high = fc*(1+a);
N_bpf = 4;           % filter order
[b,a] = butter(N_bpf,2/fs*[f_low f_high]);
AVout = filtfilt(b,a,ADC.Vout);
APert = filtfilt(b,a,ADC.AmpPert);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.75; % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display  : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_AVout] = myfft_peak(ADC.Vout,fs,NFFT,OVERLAP);
[freq, spectrum_APert] = myfft_peak(ADC.AmpPert,fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

figure(1),
plot(freq,20*log10(spectrum_AVout),freq,20*log10(spectrum_APert));grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');

%% find signals period and delay

threshold = 0; % 

t0_pos_AVout = find_zc(AdcTime,AVout,threshold);
t0_pos_APert = find_zc(AdcTime,APert,threshold);

period_AVout = mean(diff(t0_pos_AVout)); % delta time of crossing points
freq_AVout = (1./period_AVout) % signal mean frequency
period_APert = mean(diff(t0_pos_APert)); % delta time of crossing points
freq_APert = (1./period_APert) % signal mean frequency

% remove first sample (and make sure we have both arrays same size)
t0_pos_APert(t0_pos_APert<0.5*period_APert) = [];
t0_pos_AVout(t0_pos_AVout<0.5*period_AVout) = [];

delay_AVout_vs_APert = mean(t0_pos_AVout-t0_pos_APert); % time lag in seconds
delta_phase = 360*delay_AVout_vs_APert/period_AVout % in degrees

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
plot(AdcTime,AVout,'r',t0_pos_AVout,threshold*ones(size(t0_pos_AVout)),'dk','linewidth',2,'markersize',12)
hold on
plot(AdcTime,50*APert,'b',t0_pos_APert,threshold*ones(size(t0_pos_APert)),'dk','linewidth',2,'markersize',12)
legend('\DeltaVout','\DeltaPert')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[SERGIO TOLEDANO DUS]

Thank you very much!

I get results quite consistent with a lot of other data I had. So thank you very much!

The only question I have is: when applying the filter and the interpolation, do we introduce some kind of change in the phase?

Thank you!!

[Mathieu NOE]

hello again

glad it helps !

no, the filtering stuff does not modify the delta phase between signals

I have used filtfilt (and not filter) so as it works both forward and backward in time it does not create any phase distorsion

even if I had used filter (that works only in forward time) , the phase shift would have been the same for both signals as we are dealing with the same frequency peak, so once you do the delta, that shift has no more impact on the result

and the interpolation inside the zero crossing function does not either create any kind of phase distorsion