%  Simulator for AWGN channels
%  M-ary PSK
%
%  Michael Ng: sxn@ecs.soton.ac.uk

    clear;
    SNR_start=0; SNR_delta=1; SNR_stop=35; BER_stop=1e-4; 
    M=4;                  % modulation level: M=2,4,8,16,32,64      
    min_bit_errors = 200; % minimum number of bit errors
    
    % bits per modulated symbol
    bps = log2(M); 
    
    % M-PSK initialisation
    index=[0:M-1];
    psk_symbols = exp(-1i*2*pi*index/M);
    
    % Choose a file to save the results into
    filename = ['results_awgn_',num2str(M),'psk.mat'];
    
    % Doing the simulation frame-by-frame rather than bit-by-bit speeds it up. This is because we can use vector operations instead of longer loops.
    frame_length = 12000; % number of bits per frame
    symbol_frame_length = frame_length/bps; % number of symbols per frame

    % Setup the SNR for the first iteration of the loop.
    SNR_count = 1;
    SNR = SNR_start;
    BER = 1;

    % display format
    format short e;
    
    % Loop until the job is killed or until the SNR or BER target is reached.
    while SNR <= SNR_stop && BER >= BER_stop

        % Convert from SNR (in dB) to noise power spectral density.
        N0 = 1/(10^(SNR/10));
        
        % Counters to store the number of errors and bits simulated so far.
        error_count=0;
        bit_count=0;
       
        % Keep going until enough errors have been observed. This runs the simulation only as long as is required to keep the BER vs SNR curve smooth.
        while error_count < min_bit_errors
            
            % Generate some random bits.
            bits = round(rand(1,frame_length));

            % bit to integer conversion
            sym = bi2de(reshape(bits,bps,symbol_frame_length)','left-msb');
                        
            % M-PSK modulation
            symbols = psk_symbols(sym+1);
                       
            % AWGN noise
            noise = sqrt(N0/2)*(randn(1,symbol_frame_length)+1i*randn(1,symbol_frame_length));
            
            % Add some complex Gaussian distributed noise.
            rx_symbols = symbols + noise;
            
            % M-PSK demodulation
            %hard_symbol = int8(rx_symbols<0);
            for j = 1:symbol_frame_length
               index(j) = 1;
               norm_1 = abs(rx_symbols(j) - psk_symbols(index(j)));
               for m = 2:M
                  norm_2 = abs(rx_symbols(j) - psk_symbols(m));
                  if (norm_2 < norm_1)
                     index(j) = m;
                     norm_1 = norm_2;
                  end
               end
            end            
            hard_symbol = index - 1;
            
            % convert integer to bit sequence
            recovered_bits_mat = de2bi(hard_symbol,'left-msb')';
            recovered_bits = reshape(recovered_bits_mat,1,frame_length);    
            
            % Accumulate the number of errors and bits that have been simulated so far.
            error_count = error_count + sum(recovered_bits ~= bits);
            bit_count = bit_count + frame_length;
        end
        
        % Calculate the BER.
        BER = error_count/bit_count;
        
        % SNR per bit
        EbN0 = SNR - 10*log10(bps);
        
        % Store the SNR and BER in a matrix and display it.
        results(SNR_count,1) = SNR;
        results(SNR_count,2) = BER
        results(SNR_count,3) = EbN0;
        
        % Save the results into binary files. This avoids the loss of precision that is associated with ASCII files.
        save(filename, 'results', '-MAT');

        % Plot the results
        semilogy(results(:,1), results(:,2),'o-','MarkerSize',10);
        title([num2str(M) '-PSK modulation in an AWGN channel']);
        ylabel('BER');
        xlabel('SNR [dB]');
        ax = gca; ax.FontSize = 12; grid on;
        legend([num2str(M) '-PSK']);

        if BER_stop ~= 0 
            ylim([0.99*BER_stop,1]);
        end
        if SNR_stop ~= inf
            xlim([SNR_start, SNR_stop]);
        end
        
        % Setup the SNR for the next iteration of the loop.
        SNR = SNR + SNR_delta;
        SNR_count = SNR_count + 1;
    end
    


% %plotting multipole curves in one gragh
% load results_awgn_2psk.mat; BER_2 = results(:,2); SNR_2 = results(:,1);
% load results_awgn_4psk.mat; BER_4 = results(:,2); SNR_4 = results(:,1);
% load results_awgn_8psk.mat; BER_8 = results(:,2); SNR_8 = results(:,1);
% load results_awgn_16psk.mat; BER_16 = results(:,2); SNR_16 = results(:,1);
% load results_awgn_32psk.mat; BER_32 = results(:,2); SNR_32 = results(:,1);
% load results_awgn_64psk.mat; BER_64 = results(:,2); SNR_64 = results(:,1);
% 
% figure(2);
% semilogy(SNR_2, BER_2,'o-','MarkerSize',10);
% hold on
% semilogy(SNR_4, BER_4,'x-','MarkerSize',10);
% semilogy(SNR_8, BER_8,'+-','MarkerSize',10);
% semilogy(SNR_16, BER_16,'*-','MarkerSize',10);
% semilogy(SNR_32, BER_32,'s-','MarkerSize',10);
% semilogy(SNR_64, BER_64,'d-','MarkerSize',10);
% title('M-PSK modulation in an AWGN channel');
% ylabel('BER');
% xlabel('SNR [dB]');
% ax = gca; ax.FontSize = 12; grid on;
% ylim([9.99e-6,1]);
% xlim([SNR_64(1), SNR_64(end)]);
% legend('2-PSK','4-PSK','8-PSK','16-PSK','32-PSK','64-PSK');
% hold off
