%
% This simple code shows how a coherent demodulator/decoder works.
% The deciding waveform and the incoming telemetry stream are supposed
% to be synchronized (this can be done in several ways). We suppose also 
% that the Doppler shift has been removed by means of a phase-locked loop
% (PLL). 

clear all;
close all;

% Define the sequence length and parameters first.
% One can play with the number of symbols transmitted in the interval
% (Nsymb), which defines the symbol duration expressed in terms of samples
% in the sequence (symb), and the noise amplitude  (noise_ampl), relative to
% the "carrier", assumed to have unit amplitude. 
 
Ndata=2^15          % sequence length (number of samples)
t=0:1:Ndata-1;      % sequence length
Nsymb=128;          % number of symbols (bits)
Tsymb=Ndata/Nsymb;  % symbol duration
f0=1024;            % "carrier" frequency
noise_ampl=2.;     % noise amplitude (relative to signal; signal = 1)


% sequence of bits (simply, a random sequence) 
bit_seq=rand(1,Nsymb);
for ii = 1:Nsymb
    if bit_seq(ii)-0.5 > 0 
        bit_seq(ii)=0;
    else
        bit_seq(ii)=1;
    end
end

%bit_seq=[ 0 1 1 1 0 0 0 1 0 1 0 1 1 1 0 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 0 1 ];
%z= polyval([ 0 1 1 1 0 0 0 1 0 1 0 1 1 1 0 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 0 1 ],2)
%z=polyval(bit_seq,2)
% this function evaluates the polynomial whose coefficients are the vector elements,
% at x=2; assume that the number represented is an unsigned integer

% This inelegant loop applies the appropriate phase shift to each group of samples.
% Quick and dirty.
for ii = 0:Ndata-1
  jj1=ii/Tsymb;
  jj=fix(jj1);
  phi_m(ii+1)=bit_seq(jj+1); 
end

figure ('Name','Bit sequence'); 
plot (phi_m,'LineWidth',2)
axis([0 inf -0.5 1.5])

% Generating the additive noise (AWGN) and the carrier.
noise1=noise_ampl*randn(1,Ndata);
s_carrier=sin(2*pi*f0/Ndata*t);

% MODULATOR. The carrier is mixed with the modulating signal in a binary
% phase shift keying (BPSK) scheme.
for ii = 0:Ndata-1
  s_mod(ii+1)=s_carrier(ii+1)*cos(pi*phi_m(ii+1));
end

% adding noise
s=s_mod+noise1;


figure('Name','Signal and signal+noise');
%plot(t,s,t,s_mod,'LineWidth',4,)
plot(t,s)
hold all
plot(t,s_mod,'LineWidth',4)

figure ('Name','Modulated carrier + noise, and bit transitions');
plot (t,s,t,phi_m,'LineWidth',2)
% Magnify a stretch of the plot to see what happens at bit transitions.


% FFT of the modulated signal and the carrier. Do you note something interesting?
sf=fft(s);
sc=fft(s_carrier);
Psf=sf.*conj(sf)/Ndata;
Psc=sc.*conj(sc)/Ndata;


figure ('Name','FFT of the modulated signal and the carrier');
plot(t,Psf,t,Psc)


% demodulator/decoder
tt=0:1:Tsymb-1;
ss=exp(i*pi/2-2*pi*i*f0*tt/Ndata);          % for a sin carrier
%ss=exp(2*pi*i*f0*tt/Ndata);                % for a cos carrier
%figure (121);
%plot (tt,imag(ss),tt,real(ss))

for ii = 1:Nsymb
    for jj = 1:Tsymb
        buf(jj)=s((ii-1)*Tsymb+jj);
    end
%    figure (111);
%    plot (tt,buf,tt,imag(ss),tt,real(ss))
%    pause
    dot_prod(ii)=buf*ss';
    dot_prod(ii)=dot_prod(ii)/Tsymb;
    phase=atan2(imag(dot_prod(ii)),real(dot_prod(ii)))*180/pi;
    if abs(phase) < 90
        b(ii)=0;
    else
        b(ii)=1;
    end
end


tb=0:1:Nsymb-1;
%figure (12);
%plot (tb,b,'d',tb,bit_seq,'+','MarkerEdgeColor','r')


figure ('Name','Matched filter output (phase)');
plot (dot_prod,'+','MarkerEdgeColor','r')
grid on
axis square
axis equal

% I compute now the number of errors committed by the decoder.
nmatch=0;
for ii = 1:Nsymb
%    disp (bit_seq(ii))
     fprintf('%i %i %f %f\n', bit_seq(ii), b(ii), real(dot_prod(ii)), imag(dot_prod(ii)) )
     if bit_seq(ii) == b(ii)
         nmatch = nmatch+1;
     end
end 
nmatch
err=Nsymb-nmatch