Phase difference b/w two sinusoidal signals using FFT
Started by 7 years ago●8 replies●latest reply 7 years ago●2865 viewsI 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