Costas Loop oscillates in ModelSim while locks in MATLAB!

[Altera_Forum]

Hello again, 

I have designed a Costas Loop in MATLAB and it simulates very well... locking and recovering the carrier. Here are the results:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10456&stc=1   I ported the design to FPGA in VHDL, and simulated in ModelSim... however, the loop has a kinda oscillating behaviour... take a look:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10457&stc=1   This is the Loop Filter output:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10458&stc=1   Why is it oscillating? Is it the loop filter?

[Altera_Forum]

here is copy: 

  %NCO bit true model clear all; close all; %example NCO parameters n = 20000; %number of test samples A = 2^15 - 1; %max amplitude, amplitude resolution 16 bit signed M = 2^32; %NCO accumulator phase resolution inc = round(M*.001); %NCO accumulator phase increment k = 2^12; %lut phase resolution (lut size) lsb = log2(M) - log2(k); %LSBs discarded when addressing lut = round(A*exp(j*2*pi*(0:k-1)/k)); %lut, one cycle sine data ptr = 0; addr = 0; for i = 1:n y(i) = lut(addr+1); %add 1 for matlab LUT ptr = mod(ptr + inc, M); %phase accumulator(ramp) addr = round(ptr/2^lsb); %discard LSBs addr(addr >= k) = addr - k; %check address overflow end    "inc" can be chosen as Fo/Fs scaled by 2^32 but note there is bug for negative frequency, use positive increment and invert at output or avoid matlab mod function. Note also this nco is based on fixed lut resolution (no intermediate values calculated).

[Altera_Forum]

No I don't want to revert to matlab... as you said, modelsim has proper loop and delays. I am gonna focus on it. I was just making a comparison in behaviour.... But is my analysis right? 

Thanks for the code... I might need it.

[Altera_Forum]

That is good path since modelling in matlab for feedback loops needs to be dead accurate match of rtl. Yes your analysis is right if you mean your matlab model which I haven't seen its final format. The nco statement has control over phase just like rtl and this phase is used to control frequency by different jumps over sin/cos functions (luts) so it will have very high floating resolution. but you need to write as y = exp(j*2*pi*linspace(0+offset:N-1+offset,N)*.1); in order to step phase in fractions.

[Altera_Forum]

Come on Mr Kazem! The model that we spent days discussing it! We finally said yuppie remember?! Anyway here its is: 

  % Siraj Muhammad % 25/3/2015 % BPSK Demodulator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load RRC.mat fc = 0.05000001; phase_offset = pi/7; N = length(rf_signal); bb = zeros(1, N); bb_f = zeros(1, N); I_f = zeros(1, N); Q_f = zeros(1, N); I_r = zeros(1, N); Q_r = zeros(1, N); error = zeros(1, N); PhErr = zeros(1, N); error_integral = zeros(1, N); % Loop Filter Coefficients % BW = 200; % Hz % loop_theta = 2*pi*BW; % C1 = 4*(loop_theta)^2/(1+sqrt(2)*loop_theta+loop_theta^2); % C2 = 2*sqrt(2)*loop_theta/(1+sqrt(2)*loop_theta+loop_theta^2); C1 = 2^2; C2 = 1/2^8; for i = 2:N % Downconverting to Baseband bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1)); % Filtering bb_f = filter(RRC, bb(1:i)); I_f = real(bb_f); Q_f = imag(bb_f); % Error error(i) = I_f(i).*Q_f(i); % Loop Filter error_integral(i) = error(i).*C1 + error_integral(i-1); PhErr(i) = error(i).*C2 + error_integral(i); end figure; subplot 321; plot(real(bb)); title('I Channel') subplot 322; plot(imag(bb)); title('Q Channel') subplot 323; plot(I_f); title('I Channel Filter') subplot 324; plot(Q_f); title('Q Channel Filter') subplot 325; plot(error); title('Phase Error'); subplot 326; plot(PhErr); title('Loop Filter');   OK, that's great... let's stick to modelsim now. so I should expect oscillations in the loop filter output. I believe convergence time may be long and ModelSim is going to consume even longer time to complete simulation. I will carry on tests tomorrow and let you know... thanks for help Mr Kazem, I am so grateful to you.

[Altera_Forum]

The above code should prove the concept. 

 The feedback point is at this equation below  bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1));  The expression exp(-j*err) means rotate phase 90 degrees *err clockwise while frequency is unchanged.  This may suit software based dsp but in rtl the nco is controlled via inc value which controls both phase and frequency.  So keep on modelsim

[Altera_Forum]

--- Quote Start ---  The above code should prove the concept.  The feedback point is at this equation below  bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1));  The expression exp(-j*err) means rotate phase 90 degrees *err clockwise while frequency is unchanged.  This may suit software based dsp but in rtl the nco is controlled via inc value which controls both phase and frequency.  So keep on modelsim  --- Quote End ---    Wonderful! This qualifies my previous analysis.  Here is the result of a test with the loop activated and a small phase offset. Output from loop filter is shifted right by 7 and fed to NCO word.  https://www.alteraforum.com/forum/attachment.php?attachmentid=10467   https://www.alteraforum.com/forum/attachment.php?attachmentid=10468   As I shift the output with a smaller number i.e. as I divide by a smaller number, meanders increase in the specified period of time (400 microseconds), but it keeps rising. I don't know if this will lock and settle. What do you think? I am thinking of programming the FPGA with this design and capture signals using SignalTap, if it locks eventually, I should see the same input pattern... do you recommend doing that?

[Altera_Forum]

I suggest you try: 

1) single clean tone e.g. @ 0.01Fs by just setting your symbols to +1 only. This should lock first.  2)Next try single tone but add to it +1/-1 jitter (representing +1,-1,+1,-1 train of pulses) by setting symbols to +1/-1 Does it still lock?  3)Then move to full random symbols.  I have doubts that fully random BPSK will lock easily. You may need to extract a tone eventually through bandpass filter and target it for lock but let us see.

[Altera_Forum]

I am doing some tests... I will do yours as well and let you know.

[Altera_Forum]

Hello Mr Kazem, 

I did the test you suggested... maybe results are disappointing. These are all with a single tone input.  With a small amount of scaling... loop takes some time:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10470&stc=1   Without scaling:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10471&stc=1  http://www.alteraforum.com/forum/attachment.php?attachmentid=10472&stc=1   I don't know if the latter gonna lock eventually. What to do now?  By the way, I found a bug in my code... the loop filter is clocked by a 50 MHz, and its input is the phase detector (multiplier) which in turn has RRCs as input, the RRCs produce a sample every several clock cycles of the 50 MHz... I thought it would be the problem. I made an enable port for the Loop Filter so it only works when a new sample comes out of the RRCs. However, nothing changed... but I think this is how it should be.  ** UPDATE ** This is the output of the loop filter for 1 ms. http://www.alteraforum.com/forum/attachment.php?attachmentid=10473&stc=1   I and Q are sinusoidal signals!! This never gonna lock!! :'(

[Altera_Forum]

when I suggested testing the tx/rx chain with the feedback loop deactivated I wanted a sanity check of your chain. Ideally you should have testbench for every module by reading reference outputs from matlab into testbench and check for bit true comparison. That is what we do in the industry. Moreover you need to scale the outputs of each module to some practical resolution e.g 16 bits signed yet I see your dynamic range too high. 

 Repeat the chain sanity test by checking that what you send from Tx as random symbols is recovered by Rx and check that at symbol rate. I mean your overall rx chain should decimate until reaching symbol rate then pass it through a simple slicer and check pattern.  I did some experiments on this design in matlab running 1 Msamples and it locks for all cases nicely (symbols at +1 all = single tone, symbols at +1/-1 = single tone with jitter and random data. Though I used complex I/Q RF to avoid extra filter (for speed). So I am convinced it should lock eventually). It does not lock if I changed input to QPSK or 16QAM.  The main point I noticed was that I have to update the NCO word once only every say 10 error updates (ignore 9) and revert back to initial nco word. The scaling is essential to match your phase resolution word of nco. You will get a nice oscillation that converges to perfect result or so (I got zero jitter at .01Fs). I did not use proportional term.  Remember what aplies to software dsp does not apply to FPGA nco. DSP engineers do not know about fpga and are only orientated to their platform. Our nco has one control over phase and frequency (nco increment value).

[Altera_Forum]

I am doing the sanity test now. 

 Let me inquire about some points here:   --- Quote Start ---   I did some experiments on this design in matlab running 1 Msamples and it locks for all cases nicely (symbols at +1 all = single tone, symbols at +1/-1 = single tone with jitter and random data. Though I used complex I/Q RF to avoid extra filter (for speed). So I am convinced it should lock eventually). It does not lock if I changed input to QPSK or 16QAM.  The main point I noticed was that I have to update the NCO word once only every say 10 error updates (ignore 9) and revert back to initial nco word. The scaling is essential to match your phase resolution word of nco. You will get a nice oscillation that converges to perfect result or so (I got zero jitter at .01Fs). I did not use proportional term.  Remember what aplies to software dsp does not apply to FPGA nco. DSP engineers do not know about fpga and are only orientated to their platform. Our nco has one control over phase and frequency (nco increment value).  --- Quote End ---    1) You mean the Matlab design that I sent you? 2) Complex I/Q RF mean you used exp() instead of cos/sin functions and separate filters for channels right? 3) Regarding this point:   --- Quote Start ---   The main point I noticed was that I have to update the NCO word once only every say 10 error updates (ignore 9) and revert back to initial nco word. The scaling is essential to match your phase resolution word of nco.   --- Quote End ---   You did that in Matlab or in ModelSim? 4) You are right about proportional term... it seemed useless to me too.

[Altera_Forum]

Here is the chain sanity test result: 

 TX: https://www.alteraforum.com/forum/attachment.php?attachmentid=10474   RX: https://www.alteraforum.com/forum/attachment.php?attachmentid=10475   I believe the chain passes the test. So, problem is in the loop. Well... to be honest I still doubt in the NCO the problem in my point of view is in the nature of the NCO it self because as we discussed.. we are unable to control phase independently from frequency... but your point about using loop filter every 10 samples stirred some hope in me. So, what's next?

[Altera_Forum]

I did a new matlab model (private) using rtl model of nco (at tx and rx) as well as bit true modelling (except for loop filter which Ikept normalised for now).  

Input: variety of symbols (BPSK,QAM) Tx: RRC upsampling by 2 to Rs*2, then mixing with 0.01Fs and directly taking I&Q of mixer as RF (this is just to avoid having high sideband and hence no need for filter) Rx: nco and loop filter. Thus the loop works at Rs*2 rate. I update the nco word once every 5, or 10,50,100 samples reverting back to initial nco word. Hence update rate is Rs*2/5...etc. The loop filter set to 0.0001 for alpha but other values work.  I then use RRC decimating by 2 down to Rs after loop (RRC not inside loop).  It locks for clean or jittered single tone for all cases (BPSK,QPSK,16QAM) but works only for BPSK random symbols (may be it needs more research to work for all but I understand from literature it shouldn't work for QPSK and other QAMs efficiently unless aided by slicing results)

[Altera_Forum]

1) Did you use multiplier as phase detector? 

2) You didn't need RRC inside the loop because you took I/Q directly, right? 3) Do I update the nco word as you did? once every 10 samples? How do you revert back to initial word? I mean do you update nco word with error sample for one clock cycle only?

[Altera_Forum]

Yes I used same error as yours (real * imag) 

Any filter for high sideband will do as long as it does the job but I kept loop at Rs*2 rate (not at symbol rate). Updating nco word conditionaly is like this: use sample counter e.g. 0~9 (not clock counter) if counter = 9 update with error else assign initial value (do not latch error update as it will lose frequency lock)

[Altera_Forum]

1) Did you keep loop at Rs*2 because this is the rate of the input signal to loop filter (sample rate)? 

2) Updating with error is done by summing the error to the initial value or just feeding error to NCO?

[Altera_Forum]

--- Quote Start ---  1) Did you keep loop at Rs*2 because this is the rate of the input signal to loop filter (sample rate)? 2) Updating with error is done by summing the error to the initial value or just feeding error to NCO?  --- Quote End ---    I kept loop at sample rate of Rs*2 based on my past experience. Higher rate should be better in theory, loop @ symbol rate is something I haven't tried and not sure if it works as high sideband filter(branch filter) may get in trouble.  updating means summing up fraction of error to nco word. Remember how rtl nco works: its accumulator starts from some value (zero or so) then increments by nco tuning word. The update will add to accumulator value pushing its phase (nudging it) then you revert back to nco tuning word, this means you keep track of frequency lock after each update having now a different accum value but same inc step. eventually we hope the error becomes zero and keeps accum happy.  If you lose frequency lock then phase goes wild and out of control. I think software dsp don't have this problem as they control freq separately.

[Altera_Forum]

We hope error becomes zero... what about loop filter's output? do we hope it goes back to zero too? 

If we do, this is impossible because inside the loop filter there is register that is summed with the incoming error, once error is zero it accumulates with its self. If we don't, then NCO word will be updated with the steady value of LF's output (if we suppose it settles... but it actually does not).  This figure should depict what I mean:  http://www.alteraforum.com/forum/attachment.php?attachmentid=10476&stc=1   I believe we should hope for the loop filter's output to be zero or so too. Because as you mentioned we are not like DSP guys and our NCO is different.

[Altera_Forum]

I didn't set loop at sYmbol rate... I actually corrected it to be at sAmple rate... previously it was running at 50 MHz and accumulating every single sample several times.

[Altera_Forum]

as phase gets corrected then the error tends towards zero and this depends on how you have scaled it so the nudge effect gets less and less. Every nudge should reduce the amplitude of oscillation.  

 If it overshoots then the sense of error reverses and so on you have a loop error that pushes one way and pulls the other way from time to time (this is acceptable withing limits of specified jitter). In real world the loop keeps track of tx/rx oscillators freq and phase drift.   If phase is correct to begin with then error starts as zero and should stay stay zero. I notice that with single tone bu haven't tried more.

[Altera_Forum]

Error was not zero with single tone! it was sinusoidal! 

Anyway... I am really out of time and I got upset with this loop. Here is my plan and I want your advice... eventually I am not gonna use the NCO... because my RF Front-End is RTL-SDR Dongle which has its own tuner (VCO). Eventually I am going to calculate error in FPGA then convert it to something that RTL-SDR understands and feed it to shift the phase according to error... I think it has some register to control phase. As my main problem is controlling the NCO because it doesn't have a separate phase input, this is going to make my life much much easier. So, how about stop working on that and begin implementing using the RTL Dongle?