How to determine phase shift and slope for a varying waveform?

Hi,

I have 12 nos. of probes which detects and picks signals (namely obp1n to obp12n). After determining the spectogram and filtering the waveform, I want to keep obp1n as my reference probe and determine phase shift of signals from others with reference to obp1n and hence determine the slope.

Now, the problem is the amplitudes of these waveforms are varying from one probe to another, also they are shifting horizontally (along time axis). Thus I'm not getting the proper determination of phase shift.

I'm not getting the correct value for phase shift and hence slope.

Kindly advice.

rgds,

rc

f = @(filename) fullfile('D:\Data_tokamak\TT1 data',filename );
file = f('OBP6N.txt');
%filecontents = fileread('OBP1N.txt');                        % Check File Format & Contents

data = readmatrix(file, 'HeaderLines',8);                     % Skip First 8 Header Lines
data(:,1) = data(:,1) * 1E-3;                                 % Time in milliseconds  

nonfinite_rows = nnz(~any(isfinite(data),2));                 % Check To See If Any Rows Have Non-Finite Values
if nonfinite_rows > 0
    data(~isfinite(data)) = NaN;
    data = fillmissing(data, 'linear');
end

x_data = data(:,1);
y_data = data(:,2);

Ts = mean(diff(x_data));
Tsd = std(diff(x_data));

Fs = 1/Ts;

[y_data,x_data] = resample(y_data, x_data, round(Fs));          % Resample To Constant Sampling Frequency

Fs = round(Fs);
Fn = Fs/2;

figure(1)
pspectrum(y_data, Fs, 'spectrogram');                           % Plot 'pspectrum' 'spectrogram'
colormap(turbo);                                                 % Introduced: R2020b
[sp,fp,tp] = pspectrum(y_data,Fs,"spectrogram");                % Return Values, Then Plot


figure(2)
waterfall(fp,tp,sp');
set(gca,XDir="reverse",View=[60 60]);
colormap(turbo);                                                 % Introduced: R2020b
ylabel("Time (s)");
xlabel("Frequency (Hz)");

figure(3)
waterfall(fp,tp,sp');
set(gca,XDir="reverse",View=[60 60]);
Ax = gca;
Ax.ZScale = 'log';
colormap(turbo);                                                 % Introduced: R2020b
ylabel("Time (s)");
xlabel("Frequency (Hz)");

% Design filter for studying activities at particular frequency
wp  =  [39.5  40.5]/Fn;                                              % Passband frequency Normalized
ws  =  [0.98  1.02].*wp;                                        % Stopband frequency Normalized
Rp  =  1;                                                       % Passband Ripple
Rs  =  60;                                                      % Passband Ripple (Attenuation)
[n,wp] = ellipord(wp,ws,Rp,Rs);                                 % Elliptic order calculation
[z,p,k] = ellip(n,Rp,Rs,wp);                                    % Elliptic filter design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k);                                        % Second order selection for stability




figure(4)
freqz(sos, 2^18,Fs)    ;                                         % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 100])                                % OPTTIONAL
set(subplot(2,1,2), 'XLim',[0 100])                                % OPTTIONAL



[h w] = freqz(sos, 2^18,Fs)    ; 


%Implementation of filtfilt function
y_data_filt = filtfilt(sos, g, y_data);                         % Filter Signal
figure(5)
yyaxis right
plot(x_data, y_data_filt,'.')
grid
xlabel('Time')
ylabel('Amplitude')
hold on

Answer 1

[Mathieu NOE]

hello

I would suggest to compute the zero crossing time coordinates of both signals

the function used here implement linear interpolation for higher accuracy (better than one sample dt)

here a small demo of 2 signals of different amplitudes and shifted by 0.2 rad

you see the computation gives the same result as what we used in first place for the data generation

x = 0:0.1:30;
y1 = sin(x+0.25).*exp(-x/10); % reference signal 
y2 = 0.75*sin(x+0.45).*exp(-x/10);  % signal shifted by 0.2 rad

figure(1),
plot(x,y1,'r',x,y2,'b');

% find crossing points at y threshold = 0
threshold = 0;
[xc1] = find_zc(x,y1,threshold);
[xc2] = find_zc(x,y2,threshold);

% phase shift  = x distance divided by period and multiplied by 2pi
p1 = mean(diff(xc1)); % mean period of signal 1
p2 = mean(diff(xc2)); % mean period of signal 2

phase_shift = mean(xc1-xc2)/p1*2*pi % is the value used in the data generation 

phase_shift = 0.2000

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

[Rahul]

Hi Mathiew NOE,

Thanks a lot for your reply.

1) I'm getting an error as follows:

Arrays have incompatible sizes for this operation.

Error in spectrum_tt1_obp4n (line 163)

phase_shift = mean(xc1-xc2)/p1*2*pi;

2) I couldn't follow the line in the code where the function is defined as

zci = @(data) find(diff(sign(data))>0)

How does "data" receives its input?

Need your advice plz.

Sincerely,

rc

[Mathieu NOE]

hello

1/ ok I admit here is maybe a litlle improvement to do on the code as xc1 and xc2 may have different length depending of your data (the number of zero crossing points may differ between the reference and second signal)

do you have some data files where the error appears to share ? I will fix my code for that cases

2/ this is anonymous function call in matlab - see doc / examples

Anonymous Functions - MATLAB & Simulink - MathWorks France

Pass Function to Another Function - MATLAB & Simulink - MathWorks France

Anonymous Functions in Matlab | Working | Steps | Advantage (educba.com)

[Mathieu NOE]

hello again

try this updated code . i intentionnaly created a situation where the number of zero crossing points would be different between the two signals. It works now in this case

nevertheless it is always good practice to make sure the signals are clean (smooth) enough so there is no "multiple" zero crossing detections (due to zig zags in the signal (noise)) where it should be only one

x = 0:0.1:30;
y1 = sin(x+0.25).*exp(-x/10); % reference signal 
y2 = 0.75*sin(x+0.45).*exp(-x/10);  % signal shifted by 0.2 rad

% remove some data from y2 to make it shorter than y1 
y2(1:100) = NaN;
y2(end-100:end) = NaN;

figure(1),
plot(x,y1,'r',x,y2,'b');

% find crossing points at y threshold = 0
threshold = 0;
[xc1] = find_zc(x,y1,threshold);
[xc2] = find_zc(x,y2,threshold);
p1 = mean(diff(xc1)); % mean period of signal 1
p2 = mean(diff(xc2)); % mean period of signal 2
p = mean([p1 p2]); % mean period of both signals (should be anyway very close)

% as xc1 and xc2 may have different length (depending of signals) we have
% to pick the common (close by) values and use only them

dist = abs(xc1-xc2'); % distance matrix between xc1 and xc2
id = (dist<p/2);  % we want zc points that are distant not above period/2
ida = logical(sum(id,1));
idb = logical(sum(id,2)');

% the valid xc1 and xc2 are then : 
xc1 = xc1(ida);
xc2 = xc2(idb);

% phase shift  = x distance divided by period and multiplied by 2pi
phase_shift = mean(xc1-xc2)/p1*2*pi % is the value used in the data generation 

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

[Rahul]

Hi Mathiew NOE,

1. Please find attached the input txt files with data. Both files contains exactly same no. of samples.

But from workspace it seems p1 and p2 are not equal. But how??

2. I didn't get how distance matrix is evaluated. I'm again getting an error as follows:

Arrays have incompatible sizes for this operation.

Error in spectrum_tt1_obp4n (line 167)

dist = abs(xc1-xc2');

rgds,

rc

[Mathieu NOE]

hello

this is the plot of your data

seems quite different from the clean damped sinusoids I expected (similar to what you posted in first place)

what are we supposed to do with such data ??

[Rahul]

Hi Mathiew NOE,

In the passband frequency it gives the plot as under. But after that it shows the error as I mentioned before. Here, I need your suggestion plz in solving the error.

wp2  =  [39.5  40.5]/Fn2;                                              % Passband frequency Normalized
ws2  =  [0.98  1.02].*wp2;                                             % Stopband frequency Normalized
Rp2  =  1;                                                             % Passband Ripple
Rs2  =  60;                                                            % Passband Ripple (Attenuation)
[n2,wp2] = ellipord(wp2,ws2,Rp2,Rs2);                                  % Elliptic order calculation
[z2,p2,k2] = ellip(n2,Rp2,Rs2,wp2);                                    % Elliptic filter design: Zero-Pole-Gain
[sos2,g2] = zp2sos(z2,p2,k2);                                          % Second order selection for stability

[Mathieu NOE]

hello

yes I can imagine you apply this filter but what portion of the original signal are you showing here ?

this is 0.5 s duration, whereas the full record is 500 s long

so where did you extract this 0.5 s long array ?