How to determine phase shift between two waveforms?

Question

Rahul 2023-7-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1991413-how-to-determine-phase-shift-between-two-waveforms

评论： Star Strider 2023-7-18

spectrum_tt1_obp4n.m

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 个评论
显示无隐藏无

Sushma Swaraj 2023-7-4

Can you provide your data file?

Rahul 2023-7-4

Hi,

Yes.

Attached herewith the data files.

rgds,

rc

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-7-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1991413-how-to-determine-phase-shift-between-two-waveforms#answer_1266963

在 MATLAB Online 中打开

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

Peak Frequency = 39.040 Hz Phase 1 = -101.825° Phase 2 = -134.136° Phase Difference = -32.311°

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 个评论
显示 8更早的评论隐藏 8更早的评论

Star Strider 2023-7-13

The fft call concatenates the ‘sr1’ and ‘sr2’ vectors (for convenience and efficiency), subtracts their means, so that any d-c (constant) offset does not make any other peaks difficult to see, then windows that result using a Hanning window (windows correct for the fft not being infinite, improving the resolution, the Hanning window is the most commonly used), uses ‘NFFT’ as the fft length (for efficiency (powers of 2 are more appropriate because of the way the fft is calculated), and then normalises that result by dividing it by the length of the original vector, necessary again because of the way the fft is calculated.

The ‘Fv’ assignment calculates the frequency vector for the subsequent plot and the findpeaks calls. There are several ways to calculate it, this was introduced in the R2015b documentation (if I remember correctly) and since it is easy to remember, I almost always use it. (The ‘Iv’ vector is the index vector.) I used the transpose because the other vectors are also column vectors.

The fft divides the energy in the original signal between the ‘positive’ and ‘negative’ frequencies of a two-sided Fourier transforom (this is easiest to visualise using the fftshift function). In a one-sided fft, multiplying by 2 corrects for this, giving a better approximation of the fft magnitudes with respect to the original signal.

.

Star Strider 2023-7-14

在 MATLAB Online 中打开

I do not understand how you are getting frequencies in the 40 kHz range. When I re-ran my code (to be certain that I understood it again), the Nyquist frequency is only 100 Hz, and the peaks are about 40 Hz. My pspectrum plots (added to my previous code) appear to agree with those values.

In any event, find the frequency that most closely matches your frequency-of-interest to get its index, and then extract the phase values, and calculate the phase difference as I did previously.

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

Peak Frequency = 39.040 Hz Phase 1 = -101.825° Phase 2 = -134.136° Phase Difference = -32.311°

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

figure

pspectrum(s1, t1, 'spectrogram') % Original Data

colormap(turbo)

figure

pspectrum(sr1, tr1, 'spectrogram') % Resampled Data

colormap(turbo)

figure

pspectrum(s2, t2, 'spectrogram') % Original Data

colormap(turbo)

figure

pspectrum(sr2, tr2, 'spectrogram') % Resampled Data

colormap(turbo)

.

Rahul 2023-7-18

Hi Star Strider,

Thank you very much. I have one more query.

If I wish to determine the phase difference at some other value of frequency (say at 80 Hz) then how can this be done in this code.

Is it possible to have the value of frequency at which the phase difference is to be determined is given as input by the user and then the phase difference is determined for that input frequency?

with regards,

rc

Star Strider 2023-7-18

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.

.

请先登录，再进行评论。

How to determine phase shift between two waveforms?

2 个评论
显示无隐藏无

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to determine phase shift between two waveforms?

2 个评论 显示 无隐藏 无

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

10 个评论
显示 8更早的评论隐藏 8更早的评论