I have frequency sweep input signal x, and output signal y and want to know phase angle difference between them throughout the timeseries data.

In the attachment, first column is time, second is input signal x, and third is output signal y. The frequency is linearly varying in signal x from 0 to 6 Hz, and corresponding signal y may have different frequencies. I want to know the phase angle between both signals.

[Mathieu NOE]

hello

your x signal is clean but not your y signal, this is quite obvious when we look either on the raw (time domain signal) or on the spectrogram

see code and display below

now my question is do you want a phase lag info at each sinewave cycle or an averaged value for the entire sweep ?

x signal spectrogram :

it's clean for sure

y signal spectrogram :

harmonics + spurious signal between 1 and 3 Hz + noise

we can also notice that the system does not show any significant output between t = 0 and 100 s corresponding to x input

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('Phase_info.csv'); 
data = readmatrix(filename);
time = data(:,1);
signal = data(:,2:3); % select channel 2 or 3 or both
dt = mean(diff(time));
Fs = 1/dt; % sampling rate 

[samples,channels] = size(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 256*2;    % 
OVERLAP = 0.95;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 4;
if decim>1
    Fs = Fs/decim;
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)'*1/Fs;

%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);

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

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end

    % saturation of the dB range : 
    % saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

    % plots spectrogram
    figure(2+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
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

[Mathieu NOE]

and if you try to use tfestimate , the result is very poor : the coherence very low (way below 1) meaning the output signal is too much corrupted , or there is no much of the input signal present in the output signal

% 
filename = ('Phase_info.csv'); 
data = readmatrix(filename);
t = data(:,1);
x = data(:,2); %  
y = data(:,3);  % 
dt = mean(diff(t));
Fs = 1/dt; % sampling rate

figure(1),
plot(t,x,'r',t,y,'b');

NFFT = 512;
NOVERLAP = round(0.95*NFFT);  % 75 percent overlap

%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,Fs);

%% alternative with supplied function 
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,Fs,hanning(NFFT),NOVERLAP); 
% Txy = transfer function (complex), Cxy = coherence, F = freq vector

% Bode plots
% focus on the 0 -10 Hz range
ind = F<10;
figure(2),
subplot(3,1,1),plot(F(ind),20*log10(abs(Txy(ind))));
ylabel('Mag (dB)');
subplot(3,1,2),plot(F(ind),180/pi*(angle(Txy(ind))));
ylabel('Phase (°)');
subplot(3,1,3),plot(F(ind),Cxy(ind));
xlabel('Frequency (Hz)');
ylabel('Coh');


%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function  and Coherence Estimate 
    
% compute PSD and CSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x , y if length is less than the window length
    x(nwind)=0;
    y(nwind)=0;  
    n=nwind;
end
x = x(:);		% Make sure x is a column vector
y = y(:);		% Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;

Pxx = zeros(nfft,1); 
Pyy = zeros(nfft,1); 
Pxy = zeros(nfft,1); 

for i=1:k
    xw = window.*x(index);
    yw = window.*y(index);
    index = index + (nwind - noverlap);
    
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft);
    Xx2 = abs(Xx).^2;
    Yy2 = abs(Yy).^2;
    Xy2 = Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pyy = Pyy + Yy2;
    Pxy = Pxy + Xy2;
    
end

% Select first half
if ~any(any(imag([x y])~=0))   % if x and y are not complex
    if rem(nfft,2)    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end

    Pxx = Pxx(select);
    Pyy = Pyy(select);
    Pxy = Pxy(select);
    
else
    select = 1:nfft;
end
Txy = Pxy ./ Pxx;                   % transfer function estimate 
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy);    % coherence function estimate 

f = (select - 1)'*Fs/nfft;

end

Answer 1

[Mathieu NOE]

Hello again

I tried to extract the amplitude and phase of your x and y signals by "order extraction" , which here is simply a method where you apply the same logic as for a DFT computation, but instead of doing a true integral over N samples , you are simply multiplying the signals by cos and sin of the instantaneous frequency

you can see that the order extraction of the x (input) signal gives the correct amplitude (0.1) and 0 phase (which is normal as x is used to generate the cos signal from where we derive the sin signal).

There is some transient behaviou in the first 10 seconds because of the low pass filter transient behaviour for the real / imaginary part extraction of your signal , so don't pay too much attention at what happens in the first 20 seconds here.

now if you apply the same process to the y signal there is no constant phase , it rolls permanently as if y was not recorded synchronously with x .

even in the spectrogram plots you see from my earlier comments, there is like a small time shift between the x ramp up line and the corresponding y line in the second plot. can you confirm that when you record the data, x and y are recorded synchronously ?

code :

filename = ('Phase_info.csv'); 
data = readmatrix(filename);
t = data(:,1);
samples = numel(t);
dt = mean(diff(t));
Fs = 1/dt; % sampling rate 

x = data(:,2); %  
y = data(:,3);  % 


% plot x signal frequency vs time
tc1 = find_zc(t,x,max(x)/4); % find zero crossing points to determine period of signal (at each cycle)
p1 = (gradient(tc1)); %  period of signal x

% p1 = (diff(tc1)); %  period of signal x
% p1 = [p1(1);p1];
f1 = 1./p1;  % frequency of signal x (in Hz)
f = interp1(tc1,f1,t); % do some interpolation along t vector

%remove data where f is NaN
ind = isnan(f);
t(ind) = [];  % just because I want get rid of nans here
x(ind) = [];  % just because I want get rid of nans here
y(ind) = [];  % just because I want get rid of nans here
f(ind) = [];  % just because I want get rid of nans here
samples = numel(t);

figure(1),
plot(t,f);

% generate cos and sin signals for DFT computation
% % method #1
% ang = cumsum(2*pi*f)*dt;
% xr = cos(ang);
% xi = sin(ang);

% method #2
xr = x./max(abs(x)); % cos signal of amplitude 1 
xi = gradient(xr); % sin signal of non constant amplitude must be divided by instantaneaous pulsation / frequency
xi = smoothdata(xi,'gaussian',5); % a bit of smoothing to reduce noise created by signal differentiation
xi = xi./f; % amplitude correction :amplitude must be divided by instantaneaous pulsation / frequency
xi = xi./max(abs(xi(round(samples/4):end))); % sin signal of amplitude 1 

% figure(2),
% plot(t,xr,'r',t,xi,'b');

% DFT coefficients for signal : x
Mr = 2*x.*xr; % factor 2 alike in FFT : X = 2/n*sum(x*e(-jwt))
Mi = 2*x.*xi;
% low pass filter real / imag amplitudes (DFT coefficients)
[B,A] = butter(2,2/Fs*0.1);
Mr = filter(B,A,Mr);
Mi = filter(B,A,Mi);
x_amplitude = sqrt(Mr.^2 + Mi.^2);
x_phase = atan2(Mi,Mr); % in rads

figure(3),
subplot(2,1,1),plot(t,x_amplitude);
title('X signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
subplot(2,1,2),plot(t,x_phase);
xlabel('Time (s)');
ylabel('Phase (rad)');

% DFT coefficients for signal : y
Mr = 2*y.*xr; % factor 2 alike in FFT : X = 2/n*sum(x*e(-jwt))
Mi = 2*y.*xi;
% low pass filter real / imag amplitudes (DFT coefficients)
Mr = filter(B,A,Mr);
Mi = filter(B,A,Mi);
y_amplitude = sqrt(Mr.^2 + Mi.^2);
y_phase = atan2(Mi,Mr); % in rads

figure(4),
subplot(2,1,1),plot(t,y_amplitude);
title('Y signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
subplot(2,1,2),plot(t,y_phase);
xlabel('Time (s)');
ylabel('Phase (rad)');


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

[Mathieu NOE]

just to prove that the code is correct, let's assume y is now obtained by low pass fltering x

as we can compute the bode plot of that given filter , we can overlay both results

fortunately, both results are in good match

code :

filename = ('Phase_info.csv'); 
data = readmatrix(filename);
t = data(:,1);
samples = numel(t);
dt = mean(diff(t));
Fs = 1/dt; % sampling rate 

x = data(:,2); %  
% y = data(:,3);  % 
% your y data now replaced by a low pass filtered version of x 
[B0,A0] = butter(3,2/Fs*3);
y = filter(B0,A0,x);



% plot x signal frequency vs time
tc1 = find_zc(t,x,max(x)/4); % find zero crossing points to determine period of signal (at each cycle)
p1 = (gradient(tc1)); %  period of signal x

% p1 = (diff(tc1)); %  period of signal x
% p1 = [p1(1);p1];
f1 = 1./p1;  % frequency of signal x (in Hz)
f = interp1(tc1,f1,t); % do some interpolation along t vector

%remove data where f is NaN
ind = isnan(f);
t(ind) = [];  % just because I want get rid of nans here
x(ind) = [];  % just because I want get rid of nans here
y(ind) = [];  % just because I want get rid of nans here
f(ind) = [];  % just because I want get rid of nans here
samples = numel(t);

figure(1),
plot(t,f);

% generate cos and sin signals for DFT computation
% % method #1
% ang = cumsum(2*pi*f)*dt;
% xr = cos(ang);
% xi = sin(ang);

% method #2
xr = x./max(abs(x)); % cos signal of amplitude 1 
xi = gradient(xr); % sin signal of non constant amplitude must be divided by instantaneaous pulsation / frequency
xi = smoothdata(xi,'gaussian',5); % a bit of smoothing to reduce noise created by signal differentiation
xi = xi./f; % amplitude correction :amplitude must be divided by instantaneaous pulsation / frequency
xi = xi./max(abs(xi(round(samples/4):end))); % sin signal of amplitude 1 

% figure(2),
% plot(t,xr,'r',t,xi,'b');

% DFT coefficients for signal : x
Mr = 2*x.*xr; % factor 2 alike in FFT : X = 2/n*sum(x*e(-jwt))
Mi = 2*x.*xi;
% low pass filter real / imag amplitudes (DFT coefficients)
[B,A] = butter(2,2/Fs*0.1);
Mr = filter(B,A,Mr);
Mi = filter(B,A,Mi);
x_amplitude = sqrt(Mr.^2 + Mi.^2);
x_phase = atan2(Mi,Mr); % in rads

figure(3),
subplot(2,1,1),plot(t,x_amplitude);
title('X signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
subplot(2,1,2),plot(t,x_phase);
xlabel('Time (s)');
ylabel('Phase (rad)');

% DFT coefficients for signal : y
Mr = 2*y.*xr; % factor 2 alike in FFT : X = 2/n*sum(x*e(-jwt))
Mi = 2*y.*xi;
% low pass filter real / imag amplitudes (DFT coefficients)
Mr = filter(B,A,Mr);
Mi = filter(B,A,Mi);
y_amplitude = sqrt(Mr.^2 + Mi.^2);
y_phase = atan2(Mi,Mr); % in rads

figure(4),
subplot(2,1,1),plot(t,y_amplitude);
title('Y signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
subplot(2,1,2),plot(t,y_phase);
xlabel('Time (s)');
ylabel('Phase (rad)');

% now compute amplitude ratio and phase difference (like a Bode plot)
freq = 0:0.05:6;
[g,p] = dbode(B0,A0,1/Fs,2*pi*freq);
figure(5),
title('Transfer function Y/X')
subplot(2,1,1),plot(f,y_amplitude./x_amplitude,'b',freq,g,'r');
xlabel('Frequency (Hz)');
ylabel('Modulus ');
legend('order extraction method','Bode plot');
subplot(2,1,2),plot(f,unwrap(y_phase - x_phase),'b',freq,unwrap(p*pi/180),'r');
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
legend('order extraction method','Bode plot');


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

[Mathieu NOE]

problem solved ?