%
% 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=1024;                % number of symbols (bits)
Tsymb=Ndata/Nsymb;        % symbol duration
f0=1024;                  % "carrier" frequency
noise_ampl=0.1;            % 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=cos(2*pi*f0/Ndata*t);

Eb=s_carrier*s_carrier'/Nsymb     % Energy per bit
%check=Tsymb/2
N0=noise1*noise1'/Ndata           % N0, noise power (in 1 Hz bandwidth)
EB_over_N0=Eb/N0

% 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=cos(2*pi*f0*tt/Ndata);

for ii = 1:Nsymb
    for jj = 1:Tsymb
        buf(jj)=s((ii-1)*Tsymb+jj);
    end

    dot_prod(ii)=buf*ss';
    dot_prod(ii)=dot_prod(ii)/Tsymb;

    if dot_prod(ii) < 0
        b(ii)=1;
    else
        b(ii)=0;
    end
end

tb=0:1:Nsymb-1;

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

% 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 \n', bit_seq(ii), b(ii), dot_prod(ii) )
     if bit_seq(ii) == b(ii)
         nmatch = nmatch+1;
     end
end 
nmatch
err=Nsymb-nmatch
