%雙層次相位序列優化算法

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;
    l1=length(z1);    
     % divided into v sets and v sets divided into b sets
     z=zeros(d,512);
     x2=zeros(256,128);  
     x4=zeros(v,128);
     x6=zeros(256,512);
     x8=zeros(d,512);
     x10=zeros(1,512);
     w=zeros(1,128/d);
     b=ones(1,v/d);
     c=ones(1,d);
     for k=1:d;
         q=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,(l1/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);
NN=0:.1:13;
 for n=1:100;
  for l=1:131;
      if papr3(n)>NN(l);
      ccdf0(l)=ccdf0(l)+1;
      end  
  end
 end
 ccdf1=ccdf0./100;
figure(1)
 semilogy(NN,ccdf1,'r')
% title('CCDF曲線比較')
xlabel('PAPR( in dB)'),ylabel('CCDF Pr(PAPR>PAPR0')
legend('雙層次相位');
grid on