% Code for the Digital Communication 1 laboration
% 050210, 060227, 080207, 090216 Daniel Aronsson

clear
close all

channel='simulated';


%%%%%%%%%%%%%%%%%%%%%%%%%
% System variables
%%%%%%%%%%%%%%%%%%%%%%%%%

Nbits=3*4*4*11*19;          % Divisible by 3(for 8PSK),4(for BPSK,QPSK,16QAM) 
                            % another 4(for Hamming7_4),11(for Hamming15_11), Nbits=about 10000
coderate=;
Ncodedbits=round(Nbits/coderate);  % (round() removes machine rounding errors)
M=;                         % log2(M)=bits/symbol (e.g M=4 => QPSK)
Nfilter=8;                  % The length of the send/recv filter in symbols (even integer)
Nsuffix=100;                % Silence in the beginning to allow for the sound card to start up
Nprefix=100;                % Silence in the end to account for the delay in the system
Ntot=Ncodedbits/log2(M)+Nsuffix+Nprefix+32;  % Total number of symbols (32=length of training sequence)
fs=44100;                   % Sampling frequency (Hz)
R=fs/80;                    % Symbol rate (Hz). Make sure fs/R=integer.
T=Ntot/R;                   % Total time for the transmission (sec)
N=round(fs*T);              % Total number of samples (=fs/R * Ncodedbits/log2(M) i.e an integer)
                            % (round() removes machine rounding errors)
t=linspace(0,T-1/fs,N);     % Time axis
f=fftshift(linspace(-fs/2,fs/2-1/T,N)); % Frequency axis. Good for plotting the spectrum.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the baseband signal (complex).
% Create the baseband signal for the training sequence as it will
% look on the receiver side (after passing both send/recv filters).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=round(rand(Nbits,1));   % Meddelande

%%% Block 1 : Channel code
%c=some_encoder(m);  % INSERT FUNCTION CALL HERE

%%% Block 2 : Mapping onto symbols
%[I,Q]=some_modulator(c);   % INSERT FUNCTION CALL HERE

load trainingsequence  % Load the training sequency (random QPSK symbols)
            % It is needed to find the beginning of the signal on
            % the receiver side.
s=[zeros(Nprefix,1); trseq; I+sqrt(-1)*Q; zeros(Nsuffix,1)]; % Symbols to be transmitted

%%% Block 3 : Pulse shaping (square root raised cosine)
tx_pulsetrain=kron(s,[1; zeros(fs/R-1,1)]);  % Upsample by factor fs/R
pulse = srrc(fs,R,Nfilter,0.5);
tx_s_bb=fftfilt(pulse,tx_pulsetrain);

% The training signal (after passing both send and recv filters)
rx_trainingsequence=fftfilt(pulse,fftfilt(pulse,kron([trseq; zeros(2*Nfilter-1,1)],[1; zeros(fs/R-1,1)])));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create passband signal
% Send and receive passband signal.
% Send and receive "silence" (for measuring the noise).
% Mix down to baseband and filter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(channel,'soundcard')
    
    % Play and record silence for noise measurement
    pause(1)
    wavplay(zeros(N,1),fs,'async');
    rx_n_pb  = wavrecord(N+100, fs,1);
    rx_n_pb=rx_n_pb(101:end);
    rx_n_pb=rx_n_pb-mean(rx_n_pb);
    
    %%% Block 4 : Carrier modulation
    wc=2000*2*pi;
    tx_s_pb=real(tx_s_bb).*cos(wc*t.')-imag(tx_s_bb).*sin(wc*t.');
    tx_s_pb=tx_s_pb/max(abs(tx_s_pb));
    
    %%% Block 5 : Channel (send and receive)
    pause(1)
    wavplay(tx_s_pb,fs,'async');
    rx_s_pb = wavrecord(N+100, fs,1);
    rx_s_pb = rx_s_pb(101:end);
    rx_s_pb = rx_s_pb-mean(rx_s_pb);
    
    %%% Block 6 : Mix down back to baseband
    rx_s_bb=2*rx_s_pb.*cos(wc*t.')-2*i*rx_s_pb.*sin(wc*t.');
    rx_n_bb=2*rx_n_pb.*cos(wc*t.')-2*i*rx_n_pb.*sin(wc*t.');
    
elseif strcmp(channel,'simulated')
    
    % Simulated noise
    Npow=0.03; % noise power. Experiment with this to produce different Eb/N0.
    n=randn(Ntot*fs/R,1)*sqrt(Npow);
    rx_n_pb=n;
    rx_n_pb=rx_n_pb-mean(rx_n_pb);
    
    %%% Block 4 : Carrier modulation    
    wc=2000*2*pi;
    tx_s_pb=real(tx_s_bb).*cos(wc*t.')-imag(tx_s_bb).*sin(wc*t.');
    tx_s_pb=tx_s_pb/max(abs(tx_s_pb));

    %%% Block 5 : Channel (add noise)    
    rx_s_pb = tx_s_pb + n;
    rx_s_pb = rx_s_pb-mean(rx_s_pb);

    %%% Block 6 : Mix down back to baseband    
    rx_s_bb=2*rx_s_pb.*cos(wc*t.')-2*i*rx_s_pb.*sin(wc*t.');
    rx_n_bb=2*rx_n_pb.*cos(wc*t.')-2*i*rx_n_pb.*sin(wc*t.');
    
else
    error('Unknown channel (check the channel variable!)');
end

%%% Block 7 : Filter (square root raised cosine)
rx_s_bb_filtered=fftfilt(pulse,rx_s_bb);
rx_n_bb_filtered=fftfilt(pulse,rx_n_bb);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use correlation to locate the beginning of the signal. Then sample.
% Measure Eb/N0 and BER.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isequal(round(xcorr_dc1([0 1 2 3],[1 2 3])) , [0 0 3 8 14 8 3])
    error(['Your xcorr.m is old. Update to Signal Processing Toolbox' ...
            ' 5 or more.']);
end

[corr,lags] = xcorr_dc1(rx_s_bb_filtered,rx_trainingsequence);
[cmax,nmax] = max(corr);
delay       = lags(nmax);
h           = cmax/sum(abs(rx_trainingsequence).^2);
k           = delay+(32+Nfilter)*fs/R+1;

try
    fprintf('Index : %d, scaling factor : %f+i*%f\n',k,real(h),imag(h));
    %%% Sampling (denoted by a switch in the block diagram)
    rx_s_bb_sampled=rx_s_bb_filtered(k:fs/R:k+fs/R*Ncodedbits/log2(M)-1)/h;
    rx_n_bb_sampled=rx_n_bb_filtered(k:fs/R:k+fs/R*Ncodedbits/log2(M)-1)/h; 
catch
    errmsg = lasterr;
    if(strfind(errmsg, 'Index exceeds matrix dimensions'))
        fprintf('Lost sync!\n');
    else
        error(errmsg);
    end
end

%%% Block 8 : Detect and channel decode. INSERT FUNCTION CALL HERE.
%estimated_c=some_detector(real(rx_s_bb_sampled),imag(rx_s_bb_sampled));
%estimated_m=some_decoder(estimated_c);

BER=sum(estimated_m ~= m) / length(m);
EbN0=10*log10(1/coderate/log2(M) * (var(rx_s_bb_sampled)/var(rx_n_bb_sampled)-1) );
fprintf('BER : %f\n',BER);
fprintf('EbN0 : %f\n',EbN0);