How to determine phase shift between two waveforms?
206 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi,
I want to determine the phase shift between two waveforms. One of the waveform is considered as pivot/fixed and the second one is varying. Since the variation is both horizontally as well as vertical, its very difficult to determine the phase shift. Thus one of the community member told me to use the concept of zero crossing.
But it gives me an error, which shows array mismatch, which I'm unable to fix and need your suggestion.
The code is as attached herewith.
Sincerely,
rc
2 commentaires
Réponse acceptée
Star Strider
le 4 Juil 2023
I am not certain what you want.
The easiest way too analyse the phase differences is to do a Fourier transform of the signals, and plot them together. The phases are quite close together (differing by a small amount) up to about 22 Hz, and diverge significantly beyond that. There is a spike at about 39.04 Hz,, and the phase separation there is relatively constant, with a relatively constant phase difference —
% type('OBP1N.txt')
T1 = readtable('OBP1N.txt', 'HeaderLines',8, 'VariableNamingRule','preserve');
t1 = T1{:,1};
s1 = T1{:,2};
T2 = readtable('OBP2N.txt', 'HeaderLines',8, 'VariableNamingRule','preserve');
t2 = T2{:,1};
s2 = T2{:,2};
Ts1 = mean(diff(t1));
Tssd1 = std(diff(t1));
Fs1 = fix(1/Ts1);
[sr1,tr1] = resample(s1,t1,Fs1);
L1 = numel(tr1);
Ts2 = mean(diff(t2));
Tssd2 = std(diff(t2));
Fs2 = fix(1/Ts2);
[sr2,tr2] = resample(s2,t2,fix(1/Ts2));
L2 = numel(tr2);
Fn = Fs1/2;
NFFT = 2^nextpow2(L1);
FTs12r = fft(([sr1 sr2]-mean([sr1 sr2])).*hann(L1), NFFT)/L1;
Fv = linspace(0, 1, NFFT/2+1).'*Fn;
Iv = 1:numel(Fv);
figure
subplot(2,1,1)
plot(Fv, abs(FTs12r(Iv,:))*2)
grid
ylabel('Magnitude)')
subplot(2,1,2)
plot(Fv, rad2deg(unwrap(angle(FTs12r(Iv,:)))))
grid
xlabel('Frequency')
ylabel('Phase (°)')
sgtitle('Fourier Transform')
figure
subplot(2,1,1)
plot(Fv, abs(FTs12r(Iv,:))*2)
grid
xlim([0 25])
ylabel('Magnitude)')
subplot(2,1,2)
plot(Fv, rad2deg(unwrap(angle(FTs12r(Iv,:)))))
grid
xlim([0 25])
xlabel('Frequency')
ylabel('Phase (°)')
sgtitle('Fourier Transform')
[pks, locs] = findpeaks(abs(FTs12r(Iv,1))*2, 'MinPeakProminence',0.05);
Fpeak = Fv(locs);
Phs = rad2deg(unwrap(angle(FTs12r(locs,:))));
fprintf('Peak Frequency = %8.3f Hz\nPhase 1 = %8.3f°\nPhase 2 = %8.3f°\nPhase Difference = %8.3f°\n',Fpeak, Phs, diff(Phs))
figure
subplot(2,1,1)
plot(Fv, abs(FTs12r(Iv,:))*2)
grid
xlim([38.90 39.20])
ylabel('Magnitude)')
subplot(2,1,2)
plot(Fv, rad2deg(unwrap(angle(FTs12r(Iv,:)))))
grid
xlim([38.90 39.20])
xlabel('Frequency')
ylabel('Phase (°)')
sgtitle('Fourier Transform')
% return
sr1 = filloutliers(sr1,'linear','grubbs');
sr2 = filloutliers(sr2,'linear','grubbs');
% [eu,el] = envelope(sr1, 100, 'peak');
%
% figure
% plot(tr1, sr1)
% hold on
% plot(tr1, eu)
% hold off
% grid
figure
plot(tr1, sr1)
hold on
plot(tr2, sr2)
hold off
grid
That is the best I can do with this analysis.
.
10 commentaires
Star Strider
le 18 Juil 2023
As always, my pleasure!
Use the find function to get the index. (The findpeaks function returns it as ‘locs’ the way I call it.)
The 80 Hz example would be:
idx = find(Fv <= 80, 1, 'last');
or:
idx = find(Fv >= 80, 1);
then:
Fpeak = Fv(idx);
Phs = rad2deg(unwrap(angle(FTs12r(idx,:))));
fprintf('Peak Frequency = %8.3f Hz\nPhase 1 = %8.3f°\nPhase 2 = %8.3f°\nPhase Difference = %8.3f°\n',Fpeak, Phs, diff(Phs))
to get the same result as perviously, at the chosen frequency.
.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Spectral Measurements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!