SNR_dB = 15;
M = 4;             % number of points in constellation

% Q1 what does the constellation look like when M = 4?
% Q2 what does the constellation look like when M = 32?

% -------------------------------------------------------------------------
% Supported modulations
% -------------------------------------------------------------------------

% qpsk
angles = [45 135 -45 -135 ] * pi/180;
LUT_qpsk = [ cos(angles) ; sin(angles) ]';

% Q3 provide code to support M = 8, where each point is simply 22.5 + k*45 degrees


% 32-apsk
LUT_32apsk = [ 92, -38, -92, 38, 100, 20, -20, 47, 64, 0, -71, 17, -17, 20, -20, 71, -100, -47, -92, -64, 17, 38, -17,  0, 71, 47, 92, 64, -47, -38, -64, -71 ; ...
               38, 92, 38, 92, 0, 20, 20, 47, 17, 100, 71, 64, 64, -20, -20, 71,  0, 47, -38, 17, -64, -92, -64, -100, -71, -47, -38, -17, -47, -92, -17, -71 ]';


decode_8psk_string = 0;

switch M,
    case 4
        LUT = LUT_qpsk;
    case 8
        LUT = LUT_8psk;
        decode_8psk_string = 1;
    case 32
        LUT = LUT_32apsk;
    otherwise
        error 'unsupported modulation type'
end

LUT_complex = LUT * [ 1 ; j ];

% -------------------------------------------------------------------------
% Generate some data of the appropriate modulation
% -------------------------------------------------------------------------

log2M = log2(M);
N = 1000; % numSymbols;
data = round( rand( N * log2M, 1 ) );

if decode_8psk_string,
    data2 = [
        0     1     0     0     0     0     1     0     0     1     1     0     0     0     0     1     0 , ...
        1     1     1     0     0     1     0     0     1     1     1     0     1     0     0     0     0 , ...
        1     0     0     0     0     0     0     1     0     1     0     0     1     1     0     1     1 , ...
        0     1     0     0     1     0     1     1     0     1     1     0     1     0     1     1     1 , ...
        0     0     0     0     0     1     1     1     0     0     1     1     0     1     1     0     1 , ...
        1     1     1     0     1     1     0     1     1     1     0     0     0     1     0     0     0 , ...
        0     0     0     1     1     0     1     1     0     0     0     1     1     0     1     0     0 , ...
        1     0     1     1     0     1     0     1     1     0     1     1     0     0     1     0     1 , ...
        0     1     1     1     0     0     1     1     0     0     1     0     0     0     0     0     0 , ...
        0     1     1     1     0     0     0     0     0     1     0     1     1     0     1     0     1 , ...
        1     1     0     0     0     0     0     1     1     1     0     0     1     1     0     1     1 , ...
        0     1     0     1     1     0     0     1     0     0     0     0     0     0     1     1     0 , ...
        1     1     0     1     0     1     1     0     1     1     1     1     0     1     1     0     0 , ...
        1     0     0     0     1     1     1     0     1     0     1     0     1     1     0     1     1 , ...
        0     0     0     1     1     0     0     0     0     1     0     1     1     1     0     1     0 , ...
        0     0     1     1     0     1     0     0     1     0     1     1     0     1     1     1     1 , ...
        0     1     1     0     1     1     1     0     0     0     1     0     1     1     1     0
        ];
    data( 1:length( data2(:) ), 1 ) = data2(:);
end

% Q4 what alpha-numeric message pops out when decode_8psk_string = 1?  (Relies on Q3).



data_reshape = reshape( data, log2M, N )';
powers_of_2 = 2 .^ [log2M - 1 : -1 : 0 ]';
indices = data_reshape * powers_of_2;
I = LUT(indices+1,1);
Q = LUT(indices+1,2);
IQ = I + j*Q;

figure( 1 ); plot( real(IQ), imag(IQ), '-x' ); 
title('data points and transitions'); xlabel('I'); ylabel('Q');

% % Uncomment from here ................
% % -----------------------------------------------------------------------
% % Look at the data when pulse shaped
% % 
% % We could upsample, filter, modulate, send it through the channel, and
% % then reverse these operations.
% % Instead, we'll just apply them to have a look and then drop back to
% % the unfiltered data for simplicity.
% % % ---------------------------------------------------------------------
% 
% % Upsample by 4 and apply root-raised cosine filtering...
% UP = 4;
% IQup = upsample( IQ, UP );
% % complete the following function call for a root raised cosine filter
% % which is oversampled by UP = 4, alpha = 0.3, filter spans 6 symbols
% h = rcosine( fill in the proper arguments here );  
% IQrcos = filter( h, 1, IQup );
% figure( 2 ); subplot( 2, 1, 1 ); plot( real(IQrcos(end-199:end)) );
% title('pulse-shaped data');
% legend('I data');
% subplot( 2, 1, 2 ); plot( imag(IQrcos(end-199:end)) );
% legend('Q data');
% xlabel('time (samples)');
% 

%%% Q5 How does figure 2 look for M = 4
%%% Q6 How does Figure 2 look for M = 32
%%% Q7 Which of these modulations (QPSK or 32-APSK) is most robost to distortion and noise?


% % here's where you'd modulate by a carrier
% % we'll keep this baseband, add noise, and detect the data
% % Uncomment down to here ^^^^^^^^^^^^^^^^^^^

% -------------------------------------------------------------------------
% Add noise
% -------------------------------------------------------------------------

S = mean( abs(LUT_complex) .^ 2 );
No = S / 10^(SNR_dB/10);
sigma = sqrt( No / 2 );

IQreceive = IQ + sigma * ( randn( size(IQ) ) + j*randn( size(IQ) ) );

figure(3); plot( real(IQreceive), imag(IQreceive), 'b.', ...
    real(IQ), imag(IQ), 'rx' ); axis equal
title('Scatter plot'); xlabel('I'); ylabel('Q'); 
legend('received data', 'transmit values');

% Q8 How does figure 2 look for M = 4
% Q9 How does Figure 2 look for M = 32
% Q10 Which of these modulations (QPSK or 32-APSK) is most robost to distortion and noise?


% -------------------------------------------------------------------------
% Detect which constellation point is closest to each recieved symbol
% -------------------------------------------------------------------------

index = zeros(N,1);
for ii = 1:N,
    distances = abs( LUT_complex - IQreceive(ii) );
    [ m, location ] = min( distances );
    index(ii) = location;
end

% -------------------------------------------------------------------------
% De-map the bits and check for errors
% -------------------------------------------------------------------------

bitsArray = [ dec2bin( index - 1, log2M ) - '0' ]';
bits = bitsArray(:);
data = data;
errors = sum( bits ~= data );
BER = errors / (N * log2M)

% -------------------------------------------------------------------------
% Convert to ASCII to double-check your results on special data sets
% -------------------------------------------------------------------------

if decode_8psk_string,
    bits8 = reshape( bits( 1:8*36 ), 8, 36 )';
    outvals = (bits8 * ( 2.^[7:-1:0]' ))';
    outstring = char( outvals )
end