close all;
N = 120; % number of carriers
K = 1; %
S = 0; % start
h = modem.qammod(16); % create an object with default properties
x = randint(N, 1, 16); % generate input bits
y = modulate(h, x); % modulate x
% scatterplot(y)

buf1 = zeros(1536, 1);
buf1(169+S:168+S+N) = y; % subcarrier mapping
plot(abs(buf1));
dftout = fft(y);
plot(abs(dftout));

buf = zeros(4096, 1);
buf(1281+S:K:1280+S+N*K) = dftout; % subcarrier mapping
figure; plot(abs(buf));

sc_fdma = ifft(buf);
figure; plot(abs(sc_fdma));

papr = 20*log10(max(abs(sc_fdma))/mean(abs(sc_fdma)))

%function paprSCFDMA()

dataType = 'Q-PSK'; % Modulation format.
totalSubcarriers = 1024; % Number of total subcarriers.
numSymbols = 64; % Data block size.
Q = totalSubcarriers/numSymbols; % Bandwidth spreading factor of IFDMA.
Q_tilda = 31; % Bandwidth spreading factor of DFDMA. Q_tilda < Q.
subcarrierMapping = 'IFDMA'; % Subcarrier mapping scheme.
%subcarrierMapping = 'LFDMA'; % Subcarrier mapping scheme.
%subcarrierMapping = 'DFDMA'; % Subcarrier mapping scheme.

pulseShaping = 1; % Whether to do pulse shaping or not.
filterType = 'rc'; % Type of pulse shaping filter.
rolloffFactor = 0.0999999999; %Rolloff factor for the raised-cosine filter.
                              %To prevent divide-by-zero, for example, use 0.099999999 instead of 0.1.
Fs = 5e6; % System bandwidth.
Ts = 1/Fs; % System sampling rate.
Nos = 4; % Oversampling factor.
if filterType == 'rc' % Raised-cosine filter.
    psFilter = rcPulse(Ts, Nos, rolloffFactor);
elseif filterType == 'rr' % Root raised-cosine filter.
    psFilter = rrcPulse(Ts, Nos, rolloffFactor);
end
numRuns = 1e4; % Number of iterations.

papr = zeros(1,numRuns); % Initialize the PAPR results.

for n = 1:numRuns,
    % Generate random data.
    if dataType == 'Q-PSK'
        tmp0 = round(rand(numSymbols,2));
        tmp0 = tmp0*2 - 1;
        tmp1 = round(rand(numSymbols,2));
        tmp1 = tmp1*2 - 1;
        tmp2 = round(rand(numSymbols,2));
        tmp2 = tmp2*2 - 1;        
        tmp3 = round(rand(numSymbols,2));
        tmp3 = tmp3*2 - 1;        
        data0 = (tmp0(:,1) + j*tmp0(:,2))/sqrt(2);
        data1 = (tmp1(:,1) + j*tmp1(:,2))/sqrt(2);
        data2 = (tmp2(:,1) + j*tmp2(:,2))/sqrt(2);
        data3 = (tmp3(:,1) + j*tmp3(:,2))/sqrt(2);
    elseif dataType == '16QAM'
        dataSet = [-3+3i -1+3i 1+3i 3+3i ...
            -3+i -1+i 1+i 3+i ...
            -3-i -1-i 1-i 3-i ...
            -3-3i -1-3i 1-3i 3-3i];
        dataSet = dataSet / sqrt(mean(abs(dataSet).^2));
        tmp = ceil(rand(numSymbols,1)*16);
        for k = 1:numSymbols,
            if tmp(k) == 0
                tmp(k) = 1;
            end
            data0(k) = dataSet(tmp(k));
        end
        data0 = data0.';
    end

    % Convert data to frequency domain.
    X0 = fft(data0);                         % 1 абонент
    X1 = fft(data1);                         % 2 абонент
    X2 = fft(data2);                         % 3 абонент
    X3 = fft(data3);                         % 4 абонент

    % Initialize the subcarriers.
    Y = zeros(totalSubcarriers,1);    
    
    % Subcarrier mapping.
    if subcarrierMapping == 'IFDMA'
        Y(1:Q:totalSubcarriers) = X0;        
%        Y(2:Q:totalSubcarriers) = X1;
%        Y(3:Q:totalSubcarriers) = X2;
%        Y(4:Q:totalSubcarriers) = X3;
    elseif subcarrierMapping == 'LFDMA'
        Y(1:numSymbols) = X0;
    elseif subcarrierMapping == 'DFDMA'
        Y(1:Q_tilda:Q_tilda*numSymbols) = X0;
    end

    % Convert data back to time domain.
    y = ifft(Y);

    % Perform pulse shaping.
    if pulseShaping == 1
        % Up-sample the symbols.
        y_oversampled(1:Nos:Nos*totalSubcarriers) = y;
        % Perform filtering.
        y_result = filter(psFilter, 1, y_oversampled);
    else
        y_result = y;
    end
    
    % Calculate the PAPR.
    papr(n) = 10*log10(max(abs(y_result).^2) / mean(abs(y_result).^2));
end

% Plot CCDF.
[N,X] = hist(papr, 100);
semilogy(X,1-cumsum(N)/max(cumsum(N)),'b')

% Save data.
save paprSCFDMA