%this programme is used to make ccdf by using pts
%it need function of bizdec bizinter corefrmod2 fracF f_bin2gray f_DeMQAM
%f_EnGray f_EnMQAM f_GetPoints f_randint
%V=1
clear all;
close all;
papr01=5; %the first threshold
ya=zeros(1,128);
x13=zeros(1,128);
zz12=zeros(1,128);
x=[1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5]
for m=1:10
    w=0;
    r=m
    for n=1:8000;
       ya(1:128)=randint(1,128,4);%128 Qpsk
        x13(1,1:128)=f_EnMQAM(ya(1:128),4);
        zz12(1,1:128)=fft(x13(1:128));
        p03=mean(zz12.*conj(zz12));%compute mean power
        p04=max(zz12.*conj(zz12));%compute peak power
        papr(n)=10*log10(p04/p03);%compute ratio of peak to mean power
        if papr(n)>papr01;
            w=w+1;
        end
    end
    abc(m)=w/1000;%compute ratio of papr>papr0
    papr01=papr01+0.5;%add threshold of papr0 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%V=2
papr04=5; %the first threshold
ya=zeros(1,128);
x43=zeros(1,128);
x430=zeros(1,128);
x431=zeros(1,128);
zz440=zeros(1,128);
zz441=zeros(1,128);
zz442=zeros(4,128);
%zz401=zeros(4,128);
zzz4=zeros(4,128);
cccc=[1 -1 j -j];
for l=1:10
    w=0;
    r=l
    for n=1:8000
        ya=randint(1,128,4);%128 Qpsk
        x43=f_EnMQAM(ya(1:128),4);
        x430(1:64)=x43(1:64);
        x431(65:128)=x43(65:128);
        zz440(1,1:128)=fft(x430(1:128));
        zz441(1,1:128)=fft(x431(1:128));
       for p=1:4
            zz442(p,1:128)=cccc(p)*zz441;
            zzz4(p,1:128)=zz440(1,1:128)+zz442(p,1:128);
            p042(p)=mean(zzz4(p,1:128).*conj(zzz4(p,1:128)));%compute mean power
            p043(p)=max(zzz4(p,1:128).*conj(zzz4(p,1:128)));%compute peak power
        end
    papr24=10*log10(p043./p042);%compute ratio of peak to mean power
    papr246(n)=min(papr24);
    if papr246(n)>papr04;
        w=w+1;
    end
    end
hij(l)=w/8000;%compute ratio of papr>papr0
papr04=papr04+0.5;%add threshold of papr0 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% V=3
papr02=5; %the first threshold
ya=zeros(1,128);
x23=zeros(1,128);
x119=zeros(1,128);
x230=zeros(1,128);
x231=zeros(1,128);
zz220=zeros(1,128);
zz221=zeros(1,128);
zz222=zeros(4,128);
zz2201=zeros(4,128);
zzz2=zeros(16,128);
aaaa=[1 -1 j -j];
bbbb=[1 -1 j -j];
for l=1:10
    w=0;
    r=l
    for n=1:8000
        ya=randint(1,128,4);%128 Qpsk
        x23=f_EnMQAM(ya(1:128),4);
        x119(1:43)=x23(1:43);
        x230(44:86)=x23(44:86);
        x231(87:128)=x23(87:128);
        zz119(1,1:128)=fft(x119(1:128));
         zz220(1,1:128)=fft(x230(1:128));
         zz221(1,1:128)=fft(x231(1:128));
        t=0;
         for p=1:4
             for k=1:4
              t=t+1;    
              zz222(k,1:128)=aaaa(k)*zz221;
              zz2201(p,1:128)=bbbb(p)*zz220;
              zzz2(t,1:128)=zz2201(p,1:128)+zz222(k,1:128)+zz119(1,1:128);
              p032(t)=mean(zzz2(t,1:128).*conj(zzz2(t,1:128)));%compute mean power
              p022(t)=max(zzz2(t,1:128).*conj(zzz2(t,1:128)));%compute peak power
              end
          end
        papr20=10*log10(p022./p032);%compute ratio of peak to mean power
        papr2(n)=min(papr20);
        if papr2(n)>papr02;
            w=w+1;
        end
    end
        klm(l)=w/8000;%compute ratio of papr>papr0
    papr02=papr02+0.5;%add threshold of papr0     
end 
figure(2)
z=[5 5.5 6 6.5 7 7.5 8 8.5 9 9.5]
semilogy(z,abc,'Ko-',z,hij,'Ks-',z,klm,'K+-')
legend('a=-0.5\pi,V=1','a=-0.5\pi,V=2','a=-0.5\pi,V=3','location','SouthWest')
axis([1,12,0.0001,1])
grid on
xlabel('PAPR0(dB)');
ylabel('Pr(PAPR>PAPR0)');
%title(' FFT,0.02\piFRFT-PTS-OFDM系統(tǒng)的PAPR性能曲線')