I am trying to understand how to find phase difference b/w two real sinusoidal signals and then align the phases of two sinusoidal signals using this difference. In the first step, however, i have to find the phase difference first. By studying this problems on different forums, i have written a MATLAB script to simulate the problem:

f0=100e3;
samp_f=300e3;
samp_t=1/samp_f;
chunk_size=16384;

phase_shift=-170*pi/180;

NFFT=chunk_size;
chunk_t=0:samp_t:(samp_t*(chunk_size-1));

%Frequency vector calculation
 frq_vec=zeros(1,NFFT/2);
 for v=0:NFFT/2
     f(v+1)=v/(samp_t*NFFT);
 end

signal_1=sin(2*pi*f0*chunk_t);
Y1 = fft(signal_1,NFFT)/chunk_size;

signal_2=sin(2*pi*f0*chunk_t+phase_shift);
Y2 = fft(signal_2,NFFT)/chunk_size;

signal_1_fft = Y1(2:NFFT/2);
signal_2_fft = Y2(2:NFFT/2);

signal_1_phase = unwrap(angle(signal_1_fft));
signal_2_phase = unwrap(angle(signal_2_fft));

signal_phase_diff = signal_2_phase - signal_1_phase;

% Convert from bin number to frequency from 0 to Fs
desired_bin = ((length(signal_1_fft)*2)/samp_f)*f0;
desired_bin = round(desired_bin)+1 %as DC was removed

phase_diff = signal_phase_diff(desired_bin)*180/pi

%Results of code for different combination of parameters
%For samp_f=400e3;phase_shift=-100*pi/180;  >>phase_diff=1.0235e+03 Wrong
%For samp_f=300e3;phase_shift=-100*pi/180;  >>phase_diff=-100.0029  Right
%For samp_f=300e3;phase_shift=-170*pi/180;  >>phase_diff=189.9971   Wrong

The problem i am facing is that for some combinations of samp_f, chunk_size and phase shifts incorporated in signal_2, i got correct phase_diff and for some i don't. Can any one kindly guide me what am i doing wrong in this code?

Thanks

Just a quick guess:

signal_2=sin(2*pi*f0*chunk_t-100*pi/180);

shoud be

signal_2=sin(2*pi*f0*(chunk_t-100*pi/180));

some other issues:

you are entering phase in time domain as samples but fft outputs phase in radians.

The tone is going to leak to adjacent bins unless you choose full cycles

you can align in time domain without going to frequency domain

Exactly. Also keep in mind that the range of the phase you got using FFT usually wrapped between [-0.5pi +0.5pi], which means that you will not possibly measure the correct phase difference if it is larger than a half of the period.

Thanks for help . If i understand correctly, i cannot measure accurate phase difference b/w two sinusoidal signals using FFT method if the phase difference b/w them varies more than +/- 90 deg? Please correct me if i am wrong because in that case i have to go to some other technique as phase variations can definitely vary more than +/- 90 deg.

Thanks

I have this quick code to study phase computation from fft. As far as I notice I can detect phase within 1 cycle ratio(with some wrap up issue) but understandably it can't detect how many integer cycles it has changed.

I use ifft to generate exact cycles conveniently for this study.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; clc;

bin = 11;        %generate frequency at this bin

s = .3;          %shift phase by this ratio of full cycle

x = zeros(1,1000);;

x(bin) = 1;

%generate waveform

x1 = real(ifft(x));

%shift phase of above waveform

samples_per_cycle = 1000/(bin-1);

shift = floor( s * samples_per_cycle)+0*samples_per_cycle;

x2 = [x1(end-shift+1:end) x1(1:end-shift)];

figure;

plot(x1(1:200),'.-');

hold;

plot(x2(1:200),'r.--');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%calculate phase shift using fft 

y1 = fft(x1);

y2 = fft(x2);

p1 = unwrap(angle(y1)); %zero if cosine 

p2 = unwrap(angle(y2));

ph_diff = (p1(bin) - p2(bin))/(2*pi);  %fft result in samples

if ph_diff < 0, ph_diff = 1+ph_diff;end;

ph_diff

shift/samples_per_cycle               %actual entered

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thanks for help. I am trying to implement sin(wt+phi) formula. I know my frequency i.e. 100 KHz and my sampling rate e.g. 400 Ksamples/s, kindly guide me how to chose "chunk_size" corresponding to full cycles.

Thanks again

you my use windowing before fft so that cycle alignment is not required. choose the window type, for phase, hanning should be a good try. 

in addition to this you may try equivalent of parabolic interpolation to accurate phase (parabolic is used for amplitude). you need to find an equivalent one for phase. may lead to a similar approaches like cFFT described in cFFT and other by Biancacci. 

The actual problem with the code given is use of the unwrap() function.  This is not needed nor appropriate for the results of an FFT across bins. unwrap() is for use with smoothly varying phase functions that go through 2pi. The simplified code below produces the correct result.  However there are some concerns about selection of a single bin.  What happens if two bins share the energy equally because of bin size and signal frequency (spectral leakage)?  I suspect the phase difference in both bins will be essentially the same.

f0=100e3;

samp_f = 400e3;

samp_t = 1/samp_f;

chunk_size = 16384;

phase_shift = -100*pi/180;

NFFT = chunk_size;

chunk_t = 0:samp_t:(samp_t*(chunk_size-1));

signal_1 = sin(2*pi*f0*chunk_t);

Y1 = fft(signal_1,NFFT)/chunk_size;

signal_2=sin(2*pi*f0*chunk_t+phase_shift);

Y2 = fft(signal_2,NFFT)/chunk_size;

% Convert frequency to bin number

desired_bin = round(f0 / samp_f * chunk_size) + 1;

% get phase from signal bin

signal_1_phase = angle(Y1(desired_bin));

signal_2_phase = angle(Y2(desired_bin));

% calculate phase difference

phase_diff = (signal_2_phase - signal_1_phase) * 180 / pi;

% restrict phase to +- 180 degrees

if (phase_diff > 180)

phase_diff = phase_diff - 360;

elseif (phase_diff < -180)

phase_diff = phase_diff + 360;

end

% report difference

fprintf('phase difference = %8.4f\n', phase_diff);

%Results of code for different combination of parameters

%For samp_f=400e3;phase_shift=-100*pi/180; >>phase_diff= -100.0000

%For samp_f=300e3;phase_shift=-100*pi/180; >>phase_diff= -99.9986

%For samp_f=300e3;phase_shift=-170*pi/180; >>phase_diff= -169.9986