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

조회 수: 7 (최근 30일)
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.
  댓글 수: 2
Mathieu NOE
Mathieu NOE 2023년 9월 28일
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
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;
ylabel('Mag (dB)');
ylabel('Phase (°)');
xlabel('Frequency (Hz)');
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 = 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;
% 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];
select = [1:nfft/2+1]; % include DC AND Nyquist
Pxx = Pxx(select);
Pyy = Pyy(select);
Pxy = Pxy(select);
select = 1:nfft;
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;

댓글을 달려면 로그인하십시오.

답변 (1개)

Mathieu NOE
Mathieu NOE 2023년 9월 28일
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);
% 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
title('X signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
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
title('Y signal modulus / phase extraction')
xlabel('Time (s)');
ylabel('Modulus ');
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));
  댓글 수: 2

댓글을 달려면 로그인하십시오.


Help CenterFile Exchange에서 Transforms에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by