Simulating BPSK Costas Loop in Matlab

Question

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.

Answer 1

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

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 = 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.