%***************************************************
% SLP-OFDM系統PARP抑制算法仿真，BPSK調制
% N 子載波數  
% M 隨機相位序列個數
% K 輸入序列數
% X 發送符號矩陣，N*K維
% P 隨機相位矩陣，N*M維
%---------------------------------------------------
clc; clear;
echo off all;
t = cputime;
N = 128; M = 32; K = 2^15;
X = 2*randint(N,K) - ones(N,K); %隨機產生K個發送符號(N維列向量)
P = 2*randint(N,M) - ones(N,M); %隨機產生M個相位序列(N維列向量)
Y = zeros(N*4,M);   %4倍過采樣
papr = zeros(1,K);
for i = 1:K      %對每個發送符號序列循環
    for j = 1:M
        Y(1:N,j) = P(:,j).*X(:,i);  %相位旋轉
    end
    y = abs(ifft(Y)).^2;  %對列向量求IFFT
    tempv = 10*log10(max(abs(y))./mean(y)); 
    papr1(i) = min(tempv(1));    %求M=1個序列中PAPR最小的值
    papr2(i) = min(tempv(1:2));  %求M=2個序列中PAPR最小的值
    papr3(i) = min(tempv(1:4));  %求M=4個序列中PAPR最小的值
    papr4(i) = min(tempv(1:8));  %求M=8個序列中PAPR最小的值
    papr5(i) = min(tempv(1:16)); %求M=16個序列中PAPR最小的值
    papr6(i) = min(tempv);       %求M=32個序列中PAPR最小的值
    %papr7(i) = min(tempv);       %求M=64個序列中PAPR最小的值
end
[cdf1, PAPR1] = ecdf(papr1);
[cdf2, PAPR2] = ecdf(papr2);
[cdf3, PAPR3] = ecdf(papr3);
[cdf4, PAPR4] = ecdf(papr4);
[cdf5, PAPR5] = ecdf(papr5);
[cdf6, PAPR6] = ecdf(papr6);
%[cdf7, PAPR7] = ecdf(papr7);

figure(1);
semilogy(PAPR1,1-cdf1,'b',PAPR2,1-cdf2,'r',PAPR3,1-cdf3,'c',PAPR4,1-cdf4,'k',PAPR5,1-cdf5,'g',PAPR6,1-cdf6,'m')
legend('M=1','M=2','M=4','M=8','M=16','M=32');
xlabel('峰均功率比 PAPR0 [dB] ');
ylabel('CCDF (Pr[PAPR>PAPR0])');
e = cputime - t;
echo on all;