How to determine phase shift between two waveforms?

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

[Sushma Swaraj]

Can you provide your data file?

[Rahul]

Hi,

Yes.

Attached herewith the data files.

rgds,

rc

Answer 1

[Star Strider]

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.

.

[Rahul]

Hi Star Strider,

Thank you very much. Your reply always helps.

With regards,

rc

[Star Strider]

As always, my pleasure!

I very much appreciate your compliment!

[Rahul]

Hi Star Strider,

1) Can you please elaborate this step as below:

FTs12r = fft(([sr1 sr2]-mean([sr1 sr2])).*hann(L1), NFFT)/L1;

what is happening here, I didn't get it at all.

2) Fv = linspace(0, 1, NFFT/2+1).'*Fn;

Is it transpose you took here (.')?

3) Why did you multiply by factor 2 in

plot(Fv, abs(FTs12r(Iv,:))*2)

with rgds,

rc

[Star Strider]

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.

.

[Rahul]

Hi Star Strider,

Thanks a lot for your detailed explanation to my query.

I have one more query.

For the files 'OBP1N' and 'OBP2N' with your valuable help the spectrograms obtained are as below:

Now, I want to determine the phase shift only at frequency 40 kHz, (since the spectrogram shows something at this frequency value happening) and hence not for the whole system.

Do I need to put a filter again or is there any other tricks available.

Please advice.

With regards,

rc

[Star Strider]

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]

Hi Star Strider,

This is because of the following,

data = readmatrix(file, 'HeaderLines',8);

data(:,1) = data(:,1)*1E-3 ;

I multiplied 1E-3 as above. When I remove the frequency comes in Hz range.

Kindly show me how the phase shift can be determined only for 40Hz frequency.

With regards,

rahul

[Star Strider]

‘Kindly show me how the phase shift can be determined only for 40Hz frequency.’

I already did that originally:

[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°

The peak is actually at 39.04 Hz, however the phases of the two signals do not vary significantly in that region, and so should be the same for 40 Hz. The ‘Phase Difference’ can be ±32.311° depending on how you calculate it.

.

[Rahul]

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]

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.

.