%雙層次相位序列優化算法與PTS的比較

clear;
clc;
v=16;
d=4;
nl=100;

for n=1:nl;
    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;
    x1=abs(z2).^2;
    v3=max(x1);
    m3=mean(x1);
    papr4(n)=10*log10(v3/m3);
    l1=length(z1);    
    % divided into v sets and v sets divided into b sets
    z=zeros(d,512);
    x2=zeros(256,32);  
    x4=zeros(v,32);
    x6=zeros(256,512);
    x8=zeros(d,512);
    x10=zeros(1,512);
    x12=zeros(256,512);
    x14=zeros(256,512);
    w=zeros(1,128/d);
    b1=ones(1,v);
    b=ones(1,v/d);
    c=ones(1,d);
    q=zeros(v/d,d);
     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
  
    for k=1: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/d-(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);
     %   papr00(n)=papr0(n);
     %次優化搜索
     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);
end

ccdf0=zeros(1,131);
ccdf2=zeros(1,131);
ccdf4=zeros(1,131);
ccdf6=zeros(1,131);
NN=0:.1:13;
 for n=1:nl;
  for l=1:131;
      if papr4(n)>NN(l);
      ccdf0(l)=ccdf0(l)+1;
      end  
      if papr3(n)>NN(l);
      ccdf2(l)=ccdf2(l)+1;
      end  
      if papr1(n)>NN(l);
      ccdf4(l)=ccdf4(l)+1;
      end  
      if papr5(n)>NN(l);
      ccdf6(l)=ccdf6(l)+1;
      end  
  end
 end
 ccdf1=ccdf0./nl;
 ccdf3=ccdf2./nl;
 ccdf5=ccdf4./nl;
 ccdf7=ccdf6./nl;
figure(1)
 semilogy(NN,ccdf1,'gO',NN,ccdf3,'cX',NN,ccdf5,'b*',NN,ccdf7,'ms')
 title('CCDF曲線比較')
xlabel('PAPR( in dB)'),ylabel('CCDF Pr(PAPR>PAPR0')
legend('the original signal','雙層次相位','qqq雙層次相位','PTS次優化');
grid on