How to determine phase shift and slope for a varying waveform?
10 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Mathieu NOE
on 28 Jun 2023
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
7 Comments
Mathieu NOE
on 5 Jul 2023
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 ?
More Answers (0)
See Also
Categories
Find more on Single-Rate Filters in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!