% This script illustrates a number of basic ideas in wireless
% communications:
% 1 -- It translates the math found Rappaports textbook
%      (chapter 6) into simulation code similar to that found in software
%      programmable modems.
% 2 -- It illustrates what signals are present at various points in the
%      modulation of a BPSK transmitter and receiver.
% 3 -- It demonstrates the error rate of a modem in additive Gaussian
%      noise.
% 4 -- It shows that multiple users can interfere with each other without
%      some mechanism to isolate them.
%
% Compare this to what you get by calling "commdoc_mod" at the command
% prompt.

clear all
close all

Rb = 30e3;     % 30 kbits/s

% default settings ----------------
twoUsers     = false;
runBER       = false;
showSpectra  = true;
useAWGN      = false;
useFDMA      = false;
P            = 60;
Nbits        = 100;
EbNo_dB      = 7;
freq_offset  = 10 * Rb;
% end of default settings ---------

demoNum = 1;
switch demoNum,
    case 1,     % single user 60x oversampled
        % no change from default choices above
    case 2,     % single user plus noise (12 dB) to look at plots
        useAWGN      = true;
        EbNo_dB      = 14;
    case 3,     % single user plus noise (7 dB) to look at plots
        useAWGN      = true;
        EbNo_dB      = 7;
    case 4,     % single user plus noise (7 dB) to check BER
        useAWGN      = true;
        EbNo_dB      = 7;
        Nbits        = 10e3;
        showSpectra  = false;
    case 5,     % single user only 4x oversampled
        P            = 4;
        Nbits        = 5e3;
    case 6,     % two users, same frequency
        twoUsers     = true;
        useFDMA      = false;
    case 7,     % two users, distinct frequencies
        twoUsers     = true;
        useFDMA      = true;
        freq_offset  = 10 * Rb;
    case 8,     % tinker here
        twoUsers     = false;
        useFDMA      = true;
        freq_offset  = 10 * Rb;
        useAWGN      = true;
        EbNo_dB      = 10;
        Nbits        = 1e3;
end

useRRC       = false;

Tb = 1/Rb;
fc = 300e3;

EbNo_lin = 10.^(EbNo_dB/10);

modType = 'bpsk';        M = 1;
fecType = 'none';        R_fec = 1;

% generate random data
bits = round( rand( Nbits, 1 ) );
                                if showSpectra,
                                    figure(1); stem( bits(1:50) );
                                    title( 'input bit stream' );
                                end

% map data to antipodal values
mapped = 1 - 2*bits;
                                if showSpectra,
                                    figure(2); stem( mapped(1:50) );
                                    title( 'mapped bit stream' );
                                end

% upsample by N
upsampled = upsample( mapped, P );
                                if showSpectra,
                                    figure(3); stem( upsampled( 1:50*P ) );
                                    title( 'upsampled bit stream' );
                                end

% use a square pulse shape
pulse = ones( P, 1 );
shaped = filter( pulse, 1, upsampled );
                                if showSpectra,
                                    figure(4); plot( shaped( 1:50*P ) );
                                    title( 'pulse-shaped bit stream' );
                                    
                                    [ Pxx, F ] = pwelch(shaped,256,[],[],Rb*P,'twosided'); 
                                    figure(5); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
                                    title( 'power spectrum of pulse-shaped bit stream' );
                                    
                                    % figure(204); freqz(LPF);
                                end

% modulate
t = Tb/P * [0:Nbits*P-1]';
carrier = cos( 2*pi * fc * t ); 
modulated = carrier .* shaped;
                                if showSpectra,
                                    figure(6); plot( modulated( 1:50 * P ), '-x' );
                                    title( 'modulated bit stream' );
                                    
                                    [ Pxx, F ] = pwelch(modulated,256,[],[],Rb*P,'twosided'); 
                                    figure(7); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
                                    title( 'power spectrum of modulated bit stream' );                                    
                                end



if twoUsers,
    % generate a pulse-shaped, unique bit stream for second user
    bits2 = round( rand( Nbits, 1 ) );
    mapped2 = 1 - 2*bits2;
    upsampled2 = upsample( mapped2, P );
    shaped2 = filter( pulse, 1, upsampled2 );
    
    if useFDMA,
        % modulate with offset carrier
        offset_carrier = cos( 2*pi * (fc+freq_offset) * t );
        modulated2 = offset_carrier .* shaped2;
    else
        % modulate with same carrier
        modulated2 = carrier .* shaped2;
    end
    
    if showSpectra,
        [ Pxx, F ] = pwelch(modulated2,256,[],[],Rb*P,'twosided');
        figure(8); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
        title( 'power spectrum of second modulated stream' );
    end

    tx_out = modulated + 2 * modulated2;
    
                                if showSpectra,
                                    [ Pxx, F ] = pwelch(tx_out,256,[],[],Rb*P,'twosided'); 
                                    figure(9); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
                                    title( 'power spectrum of both bit streams' );
                                end
else
    tx_out = modulated;
end


if useAWGN,
%     Es = mean( abs( tx_out ) .^ 2 );
%     sigma = sqrt( Es * P / (2 * EbNo_lin * R_fec * M) );
%     noise = sigma * ( randn( size( tx_out ) ) + j * randn( size( tx_out ) ) );
%     rx_in = tx_out + noise; 
    
    % add in a tiny complex component to get complex additive noise...
    rx_in = awgn( tx_out + j*eps, EbNo_dB + 10*log10( M/P ), 'measured' );
    
else
    rx_in = tx_out;
end

% demodulate signal 1
demod = rx_in .* carrier;
                                if showSpectra,
                                    figure(10); plot( real(demod) );
                                    title( 'demodulated signal' );
                                    
                                    [ Pxx, F ] = pwelch(demod,256,[],[],Rb*P,'twosided'); 
                                    figure(11); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
                                    title( 'power spectrum of demodulated signal' );
                                end



% low pass filter (LPF)
LPF_delay = 3;
LPF = firpm( P*2*LPF_delay, [ 0 0.5*4/P 0.5*4/P+0.1 1 ], [ 1 1 0 0 ] );
LPF_filtered = filter( LPF, 1, demod );
                                if showSpectra,
                                    figure(12); plot( real(LPF_filtered(1:50*P)) );
                                    title( 'low-pass filtered, demodulated signal' );
                                    
                                    [ Pxx, F ] = pwelch(LPF_filtered,256,[],[],Rb*P,'twosided'); 
                                    figure(13); plot( F - max(F)/2, fftshift(10*log10( Pxx )) )
                                    title( 'power spectrum of low-pass filtered signal' );
                                    
                                end

% MF = matched filtered
MF = filter( fliplr( pulse ), 1, LPF_filtered );
                                if showSpectra,
                                    figure(14); plot( real(MF) );
                                    title( 'match-filtered, signal in receiver' );
                                end
                                
MF_trunc = MF( LPF_delay * P + 1:end );

% dumped
I = real( MF_trunc );
Q = imag( MF_trunc );
dumped_I = I( P:P:end );
dumped_Q = Q( P:P:end );
                                if showSpectra,
                                    figure(15); plot( reshape( I(P+1:end), P, length(I)/P-1 ) );
                                    title( 'eye diagram' );
                                    figure(16); plot( dumped_I(1:50), dumped_Q(1:50), '.' );
                                    title( 'scatter plot' ); axis equal;
                                    figure(17); stem( dumped_I( 1:50 ) );
                                    title( 'sampled MF outputs (integrate and dump)' );
                                end

% detect
detected = sign( dumped_I );
f = find( detected == 0 );
detected(f) = 1;  % arbitrarily assign indeterminate detections (equidistant 
                  %                     from both +1 and -1) the value one
%                                 if showSpectra,
%                                     % this plot correlates input, output
%                                     % data streams
%                                     figure(18); plot(xcorr(mapped, detected))
%                                 end

% demap
outBits = 0.5 - detected/2;

errs = sum( bits(1:end-LPF_delay) ~= outBits )
BER = errs/(Nbits-LPF_delay)


