clear all;
close all;

Tu=224e-6; %useful OFDM symbol period
T=Tu/2048; %baseband elementary period
G=0; %choice of 1/4, 1/8, 1/16, and 1/32
delta=G*Tu; %guard band duration
Ts=delta+Tu; %total OFDM symbol period
Kmax=1023; %number of subcarriers=1024
Kmin=0;
FS=4096; %IFFT/FFT length
q=10; %carrier period to elementary period ratio
fc=q*1/T; %carrier frequency
%Rs=4*fc; %simulation period
%t=0:1/Rs:Tu;
num=500;%number of symbols
threshold=0:.1:13;%threshold of CCDF
ccdf0=zeros(1,length(threshold));
ccdf2=ccdf0;
ccdf4=ccdf0;
u=[1 2];


%=============================================
for symbol_index=1:num
   %Data generator (A)
   M=Kmax+1;
   %rand('state',0);
   x(:,1)=randsrc(M,1,[1,-1,3,-3]);
   x(:,2)=randsrc(M,1,[1,-1,3,-3]);
   y=squeeze(x);%generate 16-QAM symbols
   A=amodce(y,fc,'qam');
   la=length(A);
   info=zeros(FS,1);
   info(1:(la/2)) = [ A(1:(la/2)).']; %Zero padding
   info((FS-((la/2)-1)):FS) = [ A(((la/2)+1):la).'];%Zero padding 
   info=info';

   %Subcarriers generation (B)
   carriers0=FS.*ifft(info,FS);
   lc=length(carriers0);
   %tt=0:T/2:Tu;


   %figure(1);
   power0=abs(carriers0).^2;
   %plot(power0);
   
   
   %figure(2);
   %f=(2/T)*(1:(FS))/(FS);
   %subplot(211);
   %plot(f,abs(fft(carriers,FS))/FS);
   %subplot(212);
   %pwelch(carriers,[],[],[],2/T);


   max0=max(power0);
   mean0=mean(power0);
   papr0(symbol_index)=10*log10(max0/mean0);
   
   
   V=ones(1,lc);
   V=V*mean(abs(carriers0));
   carriers1=(V.*sign(carriers0).*log(1+u(1)*abs(carriers0)./V))/log(1+u(1));
   carriers2=(V.*sign(carriers0).*log(1+u(2)*abs(carriers0)./V))/log(1+u(2));
   
   power1=abs(carriers1).^2;
   max1=max(power1);
   mean1=mean(power1);
   papr1(symbol_index)=10*log10(max1/mean1);
   
   power2=abs(carriers2).^2;
   max2=max(power2);
   mean2=mean(power2);
   papr2(symbol_index)=10*log10(max2/mean2);
   
   
  for l=1:length(threshold)
       if papr0(symbol_index)>threshold(l)
           ccdf0(l)=ccdf0(l)+1;
       end
       if papr1(symbol_index)>threshold(l)
           ccdf2(l)=ccdf2(l)+1;
       end
       if papr2(symbol_index)>threshold(l)
           ccdf4(l)=ccdf4(l)+1;
       end
   end
end   


ccdf1=ccdf0./num;
ccdf3=ccdf2./num;
ccdf5=ccdf4./num;
semilogy(threshold,ccdf1,'r',threshold,ccdf3,'b',threshold,ccdf5,'g')
title('壓擴變換方法的CCDF曲線')
xlabel('papr(dB)'),ylabel('ccdf')
grid on