Matlab - calculate phase shift using fft

  • Thread starter daymos
  • Start date
  • #1
8
0
hello,

I have 3 signals in the form of sampled values. When I plot them using

plot (t,vPa,t,vPb,t,vPc)

where vPa, vPb, vPc contains the values and t contains the sampling istants I get this:



tttt.jpg



when I calculate phase shift using fft I get phase angle = 0.
I have tried using my own code ( which works normally on other vectors) and also a function downloaded from Mathworks and the result is always zero.

The mathworks function is at the end.

what am I doing wrong?

Thanks
Mattia


Matlab:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Phase Difference Measurement            %
%               with MATLAB Implementation             %
%                                                      %
% Author: M.Sc. Eng. Hristo Zhivomirov        12/01/14 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function PhDiff = phdiffmeasure(x, y)

% function: PhDiff = phdiffmeasure(x, y)
% x - first signal in the time domain
% y - second signal in the time domain
% PhDiff - phase difference Y -> X, degrees

% represent x as column-vector if it is not
if size(x, 2) > 1
    x = x';
end

% represent y as column-vector if it is not
if size(y, 2) > 1
    y = y';
end

% signals length
xlen = length(x);
ylen = length(y);

% window preparation
xwin = hanning(xlen, 'periodic');
ywin = hanning(ylen, 'periodic');

% fft of the first signal
X = fft(x.*xwin);

% fft of the second signal
Y = fft(y.*ywin);

% phase difference calculation
[~, indx] = max(abs(X));
[~, indy] = max(abs(Y));
PhDiff = angle(Y(indy)) - angle(X(indx));
PhDiff = PhDiff*180/pi;

end
 

Answers and Replies

  • #2
12,811
6,693
I don't see what you did wrong.

However, have you looked at the intermediate values for X and Y and the angle(X(indx)) values that are used to compute the phase difference?
 
  • #3
8
0
I have now, angle(Y(indy)) and angle(X(indx)) are both zero.

I have attached the original 3 vectors.
 

Attachments

  • vectors.txt
    26.1 KB · Views: 589
Last edited:
  • #4
8
0
I have found a different way to calculate phase shift that can be used in this case:

h1 = hilbert(vPa);

h2 = hilbert(vPb);

phase_rad = max(angle(h1 ./ h2))

which gives back the correct value of 2/3 * pi.

However I would be couriois to know why for this specific case fft didnt work. I thought it was a pretty generic method
 
  • #5
12,811
6,693
I meant the X and Y in the function you showed.

Aren't they supposed to be mostly zeros with a few spikes?

I think you need to look at each line of the function to decide if its doing the right thing.

The code looks like it applies an FFT to the two input signals to produce a frequency spectrum. From that it determines the highest spike and from that subtracts them to determine the phase difference.
 
  • #6
8
0
yes that is what it does. Also my own code did that. ffts are mostly zero with one big spike on 50 Hz.
I am positive both my own code and the one I posted are correct. becaouse I have used them many times and they returned correct results
 
  • #7
12,811
6,693
Then I guess there's something unique about your input data and the function you listed. Perhaps if you step through the function interactively you'll be able to spot the difference. It could be a number of elements limit or datatype mismatch...
 
  • #8
12,811
6,693
Perhaps your data doesn't have an FFT with a significant frequency that can be picked out i.e. no one max value

Try plotting the FFT's of both signals to see what you get.
 
  • #9
8
0
thanks for the spectrum plot tips jedishrfu
I haven't solved it yes but I think you pointed me in the right direction.

this is the spectrum of the 3 signals:

spectrum.jpg



So there is a huge spike at zero Hz.

I have composed each signal multiplying 2 previous signals that where at 50Hz. So I expect the new signals to be at 100 Hz. And I see a spike there.
But I dont understand why the spike in zero.
 
Last edited:
  • #10
155
24
You have a DC component in you signal. max function should detected the peak of the periodic signal ?
 
  • #11
8
0
You have a DC component in you signal. max function should detected the peak of the periodic signal ?


when I use max function it gives back the index of frequency = 0 HZ
 
  • #12
155
24
Then the answer is correct, the DC components have zero phase shift :)
 
  • #13
8
0
spectrum1.jpg
Then the answer is correct, the DC components have zero phase shift :)

well, ok thats correct.

Since I need the 100Hz I subtracted the mean value of the signal to itself.



thanks a lot!
 
  • #14
105
4
The zero value of the fourier transform will always just be the mean give a or take a factor of pi and volume.
 
  • #15
boneh3ad
Science Advisor
Insights Author
Gold Member
3,262
926
In general, it's usually a good idea to make sure the signal you wish to pass through your FFT has zero mean, that way you are only dealing with nonzero frequency components.
 
  • #16
8
0
the original 2 signals had no DC component. But by multiplying them value by value the resulting signal was no longer centered on the zero for the vertical axis.
Pretty obvious after you get it..
 
  • #17
1
0
what did you change ad your coding to make that output? thanks in advance :)
 

Related Threads on Matlab - calculate phase shift using fft

  • Last Post
Replies
11
Views
53K
  • Last Post
Replies
3
Views
8K
Replies
1
Views
6K
Replies
2
Views
8K
Replies
17
Views
11K
Replies
16
Views
499
  • Last Post
Replies
5
Views
8K
  • Last Post
Replies
3
Views
3K
Replies
7
Views
12K
Replies
1
Views
5K
Top