%% 802.11b 1Mbps PHY link.
% This M code simulates DBPSK modulation and barker code spreading
% on a perfectly synchonized 802.11b link. It calculates the BER
% rate at each EsNo (EbNo) and plots the result.

%% Simulation parameters
% For each BER loop, we specify the number of packets to transmit, the packet size and
% the range of channel EsNo values to test.

EsNoRange=[0:2:10]; % Range of noice levels to calculate BER
NumPackets=2;
PacketSizeBytes=1024;
PacketSizeBits=PacketSizeBytes*8; % Here we ignore preamble and sync bits
clear BERResults;

%% System parameters and constants
% Specify a number of system constants.

% Spreading parameters
Barker=[1 -1  1  1 -1  1  1  1 -1 -1 -1]'; % Barker sequence
SpreadingRate=length(Barker);              % Spreading rate

% Upsampling rate
SamplesPerChip=8;

% Filter order and coefficients - root raised cosine
FilterOrder=40; % Set to multiple of SamplesPerChip to make delay calculation easy
h=firrcos(FilterOrder,7e6,.7,88e6,'rolloff','sqrt',FilterOrder/2,kaiser(FilterOrder+1,1));

%% Delay calculation
% Calclate (specify) the net number of bits delay in the link due
% to the filtering.
%
% * samples_delay = 2 filters x (40 coeffs / 2) = 40 samples
% * chips_delay = sample_delay/SamplesPerChip = 40/8 = 5 chips
% *  Must recalculate delay if you change any of these parameters
% We must delay the signal 6 more chips to align it with the 11 chip
% boundary. This results in an 11 chip delay or one symbol/bit delay.
% You must recalculate total and additional delay if you change any of these
% parameters

BitDelay=1;
ChipDelayAdd=6;

%% Main BER loop
% Calculates the BER for each EsNo level. 

NumEsNos=length(EsNoRange);

disp(' ');disp('Start Simulation');
for EsNoIndex =1:NumEsNos
    
    EsNo=EsNoRange(EsNoIndex);
    disp(['Simulating: EsNo=' num2str(EsNo) 'dB']);
    SNR=EsNo+10*log10(1/SpreadingRate)+10*log10(1/SamplesPerChip);

    % Initialize system and simulation measurements state

    % Bits
    TotalBits=logical(0); % Bit count for BER calculation
    ErrorBits=logical(0); % Error count for BER calculation
    LastTxSymbol=1;       % Set DBPKS Modulator state
    LastRxSymbol=1;       % Set DBPKS Demodulator state

    % Filters
    Rx_chips_delayed_store=zeros(ChipDelayAdd,1);
    Tx_bits_delayed_store=logical(zeros(BitDelay,1));

    Tx_Filter_State=h(1:end-1); % Fill filter with a +1 symbol
    Rx_Filter_State=h(1:end-1); % Fill filter with a +1 symbol

    % Main simulation loop
    % Each packet is transmitted, and the recieved bits compared with the
    % transmitted bits to calculate the BER.
    for Packet=1:NumPackets

        % Construct frame of bits
        Tx_bits=rand(PacketSizeBits,1)>.5;           % Random bits

        % Modulate
        Tx_bits_bp=(1-2*Tx_bits);                    % Convert to bipolar 0,1 --> 1, -1
        Tx_symbols=LastTxSymbol*cumprod(Tx_bits_bp); % New DBPSK symbol = previous  * 1 or -1
        LastTxSymbol=Tx_symbols(end);                % Store modulator state (last symbol)

        % Spread symbols with Barker code, upsampling by spreading rate
        Tx_chips=reshape(Barker*Tx_symbols',[],1); % Multiply by barker and reshape to a columm
        Tx_chips=complex(Tx_chips); % Make complex to ensure correct baseband transmission

        % Upsample chips by SamplesPerChip factor
        Tx_samples=zeros(length(Tx_chips)*SamplesPerChip,1); % Create empty Tx_samples
        Tx_samples(1:SamplesPerChip:end,1)=sqrt(SamplesPerChip)*Tx_chips; % Normalize power due to upsampling

        % Tx Filter
        [Tx_samples_filtered,Tx_Filter_State]=filter(h,1,Tx_samples,Tx_Filter_State); % Filter
        Tx_samples_filtered=Tx_samples_filtered*2.495; % Set output power to 1W
        var(Tx_samples_filtered); % Calculate Tx signal power, view by removing ';'

        % Transmit though AWGN Channel assuming 0dBW input power (check
        % with line above)
        Rx_samples_unfiltered = awgn(Tx_samples_filtered,SNR,0); 

        % Rx Filter
        [Rx_samples_filtered,Rx_Filter_State]=filter(h,1,Rx_samples_unfiltered,Rx_Filter_State);

        % Downsample - sample chips
        Rx_chips=Rx_samples_filtered(1:SamplesPerChip:end);

        % Add 1 chip delay to move signal to 11 chip boundary
        Rx_chips_delayed=[Rx_chips_delayed_store; Rx_chips(1:end-ChipDelayAdd)];
        Rx_chips_delayed_store=Rx_chips((end-ChipDelayAdd+1):end); % Store delayed chips

        % Despread - sample symbol
        Rx_symbols=Barker'*reshape(Rx_chips_delayed,SpreadingRate,PacketSizeBits); % Multiply by Barker
        Rx_symbols=Rx_symbols(:)/SpreadingRate;   % Make a column and normalize

        % Demodulate
        Rx_symbols_plus_last=[LastRxSymbol; Rx_symbols];
        Rx_symbols_plus_last_mult=Rx_symbols_plus_last(1:end-1).*conj(Rx_symbols_plus_last(2:end));
        Rx_bits=Rx_symbols_plus_last_mult < 0; 
        LastRxSymbol=Rx_symbols(end);         % Demodulator state

        % Calculate BER

        % Add BitDelay to Tx signal to align with Rx signal
        Tx_bits_delayed=[Tx_bits_delayed_store; Tx_bits(1:end-BitDelay)];
        Tx_bits_delayed_store=Tx_bits(end-BitDelay+1:end); % Store delayed symbol

        if Packet==1 % Ignore delayed bits on first packet
            TotalBits=TotalBits+length(Rx_bits)-BitDelay;
            ErrorBits=ErrorBits+sum(Tx_bits_delayed(BitDelay+1:end)~=Rx_bits(BitDelay+1:end));
        else
            TotalBits=TotalBits+length(Rx_bits);               % Calculate total bits
            ErrorBits=ErrorBits+sum(Tx_bits_delayed~=Rx_bits); % Compare Tx and Rx bits
        end
    end

BERResults(EsNoIndex)=ErrorBits/TotalBits; % Calculate BER

end

%% Plot BER Results
% Plot the BER results Vs EsNo.

semilogy(EsNoRange,BERResults,'*-');
grid;
title('802.11b 1Mbps DBPSK BER');
ylabel('BER')
xlabel('EsNo');