-1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ you can format code as code by selecting it and using the code format {} button in the editor; please do that to ease answering your question. Also, if you say "for some it works, for others it doesn't", trying to describe for which that is, is very very very likely key to your question... $\endgroup$ – Marcus Müller Apr 28 '18 at 11:59
  • 1
    $\begingroup$ I'll speculate your problem is likely due to where the FFT bin centers are relative to your $f_0$. But since you have provided no detailed cases that you consider good vs. bad, who can say? Anyway, why are you doing this in the frequency domain? Some simple trigonometry in the time domain is a more straightforward way to extract phase angle. $\endgroup$ – Andy Walls Apr 28 '18 at 12:16
  • $\begingroup$ @Marcus Muller Thanks for help. Like if i only change samp_f from 400e3 to 300e3 or 600e3, it will give correct phase difference. However, with samp_f of 300e3/600e3 if i induce phase shift of -170 deg in signal_2, again it shows wrong difference $\endgroup$ – naumankalia Apr 28 '18 at 12:48
  • $\begingroup$ Please edit your question to add these details in a comprehensive manner; I don't understand what you've wrote in your comment, it seems self-contradictory to me. Also, I was expecting you to also remove the superfluous empty lines while you format your code, to be completely honest and utterly frank. $\endgroup$ – Marcus Müller Apr 28 '18 at 12:51
  • $\begingroup$ @Andy Walls. Thanks for help. Will these trigonometry techniques work for noisy signals with multiple frequencies also? $\endgroup$ – naumankalia Apr 28 '18 at 12:54
0
$\begingroup$

If you measure the phase of a sinusoid that is of a frequency that is not exactly integer periodic in the FFT length, the FFT will be measuring the phase at a discontinuity in the sinusoid between the beginning and end of the FFT input aperture. This discontinuity can be either positive or negative in magnitude, causing the phase measurements to jump depending on where the FFT windows are located in the data stream.

If you do an FFTshift before the FFT (rotate the data by N/2), that will move the FFTs phase reference point to the center of your data window, where there is no discontinuity in a pure sinusoid, and thus the phase measurement will be less jumpy.

$\endgroup$
  • $\begingroup$ Thanks for help. I have tried the commands signal_1_fft = fftshift(Y1(1:end));signal_2_fft = fftshift(Y2(1:end)); after FFT but results are unpredictable i.e. for some phase difference it gives right answer and for some it does not. $\endgroup$ – naumankalia Apr 29 '18 at 1:26
0
$\begingroup$

You can get an exact measure of the phase (and amplitude) of a pure tone with non-integer frequency from DFT values using the method I describe in my blog article "Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1". In the presence of noise or other tones, you get a very good estimate. The effects of the other tones can be removed from the DFT patch once they are estimated by using the bin value formulas I give in another article "DFT Bin Value Formulas for Pure Real Tones".

If the frequency is known ahead of time, it is much simpler to frame the DFT on a whole number of cycles and calculate only the bin value of interest, then read the phase from the arg of the value.

Hope this helps.

Ced


Followup:

To figure out a frame size that holds a whole number of cycles, you merely have to figure out how many samples per cycle and use a whole multiple of that.

Call it $spc$ for samples per cycle.

The units of $f_s$ are samples per second.

The units of $f_0$ are cycles per second.

$$ spc = \frac{f_s}{f_0} $$

The units are:

$$ \frac{samples}{cycle} = \frac{\frac{samples}{second}}{\frac{cycles}{second}} $$

For $f_s = 300 KHz $ and $ f_s = 100 KHz $ means you have three samples per cycle. So, any multiple of 3 will work. If you have a steady tone, you will want several cycles to get a good read. The Nyquist frequency, the limit of the DFT, is at two samples per cycle. Obviously a multiple of 3 is never going to equal a power of 2, but don't worry about that. Modern FFT routines can do very well with non-power of two sizes as well. If you code the math yourself, you only need to evaluate one bin value, that bin number is the number of cycles per frame. Be aware that MATLAB has one based indexing and DFTs are zero based.

Since the DFT bin values are "circular" in regards to phase, if your time difference is greater than the number of samples per cycle you won't be able to tell how many cycles there are in between, only the fractional amount. The phase difference you measure using the DFT will be in terms of a fraction of the wavelength of one cycle. Thus you may want to test with a lower frequency tone which has more samples per cycle to calculate the time delay.

Sound travels at about 1130 feet per second. At 100 kHz, 100000 cycles per second:

$$ wavelength = \frac{1130 feet/second}{100,000 cycles/second} \cdot 12 \frac{inches}{foot} \approx 0.135 inches $$

If the distance between your mics is greater than that you will have more than a cycle difference.


Python Sample:

import numpy as np

N = 16384

f0=100e3
samp_f=400e3

cpf = N * f0 / samp_f

# cycles per frame = samples per frame * cycle per second 
#                  / samples per second


w = 2*np.pi * f0 / samp_f


# radians per sample = radians per cycle * cycle per second 
#                    / samples per second

# w = cpf * 2*np.pi / float( N ) # Alternative calculation

# radians per sample = cycles per frame  * radians per cycle
#                    / samples per frame

phase = .1234

signal_1 = np.zeros( N )

for n in range( N ):
  signal_1[n] = np.cos( w * n + phase )

dft = np.fft.rfft( signal_1 ) / float( N )

for d in range( -3, 4 ):
  bin = cpf + d
  print bin, dft[bin]

z = dft[cpf]
angle = np.arctan2( np.imag( z ), np.real( z ) )

print
print phase, angle

Here is the output. Notice the only bin with meaningful values is the one associated with the frequency.

4093.0 (-2.46600116419e-15-1.16625818025e-15j)
4094.0 (-1.1914697594e-14-1.77386979539e-14j)
4095.0 (-5.93014561131e-14-3.73808410544e-14j)
4096.0 (0.496197938352+0.0615435291056j)
4097.0 (6.67850950798e-14-2.19014433524e-14j)
4098.0 (1.58906441144e-14-1.41584611277e-14j)
4099.0 (2.90666468683e-15-5.09362957479e-16j)

0.1234 0.1234

The phase is recovered exactly.

$\endgroup$
  • $\begingroup$ Thanks for help. Is FFT not the fast implementation of DFT which i am doing in my code posted above ? As i know frequency i.e. fo=100 KHz and sampling frequency samp_f=300 KHz, how to chose chunk_size corresponding to whole number of cycles? $\endgroup$ – naumankalia Apr 29 '18 at 1:34
  • $\begingroup$ @naumankalia, Yes a FFT is a efficient implementation of a DFT. If you want to understand how a DFT works, search on that. Searching on "FFT" will get you results on the calculation efficiency, not the meaning of the calculations. See my follow up for more on your question. $\endgroup$ – Cedron Dawg Apr 29 '18 at 2:40
  • $\begingroup$ Thanks a lot for your help. As you explained, for fo=100 KHz and samp_f=400 Ksamples/s, chunk_size=16384 is exact multiple of 4, but still phase_difference measured has errors in it? Thanks $\endgroup$ – naumankalia Apr 29 '18 at 2:58
  • $\begingroup$ @naumankalia, You should be able to follow my Python sample I just added. $\endgroup$ – Cedron Dawg Apr 29 '18 at 3:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.