False phase lock detector for PLL (GPS application domain)

In continuation of the topic I started earlier, which is about the false ph=
ase lock in PLL (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew9=
8E).

I am trying to implement a false phase lock detector using the approach des=
cribed in Kaplan book "Understanding GPS: Principles and Applications": htt=
p://tinyurl.com/z27gl4v=20

The idea is to compute the phase angle between the in-phase and quadrature =
pair of samples (I/Q) from the prompt replica correlator within 20 ms data =
bit interval for L1 C/A GPS signal:

C =3D cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
D =3D dot product of the same vectors=20
angle =3D atan2(C, D)

The angle divided by the time 'T'. over which it was accumulated, should gi=
ve a frequency:

F =3D angle / (2 * pi * T) [Hz]

Note that the time interval 'T' is the time between (Ip-1,Qp-1) and (Ip,Qp)=
 readings and it is slightly less than 20ms. In my case, the (Ip-1,Qp-1) re=
ading is done at ~300us and the (Ip,Qp) reading is done at 20ms after the d=
ata bit edge. Therefore, T is ~19.7 ms.
The PLL discriminator sees only the sum of two vectors
(Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.

The computed correction F can be applied to PLL state to correct the false =
lock error.

This approach seems to be working in the sense that, if I apply it to the c=
orrelator readings with the PLL in the false lock condition, then it gives =
me ~12.5 Hz error estimation with a correct sign. At the same time, I know =
that the true PLL error is 25Hz. BTW, if I do the error computation with th=
e true error being 75Hz instead of 25Hz, the result is still ~12.5Hz.

The thing, which I cannot figure out is why the approach gives me only the =
half of the true error: 12.5 instead of 25 Hz. Any ideas are very much appr=
eciated.

BTW, Kaplan talks about the threshold of 15.5Hz and then the 25 correction =
with the proper sign is applied. However my results never reach that thresh=
old.

On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:

> This approach seems to be working in the sense that, if I apply it to the correlator readings with the PLL in the false lock condition, then it gives me ~12.5 Hz error estimation with a correct sign. At the same time, I know that the true PLL error is 25Hz. BTW, if I do the error computation with the true error being 75Hz instead of 25Hz, the result is still ~12.5Hz.

By the way, which value do you get as the output of false lock detector 
in case the PLL is in true lock?

Gene

When PLL is in true lock, the false lock detector returns zero. So it works as expected in this case.

On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:
> In continuation of the topic I started earlier, which is about the false phase lock in PLL (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew98E).
>
> I am trying to implement a false phase lock detector using the approach described in Kaplan book "Understanding GPS: Principles and Applications": http://tinyurl.com/z27gl4v
>
> The idea is to compute the phase angle between the in-phase and quadrature pair of samples (I/Q) from the prompt replica correlator within 20 ms data bit interval for L1 C/A GPS signal:
>
> C = cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
> D = dot product of the same vectors
> angle = atan2(C, D)
>
> The angle divided by the time 'T'. over which it was accumulated, should give a frequency:
>
> F = angle / (2 * pi * T) [Hz]
>
> Note that the time interval 'T' is the time between (Ip-1,Qp-1) and (Ip,Qp) readings and it is slightly less than 20ms. In my case, the (Ip-1,Qp-1) reading is done at ~300us and the (Ip,Qp) reading is done at 20ms after the data bit edge. Therefore, T is ~19.7 ms.
> The PLL discriminator sees only the sum of two vectors
> (Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.
>
> The computed correction F can be applied to PLL state to correct the false lock error.
>
> This approach seems to be working in the sense that, if I apply it to the correlator readings with the PLL in the false lock condition, then it gives me ~12.5 Hz error estimation with a correct sign. At the same time, I know that the true PLL error is 25Hz. BTW, if I do the error computation with the true error being 75Hz instead of 25Hz, the result is still ~12.5Hz.
>
> The thing, which I cannot figure out is why the approach gives me only the half of the true error: 12.5 instead of 25 Hz. Any ideas are very much appreciated.
>
> BTW, Kaplan talks about the threshold of 15.5Hz and then the 25 correction with the proper sign is applied. However my results never reach that threshold.
>

If I have understood you correctly, you are doing the following thing. 
Assume that one 20ms time interval starts from t_1 = 0 and ends at t_4 = 
20ms. The vector (Ip-1, Qp-1) is obtained as the correlator output for 
the time interval from t_1 = 0 to t_2 = 300us, and the vector (Ip, Qp) 
is the correlator output for the time interval from t_3 = 300 us to t_4 
= 20 ms. Right?

Now if I got that part right, there's nothing surprising in your result! 
Because the parameter T in your formula for frequency estimation is the 
time interval between CENTERS of the respectful time intervals (t_1; 
t_2) and (t_3; t_4). Which is

T = 0.5*(t_4 - t_3) - 0.5*(t_2 - t_1) = 9.7 ms

rather than the value you assumed (19.7 ms). If you use the correct 
value of T = 9.7 ms, you'd get the correct frequency reading.

So, why the centers of the time intervals, and not something else, like 
ending positions as you have assumed? It helps to understand how do you 
obtain that value. You have two vectors (Ip-1, Qp-1) and (Ip, Qp); and 
in atan2() operation you obtain the angle between these two vectors. It 
means that for purposes of frequency estimation you only care about the 
direction of each vector. What is the direction of a vector? In your 
case, the direction of each vector is the PHASE ESTIMATE obtained within 
the corresponding time interval (note that in your case you obtain phase 
estimate in the form of a phasor, i.e. a complex exponent exp^(j*angle) ).

Now the question is, if you measure phase estimate for the correlation 
interval t=(t1; t2), assuming constant frequency offset f, FOR WHICH 
MOMENT OF TIME would you get a phase estimate? Let's start with the 
system model. Assume phase evolves with time as

phase(t) = ph1 + 2*pi*(t-t1)*f,

where ph1 is the phase at the moment of time t = t1. Then our signal 
could be modeled as

z(t) = exp(j*phase(t)) = exp(j*ph1 + j*2*pi*(t-t1)*f).

And let's also define ph2 as the phase at the moment of time t = t2, i.e.

ph2 = ph1 + 2pi*(t2-t1)*f.

The correlator output (Co) is nothing else but
Co = integral (from t1 to t2) z * dt = exp(j*ph1) * integral (from t1 to 
t2) exp(j*2*pi*(t-t1)*f) * dt = exp(j*ph1) * (1/(j*2*pi*f)) * 
[exp(j*2*pi*(t2-t1)*f) - exp(j*2*pi*(t1-t1)*f)] = (1/(j*2*pi*f)) * 
[exp(j*ph2) - exp(j*ph1) ]

Now it gets a little dirty, because we have to assume that frequency 
offset was small, so that ph2-ph1 << 2*pi. In that case the correlator 
output can be approximated as
Co = (1/(j*2*pi*f)) * exp(j*ph1) * [exp(j*(ph2-ph1)) - 1] = 
(1/(j*2*pi*f)) * exp(j*ph1) * [1 + j*(ph2-ph1) + 0.5*j^2*(ph2-ph1)^2 - 
1] = (1/(j*2*pi*f)) * exp(j*ph1) * j*(ph2-ph1) * [1 + 0.5*j*(ph2-ph1)] = 
(ph2-ph1)/(2*pi*f) * exp(j*ph1) * exp(0.5*j*(ph2-ph1)) = 
(ph2-ph1)/(2*pi*f) * exp(0.5*j*(ph2+ph1))

So the angle of the correlator output is the mean of ph1 and ph2; which 
means the angle of the correlator output is the angle of z at time
tc = (t1+t2)/2. That's because
phase(tc) = ph1 + 2*pi*0.5*(t2-t1)*f = 0.5*ph1 + 0.5*(ph1 + 
2*pi*(t2-t1)*f) = 0.5*(ph1 + ph2).
Which means correlator output angle is the phase of z(t) in the center 
of the respectful time interval (t2; t1).

If I was not convincing enough, you can get that (or slightly different) 
derivation from communications textbooks like "Synchronization 
Techniques for Digital Receivers" by Mengali and D'Andrea.

Gene

On 16.06.2016 11:32, Evgeny Filatov wrote:
> On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:
>> In continuation of the topic I started earlier, which is about the
>> false phase lock in PLL
>> (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew98E).
>>
>> I am trying to implement a false phase lock detector using the
>> approach described in Kaplan book "Understanding GPS: Principles and
>> Applications": http://tinyurl.com/z27gl4v
>>
>> The idea is to compute the phase angle between the in-phase and
>> quadrature pair of samples (I/Q) from the prompt replica correlator
>> within 20 ms data bit interval for L1 C/A GPS signal:
>>
>> C = cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
>> D = dot product of the same vectors
>> angle = atan2(C, D)
>>
>> The angle divided by the time 'T'. over which it was accumulated,
>> should give a frequency:
>>
>> F = angle / (2 * pi * T) [Hz]
>>
>> Note that the time interval 'T' is the time between (Ip-1,Qp-1) and
>> (Ip,Qp) readings and it is slightly less than 20ms. In my case, the
>> (Ip-1,Qp-1) reading is done at ~300us and the (Ip,Qp) reading is done
>> at 20ms after the data bit edge. Therefore, T is ~19.7 ms.
>> The PLL discriminator sees only the sum of two vectors
>> (Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.
>>
>> The computed correction F can be applied to PLL state to correct the
>> false lock error.
>>
>> This approach seems to be working in the sense that, if I apply it to
>> the correlator readings with the PLL in the false lock condition, then
>> it gives me ~12.5 Hz error estimation with a correct sign. At the same
>> time, I know that the true PLL error is 25Hz. BTW, if I do the error
>> computation with the true error being 75Hz instead of 25Hz, the result
>> is still ~12.5Hz.
>>
>> The thing, which I cannot figure out is why the approach gives me only
>> the half of the true error: 12.5 instead of 25 Hz. Any ideas are very
>> much appreciated.
>>
>> BTW, Kaplan talks about the threshold of 15.5Hz and then the 25
>> correction with the proper sign is applied. However my results never
>> reach that threshold.
>>
>
> If I have understood you correctly, you are doing the following thing.
> Assume that one 20ms time interval starts from t_1 = 0 and ends at t_4 =
> 20ms. The vector (Ip-1, Qp-1) is obtained as the correlator output for
> the time interval from t_1 = 0 to t_2 = 300us, and the vector (Ip, Qp)
> is the correlator output for the time interval from t_3 = 300 us to t_4
> = 20 ms. Right?
>
> Now if I got that part right, there's nothing surprising in your result!
> Because the parameter T in your formula for frequency estimation is the
> time interval between CENTERS of the respectful time intervals (t_1;
> t_2) and (t_3; t_4). Which is
>
> T = 0.5*(t_4 - t_3) - 0.5*(t_2 - t_1) = 9.7 ms
>
> rather than the value you assumed (19.7 ms). If you use the correct
> value of T = 9.7 ms, you'd get the correct frequency reading.
>
> So, why the centers of the time intervals, and not something else, like
> ending positions as you have assumed? It helps to understand how do you
> obtain that value. You have two vectors (Ip-1, Qp-1) and (Ip, Qp); and
> in atan2() operation you obtain the angle between these two vectors. It
> means that for purposes of frequency estimation you only care about the
> direction of each vector. What is the direction of a vector? In your
> case, the direction of each vector is the PHASE ESTIMATE obtained within
> the corresponding time interval (note that in your case you obtain phase
> estimate in the form of a phasor, i.e. a complex exponent exp^(j*angle) ).
>
> Now the question is, if you measure phase estimate for the correlation
> interval t=(t1; t2), assuming constant frequency offset f, FOR WHICH
> MOMENT OF TIME would you get a phase estimate? Let's start with the
> system model. Assume phase evolves with time as
>
> phase(t) = ph1 + 2*pi*(t-t1)*f,
>
> where ph1 is the phase at the moment of time t = t1. Then our signal
> could be modeled as
>
> z(t) = exp(j*phase(t)) = exp(j*ph1 + j*2*pi*(t-t1)*f).
>
> And let's also define ph2 as the phase at the moment of time t = t2, i.e.
>
> ph2 = ph1 + 2pi*(t2-t1)*f.
>
> The correlator output (Co) is nothing else but
> Co = integral (from t1 to t2) z * dt = exp(j*ph1) * integral (from t1 to
> t2) exp(j*2*pi*(t-t1)*f) * dt = exp(j*ph1) * (1/(j*2*pi*f)) *
> [exp(j*2*pi*(t2-t1)*f) - exp(j*2*pi*(t1-t1)*f)] = (1/(j*2*pi*f)) *
> [exp(j*ph2) - exp(j*ph1) ]
>
> Now it gets a little dirty,

Of course there's also a simple geometrical interpretation, if you 
visualize these complex phasors as vectors on a plane. Since
Co = [1/(2*pi*f)] * exp(-j*pi/2) * [exp(j*ph2) - exp(j*ph1)],

you can say that Co is a vector orthogonal to the line connecting two 
points, exp(j*ph1) and exp(j*ph2), (which correstond to the starting and 
the ending positions of the time interval). Which, of course, is pointed 
in the direction of the mean phase, exp(j*(ph1+ph2)*0.5).

Also note that this simple interpretation did not involve making any 
approximations, so it's _exact_.

Gene

On 16/06/2016 11:32, Evgeny Filatov wrote:
> On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:
>> In continuation of the topic I started earlier, which is about the
>> false phase lock in PLL
>> (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew98E).
>>
>> I am trying to implement a false phase lock detector using the
>> approach described in Kaplan book "Understanding GPS: Principles and
>> Applications": http://tinyurl.com/z27gl4v
>>
>> The idea is to compute the phase angle between the in-phase and
>> quadrature pair of samples (I/Q) from the prompt replica correlator
>> within 20 ms data bit interval for L1 C/A GPS signal:
>>
>> C = cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
>> D = dot product of the same vectors
>> angle = atan2(C, D)
>>
>> The angle divided by the time 'T'. over which it was accumulated,
>> should give a frequency:
>>
>> F = angle / (2 * pi * T) [Hz]
>>
>> Note that the time interval 'T' is the time between (Ip-1,Qp-1) and
>> (Ip,Qp) readings and it is slightly less than 20ms. In my case, the
>> (Ip-1,Qp-1) reading is done at ~300us and the (Ip,Qp) reading is done
>> at 20ms after the data bit edge. Therefore, T is ~19.7 ms.
>> The PLL discriminator sees only the sum of two vectors
>> (Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.
>>
>> The computed correction F can be applied to PLL state to correct the
>> false lock error.
>>
>> This approach seems to be working in the sense that, if I apply it to
>> the correlator readings with the PLL in the false lock condition, then
>> it gives me ~12.5 Hz error estimation with a correct sign. At the same
>> time, I know that the true PLL error is 25Hz. BTW, if I do the error
>> computation with the true error being 75Hz instead of 25Hz, the result
>> is still ~12.5Hz.
>>
>> The thing, which I cannot figure out is why the approach gives me only
>> the half of the true error: 12.5 instead of 25 Hz. Any ideas are very
>> much appreciated.
>>
>> BTW, Kaplan talks about the threshold of 15.5Hz and then the 25
>> correction with the proper sign is applied. However my results never
>> reach that threshold.
>>
>
> If I have understood you correctly, you are doing the following thing.
> Assume that one 20ms time interval starts from t_1 = 0 and ends at t_4 =
> 20ms. The vector (Ip-1, Qp-1) is obtained as the correlator output for
> the time interval from t_1 = 0 to t_2 = 300us, and the vector (Ip, Qp)
> is the correlator output for the time interval from t_3 = 300 us to t_4
> = 20 ms. Right?
>
> Now if I got that part right, there's nothing surprising in your result!
> Because the parameter T in your formula for frequency estimation is the
> time interval between CENTERS of the respectful time intervals (t_1;
> t_2) and (t_3; t_4). Which is
>
> T = 0.5*(t_4 - t_3) - 0.5*(t_2 - t_1) = 9.7 ms
>
> rather than the value you assumed (19.7 ms). If you use the correct
> value of T = 9.7 ms, you'd get the correct frequency reading.
>
> So, why the centers of the time intervals, and not something else, like
> ending positions as you have assumed? It helps to understand how do you
> obtain that value. You have two vectors (Ip-1, Qp-1) and (Ip, Qp); and
> in atan2() operation you obtain the angle between these two vectors. It
> means that for purposes of frequency estimation you only care about the
> direction of each vector. What is the direction of a vector? In your
> case, the direction of each vector is the PHASE ESTIMATE obtained within
> the corresponding time interval (note that in your case you obtain phase
> estimate in the form of a phasor, i.e. a complex exponent exp^(j*angle) ).
>
> Now the question is, if you measure phase estimate for the correlation
> interval t=(t1; t2), assuming constant frequency offset f, FOR WHICH
> MOMENT OF TIME would you get a phase estimate? Let's start with the
> system model. Assume phase evolves with time as
>
> phase(t) = ph1 + 2*pi*(t-t1)*f,
>
> where ph1 is the phase at the moment of time t = t1. Then our signal
> could be modeled as
>
> z(t) = exp(j*phase(t)) = exp(j*ph1 + j*2*pi*(t-t1)*f).
>
> And let's also define ph2 as the phase at the moment of time t = t2, i.e.
>
> ph2 = ph1 + 2pi*(t2-t1)*f.
>
> The correlator output (Co) is nothing else but
> Co = integral (from t1 to t2) z * dt = exp(j*ph1) * integral (from t1 to
> t2) exp(j*2*pi*(t-t1)*f) * dt = exp(j*ph1) * (1/(j*2*pi*f)) *
> [exp(j*2*pi*(t2-t1)*f) - exp(j*2*pi*(t1-t1)*f)] = (1/(j*2*pi*f)) *
> [exp(j*ph2) - exp(j*ph1) ]
>

Now we just make a substitution ph_c = 0.5*(ph2 + ph1) and ph_d = 
0.5*(ph2 - ph1), to get correlator output

Co = (1/(j*2*pi*f)) * [exp(j*(ph_c + ph_d)) - exp(j*(ph_c - ph_d)) ] = 
(1/(j*2*pi*f)) * exp(j*ph_c) * [exp(j*ph_d) - exp(-j*ph_d)] = 
(1/(j*2*pi*f)) * exp(j*ph_c) * 2*j*sin(ph_d) =
[sin(ph_d) / (pi * f)] * exp(j*ph_c)

I.e. correlator output Co is a complex phasor (with some magnitude that 
we ignore for now) with the phase in the center of the time interval 
(t1; t2).

... sorry for the previous two posts. I'm a bit dumb in the morning ;-)

Gene

On 16/06/2016 13:53, Evgeny Filatov wrote:
> On 16/06/2016 11:32, Evgeny Filatov wrote:
>> On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:
>>> In continuation of the topic I started earlier, which is about the
>>> false phase lock in PLL
>>> (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew98E).
>>>
>>> I am trying to implement a false phase lock detector using the
>>> approach described in Kaplan book "Understanding GPS: Principles and
>>> Applications": http://tinyurl.com/z27gl4v
>>>
>>> The idea is to compute the phase angle between the in-phase and
>>> quadrature pair of samples (I/Q) from the prompt replica correlator
>>> within 20 ms data bit interval for L1 C/A GPS signal:
>>>
>>> C = cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
>>> D = dot product of the same vectors
>>> angle = atan2(C, D)
>>>
>>> The angle divided by the time 'T'. over which it was accumulated,
>>> should give a frequency:
>>>
>>> F = angle / (2 * pi * T) [Hz]
>>>
>>> Note that the time interval 'T' is the time between (Ip-1,Qp-1) and
>>> (Ip,Qp) readings and it is slightly less than 20ms. In my case, the
>>> (Ip-1,Qp-1) reading is done at ~300us and the (Ip,Qp) reading is done
>>> at 20ms after the data bit edge. Therefore, T is ~19.7 ms.
>>> The PLL discriminator sees only the sum of two vectors
>>> (Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.
>>>
>>> The computed correction F can be applied to PLL state to correct the
>>> false lock error.
>>>
>>> This approach seems to be working in the sense that, if I apply it to
>>> the correlator readings with the PLL in the false lock condition, then
>>> it gives me ~12.5 Hz error estimation with a correct sign. At the same
>>> time, I know that the true PLL error is 25Hz. BTW, if I do the error
>>> computation with the true error being 75Hz instead of 25Hz, the result
>>> is still ~12.5Hz.
>>>
>>> The thing, which I cannot figure out is why the approach gives me only
>>> the half of the true error: 12.5 instead of 25 Hz. Any ideas are very
>>> much appreciated.
>>>
>>> BTW, Kaplan talks about the threshold of 15.5Hz and then the 25
>>> correction with the proper sign is applied. However my results never
>>> reach that threshold.
>>>
>>
>> If I have understood you correctly, you are doing the following thing.
>> Assume that one 20ms time interval starts from t_1 = 0 and ends at t_4 =
>> 20ms. The vector (Ip-1, Qp-1) is obtained as the correlator output for
>> the time interval from t_1 = 0 to t_2 = 300us, and the vector (Ip, Qp)
>> is the correlator output for the time interval from t_3 = 300 us to t_4
>> = 20 ms. Right?
>>
>> Now if I got that part right, there's nothing surprising in your result!
>> Because the parameter T in your formula for frequency estimation is the
>> time interval between CENTERS of the respectful time intervals (t_1;
>> t_2) and (t_3; t_4). Which is
>>
>> T = 0.5*(t_4 - t_3) - 0.5*(t_2 - t_1) = 9.7 ms
>>
>> rather than the value you assumed (19.7 ms). If you use the correct
>> value of T = 9.7 ms, you'd get the correct frequency reading.
>>
>> So, why the centers of the time intervals, and not something else, like
>> ending positions as you have assumed? It helps to understand how do you
>> obtain that value. You have two vectors (Ip-1, Qp-1) and (Ip, Qp); and
>> in atan2() operation you obtain the angle between these two vectors. It
>> means that for purposes of frequency estimation you only care about the
>> direction of each vector. What is the direction of a vector? In your
>> case, the direction of each vector is the PHASE ESTIMATE obtained within
>> the corresponding time interval (note that in your case you obtain phase
>> estimate in the form of a phasor, i.e. a complex exponent
>> exp^(j*angle) ).
>>
>> Now the question is, if you measure phase estimate for the correlation
>> interval t=(t1; t2), assuming constant frequency offset f, FOR WHICH
>> MOMENT OF TIME would you get a phase estimate? Let's start with the
>> system model. Assume phase evolves with time as
>>
>> phase(t) = ph1 + 2*pi*(t-t1)*f,
>>
>> where ph1 is the phase at the moment of time t = t1. Then our signal
>> could be modeled as
>>
>> z(t) = exp(j*phase(t)) = exp(j*ph1 + j*2*pi*(t-t1)*f).
>>
>> And let's also define ph2 as the phase at the moment of time t = t2, i.e.
>>
>> ph2 = ph1 + 2pi*(t2-t1)*f.
>>
>> The correlator output (Co) is nothing else but
>> Co = integral (from t1 to t2) z * dt = exp(j*ph1) * integral (from t1 to
>> t2) exp(j*2*pi*(t-t1)*f) * dt = exp(j*ph1) * (1/(j*2*pi*f)) *
>> [exp(j*2*pi*(t2-t1)*f) - exp(j*2*pi*(t1-t1)*f)] = (1/(j*2*pi*f)) *
>> [exp(j*ph2) - exp(j*ph1) ]
>>
>
> Now we just make a substitution ph_c = 0.5*(ph2 + ph1) and ph_d =
> 0.5*(ph2 - ph1), to get correlator output
>
> Co = (1/(j*2*pi*f)) * [exp(j*(ph_c + ph_d)) - exp(j*(ph_c - ph_d)) ] =
> (1/(j*2*pi*f)) * exp(j*ph_c) * [exp(j*ph_d) - exp(-j*ph_d)] =
> (1/(j*2*pi*f)) * exp(j*ph_c) * 2*j*sin(ph_d) =
> [sin(ph_d) / (pi * f)] * exp(j*ph_c)
>
> I.e. correlator output Co is a complex phasor (with some magnitude that
> we ignore for now) with the phase in the center of the time interval
> (t1; t2).
>

Now if we do _not_ ignore the magnitude, we get, since ph_d = pi*f*(t2-t1),

Co = (t2-t1) * sinc(f*(t2-t1)) * exp(j*ph_c),

where sinc(x) = sin(pi*x)/(pi*x).

I don't understand why you are doing that, but you are clearly getting 
low-magnitude correlator outputs for _either_ interval. So if I were you 
I would stick with a couple of equal-length 10 ms intervals to ensure 
robust operation in the low SNR region.

Gene

On 16/06/2016 14:18, Evgeny Filatov wrote:
> On 16/06/2016 13:53, Evgeny Filatov wrote:
>> On 16/06/2016 11:32, Evgeny Filatov wrote:
>>> On 15.06.2016 23:59, adel.mamin.private@gmail.com wrote:
>>>> In continuation of the topic I started earlier, which is about the
>>>> false phase lock in PLL
>>>> (https://groups.google.com/forum/#!topic/comp.dsp/izuEzbew98E).
>>>>
>>>> I am trying to implement a false phase lock detector using the
>>>> approach described in Kaplan book "Understanding GPS: Principles and
>>>> Applications": http://tinyurl.com/z27gl4v
>>>>
>>>> The idea is to compute the phase angle between the in-phase and
>>>> quadrature pair of samples (I/Q) from the prompt replica correlator
>>>> within 20 ms data bit interval for L1 C/A GPS signal:
>>>>
>>>> C = cross product of vectors (Ip-1,Qp-1) and (Ip,Qp)
>>>> D = dot product of the same vectors
>>>> angle = atan2(C, D)
>>>>
>>>> The angle divided by the time 'T'. over which it was accumulated,
>>>> should give a frequency:
>>>>
>>>> F = angle / (2 * pi * T) [Hz]
>>>>
>>>> Note that the time interval 'T' is the time between (Ip-1,Qp-1) and
>>>> (Ip,Qp) readings and it is slightly less than 20ms. In my case, the
>>>> (Ip-1,Qp-1) reading is done at ~300us and the (Ip,Qp) reading is done
>>>> at 20ms after the data bit edge. Therefore, T is ~19.7 ms.
>>>> The PLL discriminator sees only the sum of two vectors
>>>> (Ip-1 + Ip, Qp-1 + Qp), which is the result of full 20ms correlation.
>>>>
>>>> The computed correction F can be applied to PLL state to correct the
>>>> false lock error.
>>>>
>>>> This approach seems to be working in the sense that, if I apply it to
>>>> the correlator readings with the PLL in the false lock condition, then
>>>> it gives me ~12.5 Hz error estimation with a correct sign. At the same
>>>> time, I know that the true PLL error is 25Hz. BTW, if I do the error
>>>> computation with the true error being 75Hz instead of 25Hz, the result
>>>> is still ~12.5Hz.
>>>>
>>>> The thing, which I cannot figure out is why the approach gives me only
>>>> the half of the true error: 12.5 instead of 25 Hz. Any ideas are very
>>>> much appreciated.
>>>>
>>>> BTW, Kaplan talks about the threshold of 15.5Hz and then the 25
>>>> correction with the proper sign is applied. However my results never
>>>> reach that threshold.
>>>>
>>>
>>> If I have understood you correctly, you are doing the following thing.
>>> Assume that one 20ms time interval starts from t_1 = 0 and ends at t_4 =
>>> 20ms. The vector (Ip-1, Qp-1) is obtained as the correlator output for
>>> the time interval from t_1 = 0 to t_2 = 300us, and the vector (Ip, Qp)
>>> is the correlator output for the time interval from t_3 = 300 us to t_4
>>> = 20 ms. Right?
>>>
>>> Now if I got that part right, there's nothing surprising in your result!
>>> Because the parameter T in your formula for frequency estimation is the
>>> time interval between CENTERS of the respectful time intervals (t_1;
>>> t_2) and (t_3; t_4). Which is
>>>
>>> T = 0.5*(t_4 - t_3) - 0.5*(t_2 - t_1) = 9.7 ms
>>>
>>> rather than the value you assumed (19.7 ms). If you use the correct
>>> value of T = 9.7 ms, you'd get the correct frequency reading.
>>>
>>> So, why the centers of the time intervals, and not something else, like
>>> ending positions as you have assumed? It helps to understand how do you
>>> obtain that value. You have two vectors (Ip-1, Qp-1) and (Ip, Qp); and
>>> in atan2() operation you obtain the angle between these two vectors. It
>>> means that for purposes of frequency estimation you only care about the
>>> direction of each vector. What is the direction of a vector? In your
>>> case, the direction of each vector is the PHASE ESTIMATE obtained within
>>> the corresponding time interval (note that in your case you obtain phase
>>> estimate in the form of a phasor, i.e. a complex exponent
>>> exp^(j*angle) ).
>>>
>>> Now the question is, if you measure phase estimate for the correlation
>>> interval t=(t1; t2), assuming constant frequency offset f, FOR WHICH
>>> MOMENT OF TIME would you get a phase estimate? Let's start with the
>>> system model. Assume phase evolves with time as
>>>
>>> phase(t) = ph1 + 2*pi*(t-t1)*f,
>>>
>>> where ph1 is the phase at the moment of time t = t1. Then our signal
>>> could be modeled as
>>>
>>> z(t) = exp(j*phase(t)) = exp(j*ph1 + j*2*pi*(t-t1)*f).
>>>
>>> And let's also define ph2 as the phase at the moment of time t = t2,
>>> i.e.
>>>
>>> ph2 = ph1 + 2pi*(t2-t1)*f.
>>>
>>> The correlator output (Co) is nothing else but
>>> Co = integral (from t1 to t2) z * dt = exp(j*ph1) * integral (from t1 to
>>> t2) exp(j*2*pi*(t-t1)*f) * dt = exp(j*ph1) * (1/(j*2*pi*f)) *
>>> [exp(j*2*pi*(t2-t1)*f) - exp(j*2*pi*(t1-t1)*f)] = (1/(j*2*pi*f)) *
>>> [exp(j*ph2) - exp(j*ph1) ]
>>>
>>
>> Now we just make a substitution ph_c = 0.5*(ph2 + ph1) and ph_d =
>> 0.5*(ph2 - ph1), to get correlator output
>>
>> Co = (1/(j*2*pi*f)) * [exp(j*(ph_c + ph_d)) - exp(j*(ph_c - ph_d)) ] =
>> (1/(j*2*pi*f)) * exp(j*ph_c) * [exp(j*ph_d) - exp(-j*ph_d)] =
>> (1/(j*2*pi*f)) * exp(j*ph_c) * 2*j*sin(ph_d) =
>> [sin(ph_d) / (pi * f)] * exp(j*ph_c)
>>
>> I.e. correlator output Co is a complex phasor (with some magnitude that
>> we ignore for now) with the phase in the center of the time interval
>> (t1; t2).
>>
>
> Now if we do _not_ ignore the magnitude, we get, since ph_d = pi*f*(t2-t1),
>
> Co = (t2-t1) * sinc(f*(t2-t1)) * exp(j*ph_c),
>
> where sinc(x) = sin(pi*x)/(pi*x).
>
> I don't understand why you are doing that, but you are clearly getting
> low-magnitude correlator outputs for _either_ interval. So if I were you
> I would stick with a couple of equal-length 10 ms intervals to ensure
> robust operation in the low SNR region.
>
> Gene
>

I think that formula appears in just _every_ normal spread-spectrum text 
(perhaps in more detailed form, but whatever). Like, for example, Don 
Torrieri's "Principles of Spread-Spectrum Communication Systems". I'm 
not sure Kaplan covers that stuff.

Gene

Gene, this explanation is beautiful! All the assumptions you've made are co=
rrect and match my case exactly. I have followed you all along and everythi=
ng you've written makes sense to me. Thank you very much for shedding light=
 on it. I really appreciate it.

About the 300us + 19.7ms split. This came from the HW correlator implementa=
tion I am dealing with at the moment. It does not have sample buffering cap=
abilities and the tracking loop corrections are taken in use only 300us aft=
er the beginning of the correlation round. So, I was going to take advantag=
e of it and use the intermediate 300us I/Q reading also for the false lock =
detector. Now I am convinced that I have to reconsider it and do an extra I=
/Q reading after 9.7ms since the last 300us reading just for the sake of th=
e false lock detector. Luckily the HW allows for that.

FWIW, it turns out that the Kaplan book has a diagram typo, where he explai=
ns the false frequency and false phase lock detector details: http://tinyur=
l.com/hz4pp37. The diagram there suggests that you should accumulate DOT an=
d CROSS products over 1 second and then compute atan2(CROSS, DOT). However,=
 since atan2's range is -pi,+pi, the overall algorithm will never result in=
 anything more than pi/(2*pi) =3D 0.5 Hz error. To fix this, the accumulati=
on step should be done after atan2 is computed.

On 16/06/2016 15:55, adel@swift-nav.com wrote:
> Gene, this explanation is beautiful! All the assumptions you've made are correct and match my case exactly. I have followed you all along and everything you've written makes sense to me. Thank you very much for shedding light on it. I really appreciate it.
>
> About the 300us + 19.7ms split. This came from the HW correlator implementation I am dealing with at the moment. It does not have sample buffering capabilities and the tracking loop corrections are taken in use only 300us after the beginning of the correlation round. So, I was going to take advantage of it and use the intermediate 300us I/Q reading also for the false lock detector. Now I am convinced that I have to reconsider it and do an extra I/Q reading after 9.7ms since the last 300us reading just for the sake of the false lock detector. Luckily the HW allows for that.
>

Adel, you are welcome!

> FWIW, it turns out that the Kaplan book has a diagram typo, where he explains the false frequency and false phase lock detector details: http://tinyurl.com/hz4pp37. The diagram there suggests that you should accumulate DOT and CROSS products over 1 second and then compute atan2(CROSS, DOT). However, since atan2's range is -pi,+pi, the overall algorithm will never result in anything more than pi/(2*pi) = 0.5 Hz error. To fix this, the accumulation step should be done after atan2 is computed.
>

I think Kaplan should be trusted here.

You can accumulate DOT and CROSS products, because if you have a couple 
of vectors 'a' and 'b', and 'angle' is angle between them, then
DOT   = abs(a) * abs(b) * cos(angle),
CROSS = abs(a) * abs(b) * sin(angle),
where abs(a) denotes magnitude of vector a.
(You can see that DOT and CROSS are actually just scalar and vector 
multiplication, respectfully.)

Accumulation is allowed at this point, because 'angle' already contains 
frequency information, as angle = 2*pi*f*delta_t, where f is frequency 
and delta_t is the time difference between centers of a pair of 
intervals for which you calculate DOT and CROSS.

Then atan2 only retrieves the angle from a point on a circle, defined by 
its Cartesian coordinates cos(angle) and sin(angle).

Gene