clear, clc, close all

PLOT = 0;

% Waterfall curves for coherent detection of PSK signals

M = [2 4 8 16 32];                              % Number of symbols in constellation
N = 50000;                                     % Number of symbols used for plots

EbNo = [0:2:20];                                  % Energy to noise ratio per bit

Results = zeros(length(M), length(EbNo));       % Matrix that stores results

for i = 1:length(M)
    i
    % Generate the symbols
    
    Symbol = zeros(M(i),2);
    for m = 0:M(i)-1
        Symbol(m+1,:) = [cos(2*pi*m/M(i)) sin(2*pi*m/M(i))];
    end
    
    % Plot the symbols
    
    if PLOT == 1
        figure, plot(Symbol(:,1), Symbol(:,2), 'o'), grid
        axis([-2 2 -2 2]);
        xlabel('real axis')
        ylabel('imaginary axis');
        title (['Constellaion for ' num2str(M(i)) '-PSK signal']);
    end
    
    for j = 1:length(EbNo)
        N0 = 1/(10^(EbNo(j)/10) * log2(M(i)));          % PSD of the thermal noise
        
        error = 0;                              % error counter
        
        for k = 1:N
            I = randperm(M(i));                 % pick a random constellation symbol
            sm = Symbol(I(1),:);
            sigma = sqrt(N0/2);                 % standard deviation of the noise
            n =  sigma * randn(1,2);            % noise sample
            r = sm + n;                         % received vector
            
            % compute distances 
            
            d = zeros(M(i),1);                  % distances vector
            for n = 1:M(i)                      % for every symbol 
                d(n) = norm(r-Symbol(n,:));     % determine Eucleadian distance
            end
            [Mn,J] = min(d);                    % find the symbol with minimum distance
            
            if J ~= I(1)                        % if the RX symbol is different than the TX symbol
                error = error + 1;              % increment the error counter
            end
        end
        Results(i,j) = error / N;               % P(e) = number of erronious symbols over
    end                                         % total number of symbols
end

% Generate the plots

figure, semilogy(EbNo, Results');
grid, axis([0 20 1e-3 1]);
xlabel('EbNo [dB]');
ylabel('probability of error');