# 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:

load RRC.mat

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:

% Siraj Muhammad
% 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