NCO and Costas Loop

Hello All,

I have been working in C++ on a Costas Loop as a the basis for a carrier
recovery for a BPSK decoder. My background is not DSP but I have been doing
quite a  bit of research trying to come up to speed. This project is not an
academic or commercial project, it's just a personal project.

Anyway, I believe that I need to use a second order loop because I want to
track both frequency and phase offsets. In my research, I found multiple
sources that suggest the following can be used to calculate the phase and
frequency offsets for an NCO:

      d_freq = d_freq + d_beta * error;
      d_phase = d_phase + d_freq + d_alpha * error;

I also have found the following relationship (as a starting point) for
alpha and beta:

alpha = ~0.01
betas = alpha^2/4

I am currently calculating error as:

error = atan2(sample.q, sample.i)

And sample.q and sample.i have been properly filters to remove the 2Wc
component (using FIR filters).

This all seems fairly simply to plug into my loop but as you might have
guessed it doesn't track, it just goes off into the weeds.

I am not currently filtering the error signal. I understand that might be
necessary.

The I and Q are generated with the following:

Inptr1->x = pIn[i] * cos (d_phase);
Inptr1->y = pIn[i] * -sin (d_phase);

Inptr1 passes thru a FIR to generate "sample.x and sample.y).

My problem is two things. While from a theoretical perspective the
equations for d_phase and d_freq make complete sense to me, I don't think I
have them implemented correctly for a discrete time implementation. I am
hoping to get some practical guidance on how to incorporate frequency and
phase offset into an NCO in a discrete time implementation of a Costas
Loop. My second issue is that I have not been able to understand how to
design a simple IIR to be used as a filter for the error signal - any
practical pointers to sources to understand this are welcome.

thanks,

Joe

On Fri, 08 May 2009 10:29:25 -0500, rfnoise wrote:

> Hello All,
> 
> I have been working in C++ on a Costas Loop as a the basis for a carrier
> recovery for a BPSK decoder. My background is not DSP but I have been
> doing quite a  bit of research trying to come up to speed. This project
> is not an academic or commercial project, it's just a personal project.
> 
> Anyway, I believe that I need to use a second order loop because I want
> to track both frequency and phase offsets. In my research, I found
> multiple sources that suggest the following can be used to calculate the
> phase and frequency offsets for an NCO:
> 
>       d_freq = d_freq + d_beta * error;
>       d_phase = d_phase + d_freq + d_alpha * error;
> 
> I also have found the following relationship (as a starting point) for
> alpha and beta:
> 
> alpha = ~0.01
> betas = alpha^2/4
> 
> I am currently calculating error as:
> 
> error = atan2(sample.q, sample.i)
> 
> And sample.q and sample.i have been properly filters to remove the 2Wc
> component (using FIR filters).
> 
> This all seems fairly simply to plug into my loop but as you might have
> guessed it doesn't track, it just goes off into the weeds.
> 
> I am not currently filtering the error signal. I understand that might
> be necessary.
> 
> The I and Q are generated with the following:
> 
> Inptr1->x = pIn[i] * cos (d_phase);
> Inptr1->y = pIn[i] * -sin (d_phase);
> 
> Inptr1 passes thru a FIR to generate "sample.x and sample.y).
> 
> My problem is two things. While from a theoretical perspective the
> equations for d_phase and d_freq make complete sense to me, I don't
> think I have them implemented correctly for a discrete time
> implementation. I am hoping to get some practical guidance on how to
> incorporate frequency and phase offset into an NCO in a discrete time
> implementation of a Costas Loop. My second issue is that I have not been
> able to understand how to design a simple IIR to be used as a filter for
> the error signal - any practical pointers to sources to understand this
> are welcome.
> 
> thanks,
> 
> Joe

Hey Joe:

You're building a closed-loop control system, that's made more difficult 
by that nasty old nonlinear phase comparator in there.

So you may want to start by taking the phase comparator out, and just 
seeing if your proposed loop servos when you use your value for phase 
instead of going through the whole thing.

If that works, then see if your loop isn't trying to servo on the spot 
where the phase comparator rolls over from -pi to +pi, instead of 
settling to the (nicely linear) zero point.

If it doesn't work, then ponder on the fact that you're designing a PI 
controller with an integrating plant, that beta is your integrator gain 
and that alpha is your proportional gain.  Then read this article: http://
www.wescottdesign.com/articles/Sampling/pidwophd.html and see if it 
doesn't help you find the right values for your application.

-- 
http://www.wescottdesign.com

Hi Tim,

Thanks for the response.

>On Fri, 08 May 2009 10:29:25 -0500, rfnoise wrote:
>
>So you may want to start by taking the phase comparator out, and just 
>seeing if your proposed loop servos when you use your value for phase 
>instead of going through the whole thing.

I tried that and the loop seems to respond as expected with manual values
of phase.

>If that works, then see if your loop isn't trying to servo on the spot 
>where the phase comparator rolls over from -pi to +pi, instead of 
>settling to the (nicely linear) zero point.

When I hook the phase detector back up and examine the error signal
periodically across the sample space I see values typically less than 2 and
greater then -2. Some values slightly above and below. The error value
never reaches zero, it's gets close (0.06) occasionally within sample
space.

>If it doesn't work, then ponder on the fact that you're designing a PI 
>controller with an integrating plant, that beta is your integrator gain 
>and that alpha is your proportional gain.  Then read this article:
http://
>www.wescottdesign.com/articles/Sampling/pidwophd.html and see if it 
>doesn't help you find the right values for your application.

Thanks. Very well written and understandable. 

My BETA must be very small (and smaller than ALPHA) before ALPHA has any
affect. By adjusting the values following your article I see mostly what I
would call oscillation - that is, never finding the correct frequency and
phase (hunting). As I get close to my original values (ALPHA = 0.00001 and
BETA = ALPHA*ALPHA/4.0) I see that sigma goes to 0 but unfortunately the
error value never settles in close to zero.

These observations are all with a real BPSK signal where I know the
approximate center frequency (which is what I seed the loop with). If I
switch over to a basic sine wave (same center frequency) the loop behaves
better in that error signal goes much much closer to zero, sigma goes to
zero and the loop tracks small changes in frequency offsets in the sine
wave input. Also, the Q values go to zero, which I think indicates there is
no phase rotation - I am locked.

Perhaps the lack of proper filtering on the error signal is the culprit
with the real BPSK signal?

Any other tips are appreciated.

Thanks,

Joe (wishing noncoherent recovery would be feasible)

On Sat, 09 May 2009 18:07:03 -0500, rfnoise wrote:

> Hi Tim,
> 
> Thanks for the response.
> 
>>On Fri, 08 May 2009 10:29:25 -0500, rfnoise wrote:
>>
>>So you may want to start by taking the phase comparator out, and just
>>seeing if your proposed loop servos when you use your value for phase
>>instead of going through the whole thing.
> 
> I tried that and the loop seems to respond as expected with manual
> values of phase.
> 
>>If that works, then see if your loop isn't trying to servo on the spot
>>where the phase comparator rolls over from -pi to +pi, instead of
>>settling to the (nicely linear) zero point.
> 
> When I hook the phase detector back up and examine the error signal
> periodically across the sample space I see values typically less than 2
> and greater then -2. Some values slightly above and below. The error
> value never reaches zero, it's gets close (0.06) occasionally within
> sample space.
> 
>>If it doesn't work, then ponder on the fact that you're designing a PI
>>controller with an integrating plant, that beta is your integrator gain
>>and that alpha is your proportional gain.  Then read this article:
> http://
>>www.wescottdesign.com/articles/Sampling/pidwophd.html and see if it
>>doesn't help you find the right values for your application.
> 
> Thanks. Very well written and understandable.
> 
> My BETA must be very small (and smaller than ALPHA) before ALPHA has any
> affect. By adjusting the values following your article I see mostly what
> I would call oscillation - that is, never finding the correct frequency
> and phase (hunting). As I get close to my original values (ALPHA =
> 0.00001 and BETA = ALPHA*ALPHA/4.0) I see that sigma goes to 0 but
> unfortunately the error value never settles in close to zero.
> 
> These observations are all with a real BPSK signal where I know the
> approximate center frequency (which is what I seed the loop with). If I
> switch over to a basic sine wave (same center frequency) the loop
> behaves better in that error signal goes much much closer to zero, sigma
> goes to zero and the loop tracks small changes in frequency offsets in
> the sine wave input. Also, the Q values go to zero, which I think
> indicates there is no phase rotation - I am locked.
> 
> Perhaps the lack of proper filtering on the error signal is the culprit
> with the real BPSK signal?
> 
> Any other tips are appreciated.

Just to review, with a Costas loop you integrate the phase error over one 
bit period, then multiply it by +1 or -1 depending on the value that you 
decided for that bit.  You _are_ doing this, right?  Note that this 
requires an "inphase" signal (i.e. cos(phase)) for bit detection and a 
"quadrature" signal (i.e. sin(phase)) for carrier sync.

An easier loop to implement, assuming you're bringing in sinusoidal data, 
is one where you square the incoming signal (or multiply it by it's 
absolute value), run it through a DC-block filter, then phase lock to the 
resulting 2nd harmonic of your carrier.

Both of these loops have more or less the same performance potential, but 
one may be easier than the other to implement for your particular 
receiver.

-- 
http://www.wescottdesign.com

Tim, thanks again for responding. My comments are below.

>On Sat, 09 May 2009 18:07:03 -0500, rfnoise wrote:
>Just to review, with a Costas loop you integrate the phase error over one

>bit period, then multiply it by +1 or -1 depending on the value that you

>decided for that bit.  You _are_ doing this, right?  Note that this 
>requires an "inphase" signal (i.e. cos(phase)) for bit detection and a 
>"quadrature" signal (i.e. sin(phase)) for carrier sync.

Sorry, Tim, i'm afraid I am not following you on the multiplication of +1
or -1. Are you referring to data recovery? I am not dealing with that yet,
I'm just trying to recover the carrier (however, I believe I will have the
bit stream available in the I arm output once the loop is locked). Perhaps
you are referring to the 180 degree phase ambiguity that I will need to
deal with once I get carrier recovered? Then I will multiply by +1 or -1
once I recognize the particular unique word in the bit stream (or it's
complement). So, it's possible I am completely off track here.... Let me
give you the background that I am working from:

The best description of a costas loop that I have is the following
articles:

http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010315F1.pdf
http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010531F1.pdf
http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010419F1.pdf
http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010628F1.pdf

I believe I have implemented a costas loop as described in the article
above.

I also have referenced this article:

http://rfdesign.com/images/archive/0102Feigin20.pdf

In my source BPSK (suppressed carrier) signal there is some carrier
frequency uncertainty and the usual phase offset to deal with. So, I
believe based on this that I needed to synchronize both frequency and phase
that I need a second order loop (proportional-integration loop). Here is my
specific equations for the controls:

	sigma += (BETA * error);
        loopfilter = sigma + (ALPHA * error);
	phi += PI2*SeedFreq/Fs + loopfilter;

The phase detector is simply:

	phase_detector = atan2(acc.y,acc.x); 

        error = phase_detector; // no filtering applied yet....

acc.y is the imaginary and acc.x is the real output from the NCO mixer
post FIR filtering (3dB half symbol width).

Complex sample mixing is done like this:

	Inptr1->x = pIn[i] *  cos( phi );	
	Inptr1->y = pIn[i] * -sin( phi );

Inptr1 is fed to FIR filters in the I and Q arms.

All this is in a for loop (i) that steps thru the current sample. 

I am monitoring the input, post FIR filtering and error spectra using FFT
display and spitting out phase, sigma, error, Q, I across the sample space
to see what's going on.

No doubt I am missing something critical. I am happy to follow any
specific pointers to papers/books etc. to get me going in the right
direction. 

>
>An easier loop to implement, assuming you're bringing in sinusoidal data,

>is one where you square the incoming signal (or multiply it by it's 
>absolute value), run it through a DC-block filter, then phase lock to the

>resulting 2nd harmonic of your carrier.

Yes, I did run across this method in my research. And I may end up doing
this if I can figure out the costas loop - seems like costas work should
work and I hate to give up simply because I don't understand it yet. 

>Both of these loops have more or less the same performance potential, but

>one may be easier than the other to implement for your particular 
>receiver.
>-- 
>http://www.wescottdesign.com
>

Joe

On Sun, 10 May 2009 00:45:59 -0500, rfnoise wrote:

> Tim, thanks again for responding. My comments are below.
> 
>>On Sat, 09 May 2009 18:07:03 -0500, rfnoise wrote: Just to review, with
>>a Costas loop you integrate the phase error over one
> 
>>bit period, then multiply it by +1 or -1 depending on the value that you
> 
>>decided for that bit.  You _are_ doing this, right?  Note that this
>>requires an "inphase" signal (i.e. cos(phase)) for bit detection and a
>>"quadrature" signal (i.e. sin(phase)) for carrier sync.
> 
> Sorry, Tim, i'm afraid I am not following you on the multiplication of
> +1 or -1. Are you referring to data recovery? I am not dealing with that
> yet, I'm just trying to recover the carrier (however, I believe I will
> have the bit stream available in the I arm output once the loop is
> locked). Perhaps you are referring to the 180 degree phase ambiguity
> that I will need to deal with once I get carrier recovered? Then I will
> multiply by +1 or -1 once I recognize the particular unique word in the
> bit stream (or it's complement). So, it's possible I am completely off
> track here.... Let me give you the background that I am working from:
> 
> The best description of a costas loop that I have is the following
> articles:
> 
> http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010315F1.pdf
> http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010531F1.pdf
> http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010419F1.pdf
> http://i.cmpnet.com/chipcenter/dsp/images/dspsourced/DSP010628F1.pdf
> 
> I believe I have implemented a costas loop as described in the article
> above.
> 
> I also have referenced this article:
> 
> http://rfdesign.com/images/archive/0102Feigin20.pdf
> 
-- snip --

> acc.y is the imaginary and acc.x is the real output from the NCO mixer
> post FIR filtering (3dB half symbol width).
> 
-- snip --

Based on your description of the problem and the fact that just the 
control loop part works correctly, I think you need to focus on the phase 
detector.  Now that you've tried just the loop, try feeding the loop with 
a constant tone; think of it as BPSK with all ones or all zeros.  See if 
that works.  If it does, then try just one phase shift and see how the 
loop behaves -- if it's working right there should be no big transient 
during the change.

One thing I can say from experience is that when you're commissioning a 
system like this, you want to be able to break it down to bits and figure 
out how to put it together one bit at a time, verifying functionality as 
you go -- otherwise you'll be in for a lot of flailing around trying to 
get things to work.

-- 
http://www.wescottdesign.com

>On Sun, 10 May 2009 00:45:59 -0500, rfnoise wrote:
>
>Based on your description of the problem and the fact that just the 
>control loop part works correctly, I think you need to focus on the phase

>detector.  Now that you've tried just the loop, try feeding the loop with

>a constant tone; think of it as BPSK with all ones or all zeros.  See if

>that works.  If it does, then try just one phase shift and see how the 
>loop behaves -- if it's working right there should be no big transient 
>during the change.

Sounds good. I re-tuned with single tone and was able to speed up tracking
and precision. I also added a hard limiter on the output of the phase
detector (this seemed to help quite a bit).

A couple of observations that would be nice to have comment on are the
following:

When the loop is in lock it has been noted in the literature that most of
the energy is in the I path and very little energy is in the Q path. With a
single tone well above the noise floor I see this as an always positive I
values and very close to zero Q values. I see the error value very near
zero. When I increase the noise I see the Q values come up quite a bit - I
am assuming this is normal, is this correct?

When I turned on my real BPSK signal I did not see what I would call a
lock condition. Perhaps a false lock? The loop servo'd to a point and then
remained at that position, however, error signal was not near zero, and
clearly there was rotation occurring. Why does the loop not servo the real
lock point? I can move the seed frequency and find another false lock point
- they seem to be all over the place in my implementation. 

I changed my atan2(q, i) to atan(q/i) and this caused my loop when exposed
to the real BPSK signal to "lock" onto the real carrier frequency - at
least it appears that way. I has much of the energy, error is near zero,
etc. I would call it a rough lock. Now, I know the difference between atan
and atan2 but the literature is conflicting - i see both being used.  

Taming this loop is much like herding cats....

>One thing I can say from experience is that when you're commissioning a 
>system like this, you want to be able to break it down to bits and figure

>out how to put it together one bit at a time, verifying functionality as

>you go -- otherwise you'll be in for a lot of flailing around trying to 
>get things to work.

What I need is the ability to have both frequency and time display of
various points in the system. Maybe I need to bite the bullet and get
MatLab or the simulation package.

Joe

Hello joe

i am trying to make costas loop and i guess i followed the same procedure
and articles as mentioned by you . 
i did get similar problem , the error signal never goes zero, but
fluctuates between -3 and 3 . 
Can you tell me , how you were able to solve the problem . 

Regards 
Pasa

On Fri, 23 Oct 2009 07:18:40 -0500, pasa wrote:

> Hello joe
> 
> i am trying to make costas loop and i guess i followed the same
> procedure and articles as mentioned by you .
> i did get similar problem , the error signal never goes zero, but
> fluctuates between -3 and 3 .
> Can you tell me , how you were able to solve the problem .
> 
> Regards
> Pasa

Context?  This is USENET.

In general, you find the problem and fix it.  In specific -- you haven't 
given enough detail.  If the range of your error signal is +/- 2^15, then 
+/- 3 is pretty good.

-- 
www.wescottdesign.com

Hello , 
thank you for your reply . 

I used  

fc=150;  /* signal frequency  in Hz */
fs=8000;  /* sampling frequency  */
n=0:N-1;   /* N =  2048    */ 
x1=sin(2*pi*(fc/fs)*n); /* input signal generation */

This input signal is multiplied with cosine signal and  sine signal from
NCO . 


which gives signal as given by formula 
    
    v1 = 1/2  A(t) cos(phi - theta)+ 1/2 A(t) cos(4*pi*fc*t+phi+theta)
    v2 = 1/2  A(t) sin(phi - theta)+ 1/2 A(t) sin(4*pi*fc*t+phi+theta)

where phi is phase of received signal and theta is phase of signal from
nco 



So inoder to remove 2*fc frequency of the resulting signal i used low pass
FIR filter 
and When i checked the output signal with N point fft .the signal with
2*fc frequency is attenuated to magnitude of 
12  from 195. and zero frequency magnitude  to 54 from 55. It seems that
filter is working fine. 


The output from filters signal.x and signal.y are then multiplied to give
phase difference (error ) signal, which is used to 
drive NCO .The signal from the NCO have 150 Hz (starting frequency)
i havenot used loop filter  betwen NCO and the output from multiplier.
For NCO i used the same formula as Joe used. and i checked NCO seperately
changing the input signal to NCO and checked the output frequency of the
NCO. 
The NCO seems working fine.  
 

When i plot the error signal i found that the error signal starts from
amplitude 1.74 and jumps to -2.89 and then increases 
to 1.74 again very quickly and forms a ripple between max value of 1.74
and minimum value of 1.62. 


I checked the  frequency of signal from the NCO and its frequency  = 150
Hz.


I repeated the simulation changing the input signal frequency to 180 Hz. 

When i plot the error signal i found that the error signal starts from
1.74 and jumps to -2.87 and then increases 
to 1.28 again very quickly and decreases to -3.125 then it forms a saw
tooth with maximum value of 3.121 and minimum of -3.125
and the NCO frequency is still 150 Hz. 

I must have missed something. 
Please correct me if i wrong . I was hoping that this error signal should
tend towards zero 

Any help is appreciated. 


Regards
Pasa


>On Fri, 23 Oct 2009 07:18:40 -0500, pasa wrote:
>
>> Hello joe
>> 
>> i am trying to make costas loop and i guess i followed the same
>> procedure and articles as mentioned by you .
>> i did get similar problem , the error signal never goes zero, but
>> fluctuates between -3 and 3 .
>> Can you tell me , how you were able to solve the problem .
>> 
>> Regards
>> Pasa
>
>Context?  This is USENET.
>
>In general, you find the problem and fix it.  In specific -- you haven't

>given enough detail.  If the range of your error signal is +/- 2^15, then

>+/- 3 is pretty good.
>
>-- 
>www.wescottdesign.com
>