%雙層次相位序列優化算法（不同的分塊取值）

clear;
clc;
v=16;
v1=8;
nl=100;

for n=1:nl;
    k2=1;
    b1=ones(1,v);
    b2=ones(1,v1);
    z3=zeros(v1,512);
    iz1=zeros(v1,512);
    x12=zeros(1,512);
    x14=zeros(v,512);
    x16=zeros(1,512);
    x18=zeros(v1,512);
    
    a=rand(1,256)>0.5;
    C=(a*2-1)/sqrt(2);
    sreal=C(1:2:256);
    simage=C(2:2:256);
    D=sreal+j*simage; 
    % oversample 4
    z1=[D(1:128/2),zeros(1,3*128),D(128/2+1:128)];
    z2=ifft(z1)*512;
    x0=abs(z2).^2;
    v3=max(x0);
    m3=mean(x0);
    papr4(n)=10*log10(v3/m3);
    papr4(n)=papr4(n)+3;
    l1=length(z1);    
    
    z3(1,:)=[D(1:16),zeros(1,496)];
    z3(2,:)=[zeros(1,16),D(17:32),zeros(1,480)];
    z3(3,:)=[zeros(1,32),D(33:48),zeros(1,464)];
    z3(4,:)=[zeros(1,48),D(49:64),zeros(1,448)];
    z3(5,:)=[zeros(1,64),D(65:80),zeros(1,432)];
    z3(6,:)=[zeros(1,80),D(81:96),zeros(1,416)];
    z3(7,:)=[zeros(1,96),D(97:102),zeros(1,424)];
    z3(8,:)=[zeros(1,102),D(103:1228),zeros(1,496)];
    
    % divided into v sets
    for k1=1:v;
        z(k1,:)=[zeros(1,(k1-1)*128/v),D((k1-1)*128/v+1:k1*128/v),zeros(1,(l1-(k1-1)*128/v-128/v))];
        iz(k1,:)=ifft(z(k1,:));  
        iz(k1,:)=iz(k1,:)*512;
        x12(1,:)=x12(1,:)+b1(k1)*iz(k1,:);
    end
    x13=abs(x12).^2;
    m10=mean(x13);
    v10=max(x13);
    papr5(n)=10*log10(v10/m10);
   
  %次優化搜索
    index=1; 
    while index<v+1
          b1(index)=-1;
          for m=1:v;
              x14(index,:)=x14(index,:)+b1(m)*iz(m,:);
              m+1;
          end
          x15=abs(x14(index,:)).^2;
          m1=mean(x15);
          v1=max(x15);
          papr0000(n)=10*log10(v1/m1);
          if papr0000(n)<papr5(n)       
             papr5(n)=papr0000(n);
          else
             b1(index)=1;
          end  
          index=index+1;
    end

    % divided into v1 sets
    z3(1,:)=[D(1,16),zeros(17,512)];
    while k2<v1+1
        z3(k2,:)=[zeros(1,(k2-1)*128/v1),D((k2-1)*128/v1+1:k2*128/v1),zeros(1,(l1-(k2-1)*128/v1-128/v1))];
        iz1(k2,:)=ifft(z3(k2,:));  
        iz1(k2,:)=iz1(k2,:)*512;
        x16(1,:)=x16(1,:)+b2(k2)*iz1(k2,:);
        k2=k2+1;
    end
    x17=abs(x16).^2;
    m10=mean(x17);
    v10=max(x17);
    papr10(n)=10*log10(v10/m10);
   
  %次優化搜索
    index=1; 
    while index<v1+1
          b2(index)=-1;
          for m=1:v1;
              x18(index,:)=x18(index,:)+b1(m)*iz1(m,:);
              m+1;
          end
          x19=abs(x18(index,:)).^2;
          m1=mean(x19);
          v1=max(x19);
          papr100(n)=10*log10(v1/m1);
          if papr100(n)<papr10(n)       
             papr10(n)=papr100(n);
          else
             b2(index)=1;
          end  
          index=index+1;
    end
    
     % divided into v sets and v sets divided into b sets  
    t=1;
    d=2;  % 2,4,8
    while d<9          

          x2=zeros(1,128);  
          x4=zeros(v,128);
          x6=zeros(1,l1);
          x8=zeros(d,l1);
          x10=zeros(1,l1);
          w=zeros(1,128/d);
          b=ones(1,v/d);
          c=ones(1,d);
          for k=1:d;
              q=zeros(v/d,d);
              p=zeros(v/d,d);
              z(k,:)=[zeros(1,(k-1)*128/d),D((k-1)*128/d+1:k*128/d),zeros(1,(l1-(k-1)*128/d-128/d))];
              w=z(k,(k-1)*128/d+1:k*128/d);%長度32
              for m=1:v/d;
                  s(m,:)=[zeros(1,(m-1)*128/v),w((m-1)*128/v+1:m*128/v),zeros(1,(128-(m-1)*128/v-128/v))];
                  is(m,:)=ifft(s(m,:));  
                  is(m,:)=is(m,:)*128;
                  x2(1,:)=x2(1,:)+b(m)*is(m,:);
              end
              x3=abs(x2).^2;
              m0=mean(x3);
              v0=max(x3);
              papr0(n)=10*log10(v0/m0);
         
              index=1; 
              while index<v/d+1
                    b(index)=-1;
                    for m=1:v/d;
                        x4(index,:)=x4(index,:)+b(m)*is(m,:);
                        m+1;
                    end
                    x5=abs(x4(index,:)).^2;
                    m1=mean(x5);
                    v1=max(x5);
                    papr1(n)=10*log10(v1/m1);
                    if papr1(n)<papr0(n)       
                       papr0(n)=papr1(n);
                       q(index,k)=-1;
                    else
                    b(index)=1;
                    q(index,k)=1;
                    end  
                    index=index+1;
              end   
              iz(k,:)=ifft(z(k,:));  
              iz(k,:)=iz(k,:)*512;
              x6(1,:)=x6(1,:)+c(k)*iz(k,:);
          end
          x7=abs(x6).^2;
          m1=mean(x7);
          v1=max(x7);
          papr00(n)=10*log10(v1/m1);

          %次優化搜索
          index=1; 
          while index<d+1
                c(index)=-1;
                for m=1:d;
                    x8(index,:)=x8(index,:)+c(m)*iz(m,:);
                    m+1;
                end
                x9=abs(x8(index,:)).^2;
                m1=mean(x9);
                v1=max(x9);
                papr2(n)=10*log10(v1/m1);
                if papr2(n)<papr00(n)       
                   papr00(n)=papr2(n);
                   c(1,index)=-1;
                else
                c(index)=1;
                c(1,index)=-1;
                end  
                index=index+1;
          end
     
          %矩陣相乘再計算峰均比  
          p(:,k)=q(:,k)*c(k);
          P=reshape(p,1,v);
          % oversample 4
          for w=1:v;
              z3(w,:)=[zeros(1,(w-1)*128/v),D((w-1)*128/v+1:w*128/v),zeros(1,(l1-(w-1)*128/v-128/v))];
              iz(w,:)=ifft(z3(w,:));  
              iz(w,:)=iz(w,:)*512;
              x10(1,:)=x10(1,:)+P(w)*iz(w,:);
          end
          x11=abs(x10).^2;
          v3=max(x11);
          m3=mean(x11);
          papr3(n)=10*log10(v3/m3);
          if t==1
              papr7(n)=papr3(n);
          elseif t==2
              papr8(n)=papr3(n);
          else
              papr9(n)=papr3(n);
          end
          t=t+1;
          d=2*d;
    end
end

ccdf0=zeros(1,131);
ccdf2=zeros(1,131);
ccdf4=zeros(1,131);
ccdf6=zeros(1,131);
ccdf8=zeros(1,131);
ccdf10=zeros(1,131);
NN=0:.1:13;
 for n=1:nl
     for l=1:131;
         if papr5(n)>NN(l);
            ccdf0(l)=ccdf0(l)+1;
         end  
         if papr4(n)>NN(l);
            ccdf2(l)=ccdf2(l)+1;
         end  
         if papr7(n)>NN(l);
            ccdf4(l)=ccdf4(l)+1;
         end
         if papr8(n)>NN(l);
            ccdf6(l)=ccdf6(l)+1;
         end
         if papr9(n)>NN(l);
            ccdf8(l)=ccdf8(l)+1;
         end     
         if papr10(n)>NN(l);
            ccdf10(l)=ccdf10(l)+1;
         end     
     end
 end
 ccdf1=ccdf0./nl;%迭代 
 ccdf3=ccdf2./nl;%
 ccdf5=ccdf4./nl;%d=2
 ccdf7=ccdf6./nl;%d=4
 ccdf9=ccdf8./nl;%d=8
 ccdf11=ccdf10./nl;%d=8
figure(1)
semilogy(NN,ccdf1,'gO-',NN,ccdf3,'k-',NN,ccdf5,'cX-',NN,ccdf7,'rs-',NN,ccdf9,'b*-',NN,ccdf10,'m+-')
title('CCDF曲線比較')
xlabel('PAPR( in dB)'),ylabel('CCDF Pr(PAPR>PAPR0')
legend('PTS次優化算法v=16','原始曲線','雙層次相位序列優化v=16 d=2','v=16 d=4','v=16 d=8');
grid on