% The Demodulation Function:
%Costas Loop Implementation for an BPSK signal
function y=demod(input)
len = length(input); %number of input samples
load b.txt;
load alpha_BPSK.txt;
load beta_BPSK.txt;
num=b;
hilbert_output = filter(num, 1, input);
%groupdelay should be accounted for generating a analytic signal
[gd,f] = grpdelay(num,1,512,'whole',8000);
groupdelay=gd(10);
analytic = input(1:len-groupdelay) + i*hilbert_output(groupdelay+1:len);
%Initialization
out = [];
phi = [];
phi(1) = 0;
temp_out1=0;
temp_pre_out1=0;
temp_out2=0;
temp_out3=0;
alpha=alpha_BPSK;
beta=beta_BPSK;
fc=2000;
%fs=8;
%Simulation
for I=1:len-groupdelay
phi(I)= temp_out3;
phase(I) = exp(-j*phi(I));

c1(I) = real(analytic(I)*phase(I));
c2(I) = imag(analytic(I)*phase(I));
out(I)= c1(I);
q(I) = c1(I)*c2(I);
temp_out1=temp_pre_out1+q(I)*beta;
temp_out2=alpha*q(I)+ temp_out1;
temp_out3= 2*pi/8+phi(I)+temp_out2;
temp_pre_out1=temp_out1;
end;
y=out;