Simulating BPSK Costas Loop in Matlab

I am trying to implement BPSK demodulator in FPGA. First, I am modeling the system in MATLAB. Now, I am trying to design Costas Loop for carrier synchronization. This is the block diagram I'm trying to achieve: I'm using RRC as a low pass filter. This is the MATLAB code:

fc = 0.05; % 0.11  %0.09
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);
Theta = zeros(1, N);

% Loop Filter Coefficients
BW = 100;  % 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);

for i = 2:N
% Downconverting to Baseband
bb(i) = rf_signal(i).*exp(j*2*pi*fc*i).*exp(j*Theta(i-1));

% Filtering
bb_f = filter(RRC, bb(1:i));
I_f(i) = real(bb_f(i));
Q_f(i) = imag(bb_f(i));

% 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);

% Phase Accumulator
Theta(i) = Theta(i-1) + PhErr(i);

end

figure; subplot 221; plot(real(bb)); title('I Channel')
subplot 222; plot(imag(bb)); title('Q Channel')
subplot 223; plot(I_f); title('I Channel Filter')
subplot 224; plot(Q_f); title('Q Channel Filter')

figure; subplot 311; plot(error); title('Phase Error');
subplot 312; plot(PhErr); title('Loop Filter');
subplot 313; plot(Theta); title('Control Signal \theta');

rf_signal comes from another script that generates BPSK signal. However, when I close the loop and set frequency exactly as the transmitter's I get these waveforms:  Although the signal demodulates very well when I open the loop. I have tested the loop with a cosine wave and it works just fine by producing a DC signal.

Any advice would be appreciated.

• I'm wondering, why are you using a quadrature demodulator for receiving a BPSK signal? – MBaz Mar 25 '15 at 15:09
• This is how Costas Loop is designed. Plus, you can use IQ demodulator for any kind of signals (to the best of my knowledge). Here in BPSK, I will ignore the Q channel after the carrier is synchronized since it carries no data. – Siraj Muhammad Mar 26 '15 at 7:03
• Fair enough -- I didn't realize at first sight that you're not using the quadrature channel for data. – MBaz Mar 26 '15 at 15:27
• Thanks for your kindness. That's what I looking for the code of Costas loop example. But, I have some problems of running the code. Could I get the RRC.mat and BPSK signal? I'm looking forward to your reply. – user16108 Jun 5 '15 at 6:17

I figured it out. It turned out that I don't really need that 'Phase Accumulator', and I should reverse the control signal in the NCO. So, I use -PhErr(i-1) as a control signal to NCO. Here is the modified code:

% 25/3/2015
% BPSK Demodulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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 = 100;  % 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);

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(i) = real(bb_f(i));
Q_f(i) = imag(bb_f(i));

% 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');

I added a variable phase_offset to test the loop with different offsets, and I can drift the frequency a little bit and the loop will lock very nicely.

I hope this helps somebody as I couldn't find any MATLAB code for Costas Loop online.

• There is an analysis and Matlab code in the book "Telecommunications Breakdown" by Johnson and Sethares. A pre-print of the book is available for free at sethares.engr.wisc.edu/telebreak.html. Sometimes it's better to search in books than online :) – MBaz Mar 26 '15 at 15:26
• Brilliant! It's good to be recommended a book. I was looking for a similar one. Thanks. – Siraj Muhammad Mar 26 '15 at 20:15

protected by jojek♦Jun 5 '15 at 8:13

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?