src code

BerOFDM.m

+function [i,y] = BerOFDM(n)
+    len = length(n);
+    ber = zeros(1,len);
+    index = zeros(1,len);
+    for i=0:1:n
+        BER=process_ber(6,i);
+        ber(i+1)=BER;
+        index(i+1)=i;
+    end
+    y=ber;
+    x2=linspace(min(index),max(index));
+    y2=interp1(index,ber,x2,'cubic');
+%     semilogy(index,ber);
+    semilogy(x2,y2);
+    i=index;
+end
+

BerPlot.m

+function [i,y] = BerPlot(n,method)
+    len = length(n);
+    ber = zeros(1,len);
+    index = zeros(1,len);
+    for i=0:1:n
+        [BER,v1,v2,v3,p]=process(6,i,method);
+%         ber(int64(i*500)+1)=BER;
+%         index(int64(i*500)+1)=i;
+        ber(i+1)=BER;
+        index(i+1)=i;
+    end
+    y=ber;
+    x2=linspace(min(index),max(index));
+    y2=interp1(index,ber,x2,'cubic');
+%     semilogy(index,ber);
+    semilogy(x2,y2);
+    hold on;
+    i=index;
+end
+

HPA_3.m

+function [ya,yp] = HPA_3(xa,xp)
+%DHPA Summary of this function goes here
+%   Detailed explanation goes here
+
+% 
+% ya=xa;
+% yp=xp;
+ya=(6*xa)./(xa.^2+3);
+yp=(pi/3)*((xa.^2)./(xa.^2+1))+xp;

SCFDE_ccdf.m

+len=50000;
+p_scfde=zeros(1,len);
+for n=1:len
+    x_scfde=abs(scfde());
+    p_scfde(n)=10*log10(max(abs(x_scfde)).^2/(mean(abs(x_scfde).^2)));
+end
+[papr_scfde,prob_scfde]=ccdf1(p_scfde,len);
+figure(4) 
+axis([0,20,10^(-5),1])
+semilogy(papr,prob,'-o','linewidth',2,'markersize',6)
+hold on
+semilogy(papr_scfde,prob_scfde,'-d','linewidth',2,'markersize',6)
\ No newline at end of file

TOFDMtest.m

+len=50000;
+p_Tofdm=zeros(1,len);
+
+for n=1:len
+   x_Tofdm=Tofdm_signal(256);
+   xa_Tofdm=abs(x_Tofdm);
+   p_Tofdm(n)=10*log10(max(xa_Tofdm).^2/(mean(xa_Tofdm.^2)));
+end
+
+[papr_Tofdm,prob_Tofdm]=ccdf1(p_Tofdm,len);
+
+figure(1)
+semilogy(papr_Tofdm,prob_Tofdm,'-o')
+axis([0,20,10^(-5),1])
\ No newline at end of file

Tofdm_signal.m

+%function [ F,W,T,T1 ] = matrixWF( N )
+function y=Tofdm_signal(N)
+a=(0:N-1)'*(0:N-1);
+F=exp(-2*pi*1i*a/N)/sqrt(N);%DFT󡪡F
+
+W=hadamard(N);
+
+T=(1/N)*F*W;
+%T1=(1/N)*W*F';
+
+x_ofdm=ofdm_signal();
+
+x=x_ofdm'*T;
+
+y=x;
+
+%E_TOFDM=(1/N)*()
+
+end
+

Untitled.m

+plot(snr1,p1,'-v')
+hold on
+plot(snr2,p2,'-o')
+hold on
+plot(snr3,p3,'-squre')
+hold on
+plot(snr4,p4u,'-g')
\ No newline at end of file

Untitled1.m

+for n=1:length(x1)
+    if (x1(n)>0.5)&&(x1(n)<1.5)
+        x1(n)=x1(n);
+    else
+        x1(n)=0;
+    end
+end
+
+t=ze(x1);
+
+figure(1)
+subplot(3,1,1);
+plot(abs(x))
+axis([0,256,0,2])
+title('(a) Original Signal')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(3,1,2);
+plot(x1)
+axis([0,256,0,2])
+title('(b) Amplitudes of signal among the interval [r_a, r_b]')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(3,1,3)
+plot(abs(t))
+axis([0,256,0,2])
+title('(c) S_{new}')
+xlabel('t')
+ylabel('r(t)')
+

Untitled2.m

+
+n=1
+for ibo=10:-2:-10
+    [ber(n)]=tested(ibo);
+    n=n+1
+end

Untitled3.m

+% 
+
+
+% yitayita1=mean(abs(y1).^2);
+% yitayita2=mean(abs(y2).^2);
+% yitayita3=mean(abs(y3).^2);
+
+for n=1:length(y1)
+    if abs(y1(n))>1
+        y1e(n)=1;
+    else
+        y1e(n)=y1(n);
+    end
+end
+
+
+for n=1:length(y1)
+    if abs(y2(n))>1
+        y2e(n)=1;
+    else
+        y2e(n)=y2(n);
+    end
+end
+
+
+for n=1:length(y1)
+    if abs(y1(n))>1
+        y3e(n)=1;
+    else
+        y3e(n)=y3(n);
+    end
+end
+
+yitayita1=mean(abs(y1e).^2);
+yitayita2=mean(abs(y2e).^2);
+yitayita3=mean(abs(y3e).^2);
+
+m=0;
+count=0;
+for n=1:length(y3)
+    if abs(y3e(n))>0.0001
+        m=m+abs(y3(n))^2;
+        count=count+1;
+    end
+end
+
+m/count
+% 
+% subplot(3,1,1)
+% plot(abs(y1e))
+% subplot(3,1,2)
+% plot(abs(y2e))
+% subplot(3,1,3)
+% plot(abs(y3e))
\ No newline at end of file

a.m

+function y=a(x,r)
+
+y=sqrt((abs(x)/r)^2-1);
\ No newline at end of file

average.m

+function [y1,y2,y3,p1,p2,p3]=average(x,V1,V2)
+r=(1/(sqrt(3)))*(0.3660+0.8660*1i);
+
+for n=1:length(x)
+    if abs(x(n))<V1
+        y1(n)=x(n);
+        y2(n)=0;
+        y3(n)=0;
+    elseif (abs(x(n))<sqrt(V1^2+V2^2))&&(abs(x(n))>V1)
+        y1(n)=(1+1i)*x(n)/2;
+        y2(n)=(1-1i)*x(n)/2;
+        y3(n)=0;
+    else%if (abs(x(n))<1.732)&&(abs(x(n))>1.414)
+        y3(n)=x(n)/sqrt(3);
+        y2(n)=x(n)*(0.2113-0.5373*1i);
+        y1(n)=x(n)*(0.2113+0.5373*1i);
+    end
+end
+
+p1=funcpapr(y1);
+p2=funcpapr(y2);
+p3=funcpapr(y3);

ccdf.m

+function [papr0,pp]=ccdf(p)
+
+N=length(p);
+papr0=0:0.5:20;
+prob=zeros(1,length(papr0));
+for n=1:N
+    for papr=0:0.5:20
+        if p(n)<papr
+            prob(papr*2+1)=prob(papr*2+1)+1;
+        end
+    end
+end
+
+pp=prob/5000;
+

ccdf1.m

+function [papr0,pp]=ccdf1(p,len)
+
+papr0=0:0.25:20;
+prob=zeros(1,length(papr0));  %%41
+for n=1:len
+    for papr=0:0.25:20
+        if p(n)>papr
+            prob(papr*4+1)=prob(papr*4+1)+1;
+        end
+    end
+end
+
+pp=prob/len;
+

coutnn.m

+cy1=0;
+cy2=0;
+cy3=0;
+
+x=abs(d);
+for n=1:length(x)
+    if (x(n)<1)||(x(n)==1)
+        cy1=cy1+1;
+    elseif (x(n)<sqrt(2))||(x(n)==sqrt(2))
+       cy2=cy2+1;
+    else
+        cy3=cy3+1;
+    end
+end
\ No newline at end of file

cubic.m

+function [y1,y2,y3,p1,p2,p3]=cubic(x,v1,v2)
+
+for n=1:length(x)
+    if abs(x(n))<v1
+        y1(n)=x(n);
+        y2(n)=0;
+        y3(n)=0;
+    elseif (abs(x(n))<sqrt(v1^2+v2^2))&&(abs(x(n))>v1)
+        y1(n)=x(n)/(1+a(x(n),1)*1i);
+        y2(n)=x(n)-y1(n);
+        y3(n)=0;
+    else%if (abs(x(n))<1.732)&&(abs(x(n))>1.414)
+        yy1(n)=x(n)/(1+a(x(n),sqrt(2))*1i);
+        yy2(n)=x(n)-yy1(n);
+
+        y1(n)=yy1(n)/(1+a(yy1(n),1)*1i);
+        y2(n)=yy1(n)-y1(n);
+        y3(n)=yy2(n);
+    end
+end
+
+% g1=ze(y1);
+% g2=ze(y2);
+% g3=ze(y3);
+p1=funcpapr(y1);
+p2=funcpapr(y2);
+p3=funcpapr(y3);

cubicNcout.m

+function [ power1,power2,power3,n1,n2,n3 ] = cubicNcout( x,v1,v2 )
+    n1=0;
+    n2=0;
+    n3=0;
+    for n=1:length(x)
+        if abs(x(n))<v1
+            power1(n)=abs(x(n))^2;
+            n1=n1+1;
+        elseif (abs(x(n))<sqrt(v1^2+v2^2))&&(abs(x(n))>v1)
+            power2(n)=abs(x(n))^2;
+            n2=n2+1;
+        else%if (abs(x(n))<1.732)&&(abs(x(n))>1.414)
+            power3(n)=abs(x(n))^2;
+            n3=n3+1;
+        end
+    end
+
+end
+
+
+
+

cubicPAPR.m

+len=50000;
+
+for n=1:len
+    x=ofdm_signal();
+    xa=abs(x);
+    p(n)=10*log10(max(xa).^2/(mean(xa.^2)));%%PAPRļ
+end
+
+pow=mean(xa.^2);
+
+[papr,prob]=ccdf1(p,len);
+
+% 
+% x=ofdm_signal();
+% xa=abs(x);
+% pow=mean(xa.^2);
+f1=(1-exp(-1/pow));
+f2=(exp(-1/pow)-exp(-2/pow));
+f3=exp(-2/pow);
+p1=pow-(1+pow)*exp(-1/pow);
+p2=(1+pow)*exp(-1/pow)-(2+pow)*exp(-2/pow);
+p3=(2+pow)*exp(-2/pow);
+
+papr_cubic_hds=10*log10(1/p1+f2+f3);
+papr_cubic_mds=10*log10(f2+f3/p2+f3-f2);
+papr_ofdm=10*log10(1/pow);
+deltaPAPR_hds=papr_cubic_hds-papr_ofdm;
+deltaPAPR_mds=papr_cubic_mds-papr_ofdm;
+
+
+figure(1)
+hold on
+semilogy(papr-papr_cubic_hds,prob,'o')
+hold on
+semilogy(papr-papr_cubic_mds,prob,'*')
+hold on
+% semilogy(papr-papr3,prob,'s')
+% hold on 
+%semilogy(ofdmpr+3,prob1,'d')
+semilogy(papr,prob,'d')
+
+axis([0,20,10^(-4),1]);
\ No newline at end of file

cubic_delta.m

+function [y1,y2,y3,p1,p2,p3]=cubic_delta(x)
+d=0.2;
+R1=1+d;
+R2=1+d;
+R3=1+d;
+for n=1:length(x)
+    if abs(x(n))<R1
+        y1(n)=x(n);
+        y2(n)=0;
+        y3(n)=0;
+    elseif (abs(x(n))<sqrt(R1^2+R2^2))&&(abs(x(n))>R2)
+        y1(n)=x(n)/(1+a(x(n),R1)*1i);
+        y2(n)=x(n)-y1(n);
+        y3(n)=0;
+    else%if (abs(x(n))<1.732)&&(abs(x(n))>1.414)
+        yy1(n)=x(n)/(1+a(x(n),sqrt(R1^2+R2^2))*1i);
+        yy2(n)=x(n)-yy1(n);
+
+        y1(n)=yy1(n)/(1+a(yy1(n),R1)*1i);
+        y2(n)=yy1(n)-y1(n);
+        y3(n)=yy2(n);
+    end
+end
+
+% g1=ze(y1);
+% g2=ze(y2);
+% g3=ze(y3);
+p1=funcpapr(y1);
+p2=funcpapr(y2);
+p3=funcpapr(y3);

cubicmethod.m

+% function [z,BER,SER,e]=cubicmethod(ibo,snr)
+tic
+clear
+snr='none';
+ibo=4.1;
+method='cubic';
+
+
+
+dd=1;
+
+
+
+loop=5000;
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+cap=N;
+num=sqrt((10/dd)*10^(ibo/10));
+aa=0.1;
+bb=0.1;
+
+Tab_a_a=ones(1,cap); %时域预失真幅度查询表
+Tab_p_a=zeros(1,cap); %时域预失真相位查询表
+address_a=zeros(1,N); %索引表
+
+Tab_a_b=ones(1,cap); %时域预失真幅度查询表
+Tab_p_b=zeros(1,cap); %时域预失真相位查询表
+address_b=zeros(1,N); %索引表
+
+Tab_a_c=ones(1,cap); %时域预失真幅度查询表
+Tab_p_c=zeros(1,cap); %时域预失真相位查询表
+address_c=zeros(1,N); %索引表
+%--------------------------------------------以上为策略定义,之后为增加速度定义
+Ina1_a=ones(1,N);
+Outa_a=ones(1,N);
+Outp_a=zeros(1,N);
+Ina1_b=ones(1,N);
+Outa_b=ones(1,N);
+Outp_b=zeros(1,N);
+Ina1_c=ones(1,N);
+Outa_c=ones(1,N);
+Outp_c=zeros(1,N);
+%------------------------------------------------
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=addzeros(xo,dd);
+    d=ifft(xT)*sqrt(N);
+    %------------------------------------------
+    switch method
+        case 'average'
+            [y1,y2,y3,p1,p2,p3]=average(d);
+        case 'cubic'
+            [y1,y2,y3,p1,p2,p3]=cubic(d);
+    end
+    %------------------------------------------
+    %--------------------the first layer
+    Ina_a=abs(y1);
+    Inp_a=angle(y1);
+    for uj=1:length(Ina_a)                          %时域预失真
+        [address_a(uj),Ina1_a(uj)]=Q(Ina_a(uj),1,cap);% Q: adrress function
+        Outa_a(uj)=Ina1_a(uj)*Tab_a_a(address_a(uj));
+        Outp_a(uj)=Inp_a(uj)+Tab_p_a(address_a(uj));
+    end
+    
+    [ya_a,yp_a]=HPA(Outa_a,Outp_a);
+    y_a=p2c(ya_a,(yp_a));
+    
+    for ui=1:length(address_a)    %时域预失真更新公式
+        Tab_a_a(address_a(ui))=Tab_a_a(address_a(ui))-aa*(ya_a(ui)-Ina1_a(ui));
+        Tab_p_a(address_a(ui))=Tab_p_a(address_a(ui))-bb*Ep(yp_a(ui)-Inp_a(ui));
+    end
+    %-----------------------the second layer
+    Ina_b=abs(y2);
+    Inp_b=angle(y2);
+    for uj=1:length(Ina_b)                          %时域预失真
+        [address_b(uj),Ina1_b(uj)]=Q(Ina_b(uj),1,cap);% Q: adrress function
+        Outa_b(uj)=Ina1_b(uj)*Tab_a_b(address_b(uj));
+        Outp_b(uj)=Inp_b(uj)+Tab_p_b(address_b(uj));
+    end
+    
+    [ya_b,yp_b]=HPA(Outa_b,Outp_b);
+    y_b=p2c(ya_b,(yp_b));
+    
+    for ui=1:length(address_b)    %时域预失真更新公式
+        Tab_a_b(address_b(ui))=Tab_a_b(address_b(ui))-aa*(ya_b(ui)-Ina1_b(ui));
+        Tab_p_b(address_b(ui))=Tab_p_b(address_b(ui))-bb*Ep(yp_b(ui)-Inp_b(ui));
+    end
+    %-----------------------the third layer
+    Ina_c=abs(y3);
+    Inp_c=angle(y3);
+    for uj=1:length(Ina_c)                          %时域预失真
+        [address_c(uj),Ina1_c(uj)]=Q(Ina_c(uj),1,cap);% Q: adrress function
+        Outa_c(uj)=Ina1_c(uj)*Tab_a_c(address_c(uj));
+        Outp_c(uj)=Inp_c(uj)+Tab_p_c(address_c(uj));
+    end
+    
+    [ya_c,yp_c]=HPA(Outa_c,(Outp_c));
+    y_c=p2c(ya_c,yp_c);
+    
+    for ui=1:length(address_c)    %时域预失真更新公式
+        Tab_a_c(address_c(ui))=Tab_a_c(address_c(ui))-aa*(ya_c(ui)-Ina1_c(ui));
+        Tab_p_c(address_c(ui))=Tab_p_c(address_c(ui))-bb*Ep(yp_c(ui)-Inp_c(ui));
+    end
+    %-----------------------------------------------开始做反馈
+    y=y_a+y_b+y_c;
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+    
+    Go=Yreceive;     %反馈
+    
+    Gone=fft(Go)/sqrt(N);  %反馈做fft
+    Gone=reshape(Gone,1,length(Gone));  %串并变换
+    % receive=Gone*num;
+    
+    Gone1=quitzero(Gone,N);
+    Result=demodulate(de,Gone1*num);
+    z=Gone1*num;
+    e(start)=mean(abs(Gone1-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+    [ses,ser(start)]=symerr(Result,data);
+end
+BER=mean(ber(loop-2500:loop));
+SER=mean(ser(loop-2500:loop));
+figure(2)
+plot(z,'*')
+toc
\ No newline at end of file

cut.m

+function [BER,v]=cut(ibo,snr)
+% snr='none';
+% ibo=0;
+dd=1;
+% nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+
+loop=10000;
+ber=zeros(1,loop);
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+
+num=sqrt((10/dd)*10^(ibo/10));
+
+
+
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=xo;
+    d=ifft(xT)*sqrt(N);
+    Ina=abs(d);
+    Inp=angle(d);
+    
+    [ya,yp]=idear(Ina,Inp,sqrt(2));
+    y=p2c(ya,yp);
+    
+    %-----------------------------------------------ʼ
+    
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+    
+    Go=Yreceive;     %
+    
+    Gone=fft(Go)/sqrt(N);  %fft
+    Gone=reshape(Gone,1,length(Gone));  %任
+    Result=demodulate(de,Gone*num);
+    z=Gone*num;
+    e(start)=mean(abs(Gone-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+%     [ses,ser(start)]=symerr(Result,data);
+end
+BER=mean(ber);
+% SER=mean(ser(loop-2500:loop));
+v=mean(abs(y)/sqrt(3));
+% plot(z,'*')

dbplot.m

+function y=dbplot()
+    leng=500;
+    numX = zeros(1,41);
+%     xo=ofdm_signal();
+%     xoA=abs(10*log10(xo.^2));
+%     for i=1:256
+%         for n=0:1:16
+%             if (xoA(i)>=n && xoA(i)<(n+0.2))
+%                 numX(int64(n)+1)=numX(int64(n)+1)+1;
+%             end
+%         end
+%     end
+
+%% OFDM
+%     for n=1:leng
+%         xo=ofdm_signal();
+%         xoA=abs(xo);
+%         p(n)=10*log10(max(xoA).^2/(mean(xoA.^2)));%ԭʼOFDMźxPAPR
+%     end
+
+%% T-OFDM
+%     for n=1:leng
+%         xo=Tofdm_signal(256);
+%         xoA=abs(xo);
+%         p(n)=10*log10(max(xoA).^2/(mean(xoA.^2)));
+%     end
+
+%% author1
+%     for m=1:leng
+%         x=ofdm_signal();
+%         [y1,y2,y3,p1s(m),p2s(m),p3s(m)]=average(x,1,1);
+%     end
+
+%% author2
+    for m=1:leng
+        x=ofdm_signal();
+        [y1,y2,y3,p1s(m),p2s(m),p3s(m)]=cubic(x,1,1);
+    end
+
+    for i=1:leng
+        for n=4:0.2:12
+            if (p1s(i)>=n && p1s(i)<(n+0.2))
+                numX(int64((n-4)*5)+1)=numX(int64((n-4)*5)+1)+1;
+            end
+        end
+    end
+    
+    y=numX;
+end
+

e.m

+function r=e(x,p)
+
+r=exp(-x/p);
\ No newline at end of file

eta.m

+function y=prob(x1,x2,t)
+
+%y=-t*x2*exp(-(t*x2)^2)+t*x1*exp(-(t*x1)^2)+sqrt(pi)/2*(erf(t*x2)-erf(t*x1));
+
+y=t*x1*exp(-(t*x1)^2)-t*x2*exp(-(t*x2)^2)+(sqrt(pi)/2)*(erf(t*x2)-erf(t*x1));
+

excute.m

+tic
+m=1;
+BER1=zeros(5,11);
+BER2=zeros(5,11);
+BER3=zeros(5,11);
+BER4=zeros(5,11);
+BER5=zeros(5,11); 
+
+p1=zeros(5,11);
+
+p21=zeros(5,11);
+p22=zeros(5,11);
+p23=zeros(5,11);
+p2x=zeros(5,11);
+
+p31=zeros(5,11);
+p32=zeros(5,11);
+p33=zeros(5,11);
+p3x=zeros(5,11);
+
+p41=zeros(5,11);
+p42=zeros(5,11);
+p43=zeros(5,11);
+p4x=zeros(5,11);
+p4y=zeros(5,11);
+
+p5=zeros(5,11);
+
+snr5=zeros(1,5);
+snr4=zeros(1,5);
+snr3=zeros(1,5);
+snr2=zeros(1,5);
+snr1=zeros(1,5);
+%---------------------------------------------
+ber0=1E-4;
+snr='none';
+for ibo=4:-1:0
+    n=1;
+    for snr=0:2:20
+        [BER3(m,n),p31(m,n),p32(m,n),p33(m,n),p3x(m,n)]=process(ibo,snr,'average');
+        [BER2(m,n),p21(m,n),p22(m,n),p23(m,n),p2x(m,n)]=process(ibo,snr,'cubic');
+        [BER1(m,n),p1(m,n)]=tdpd(ibo,snr);
+%         [BER5(m,n),p5(m,n)]=cut(ibo,snr);
+        n=n+1;
+    end
+    snr3(m)=forSNR(BER3(m,1:11),0:2:20,ber0);
+    snr2(m)=forSNR(BER2(m,1:11),0:2:20,ber0);
+    snr1(m)=forSNR(BER1(m,1:11),0:2:20,ber0);
+%     snr5(m)=forSNR(BER5(m,1:11),0:2:20,ber0);
+    m=m+1;
+end
+% 
+% m=1;
+% for ibo=0:-1:-8
+%     n=1;
+%     for snr=0:2:20
+%         [BER4(m,n),p41(m,n),p42(m,n),p43(m,n),p4x(m,n),p4y(m,n)]=segmentmain(ibo,snr);
+%         n=n+1;
+%     end
+%     snr4(m)=forSNR(BER4(m,1:11),0:2:20,ber0);
+%     m=m+1;
+% end
+p1=(pi/4)*(mean(p1,2));
+p2=(pi/4)*(mean(p2x,2));
+p3=(pi/4)*mean(p3x,2);
+% p4=(pi/4)*(mean(p4x,2));
+% p4u=(pi/4)*(mean(p4y,2));
+% p5=(pi/4)*(mean(p5,2));
+save 20100518_curve.mat
+toc
\ No newline at end of file

f.m

+
+
+n=1;
+
+% p=1.3;
+N=256;
+p=1;
+po1=p-(1+p)*e(1,p);
+po2=(1+p)*e(1,p)-(2+p)*e(2,p);
+po3=(2+p)*e(2,p)-(3+p)*e(3,p);
+
+p1=1-exp(-1);
+p2=exp(-1)-exp(-2);
+p3=exp(-2)-exp(-3);
+
+
+papr1(n)=10*log10(3*(po1+po2/2+po3/3)/p);
+papr2(n)=10*log10(3*(po2/2+po3/3)/(p2+p3)/p);
+papr3(n)=10*log10(po3/p3/p);
+
+% y1(n)=10*log10(1/(p1+f2+f3));
+% y2(n)=10*log10((f2+f3)/(p2-f2+f3));
+% for s=0.6:0.01:0.8
+% y3(n)=10*log10(10^s/(10^s-2)*(p3/f3-2));
+% n=n+1;
+% end
+% n=n+1;
+eta1=po1+po2/2+po3/3;
+eta2=(po2/2+po3/3)/(p2+p3);
+eta3=(po3/3)/p3;
+
+yita1=po1+po2/2+po3/3;
+yita2=(po2/2+po3/3);
+yita3=(po3/3);
+
+yita=(po1+po2+po3)/3;
+
+
+
+
+

forPLOT.m

+x1=mean(p1(1:5,1:11)')';
+x2=mean(p2(1:5,1:11)')';
+xxx=10*log10(x2'./x1')-(snr2(1:5)-snr1(1:5));
+figure(2)
+plot(xxx,'-o')
\ No newline at end of file

forSNR.m

+function y=forSNR(ber,snr,x)
+
+y=interp1(ber,snr,x,'spline');
\ No newline at end of file

funcpapr.m

+%papr
+function y=funcpapr(x)
+g=ones(1,1);
+t=ones(1,1);
+xmax=max(abs(x).^2);%POWmax
+
+j=1;
+for n=1:length(x)
+    if (real(x(n))~=0)||(imag(x(n)~=0))
+        g(j)=x(n);
+        j=j+1;
+    end
+end
+
+if g==t
+    y=0;
+else
+    xmean=mean(abs(g).^2);
+    y=10*log10(xmax/xmean);
+end

huatu.m

+figure(2)
+subplot(4,1,1)
+plot(abs(x))
+axis([0,256,0,2])
+title('Original Signal')
+xlabel('t')
+ylabel('r(t)')
+
+
+subplot(4,1,2)
+plot(abs(y1))
+axis([0,256,0,1])
+title('HDS')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(4,1,3)
+plot(abs((y2)))
+axis([0,256,0,1])
+title('MDS')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(4,1,4)
+plot(abs((y3)))
+axis([0,256,0,1])
+title('LDS')
+xlabel('t')
+ylabel('r(t)')
+
+
+figure(3)
+subplot(4,1,1)
+plot(abs(x))
+axis([0,256,0,2])
+title('Original Signal')
+xlabel('t')
+ylabel('r(t)')
+
+
+subplot(4,1,2)
+plot(abs(y1))
+axis([0,256,0,1])
+title('HDS')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(4,1,3)
+plot(abs(ze(y2)))
+axis([0,256,0,1])
+title('S(MDS)')
+xlabel('t')
+ylabel('r(t)')
+
+subplot(4,1,4)
+plot(abs(ze(y3)))
+axis([0,256,0,1])
+title('S(LDS)')
+xlabel('t')
+ylabel('r(t)')
\ No newline at end of file

idear.m

+function [ya,yp]=idear(xa,xp,maxa)
+
+for n=1:length(xa)
+    if (xa(n)<maxa)||(xa(n)==maxa)
+        ya(n)=xa(n);
+        yp(n)=xp(n);
+    else
+        ya(n)=maxa;
+        yp(n)=xp(n);
+    end
+end
\ No newline at end of file

main.m

+clear
+clc
+% load precode.mat
+len=100000;
+% len=5000;
+% p1s=zeros(1,len);
+% p2s=zeros(1,len);
+% p3s=zeros(1,len);
+% p_scfde=zeros(1,len);
+% p_pre=zeros(1,len);
+% p=zeros(1,len);
+% p_Tofdm=zeros(1,len);
+
+
+for n=1:len
+    x=ofdm_signal();
+    xa=abs(x);
+%     x_pre=pre_ofdm_signal();
+%     [y1,y2,y3,p1s(n),p2s(n),p3s(n)]=average(x,0.5,1);%p1s(n)为y1的PAPR，以此类推
+     [yy1,yy2,yy3,pp1(n),pp2(n),pp3(n)]=cubic(x,1,1);
+%     p(n)=10*log10(max(xa).^2/(mean(xa.^2)));%原始OFDM信号x的PAPR
+%      p_pre(n)=10*log10(max(abs(x_pre)).^2/(mean(abs(x_pre).^2)));
+end
+
+% ya1=abs(y1);
+% ya2=abs(y2);
+% ya3=abs(y3);
+
+yya1=abs(yy1);
+yya2=abs(yy2);
+yya3=abs(yy3);
+
+% [papr,prob]=ccdf1(p,len);
+% [papr_pre,prob_pre]=ccdf1(p_pre,len);
+% for n=1:len
+% %     x_scfde=abs(scfde());
+% %     p_scfde(n)=10*log10(max(abs(x_scfde)).^2/(mean(abs(x_scfde).^2)));
+%     x_Tofdm=Tofdm_signal(256);
+%     xa_Tofdm=abs(x_Tofdm);
+%     p_Tofdm(n)=10*log10(max(xa_Tofdm).^2/(mean(xa_Tofdm.^2)));
+% end
+
+% [papr_scfde,prob_scfde]=ccdf1(p_scfde,len);
+% [papr_Tofdm,prob_Tofdm]=ccdf1(p_Tofdm,len);
+% 
+% [papr1,prob1]=ccdf1(p1s,len);
+% [papr2,prob2]=ccdf1(p2s,len);
+% [papr3,prob3]=ccdf1(p3s,len);
+% 
+% [papr_cubic1,prob_cubic1]=ccdf1(pp1,len);
+% [papr_cubic2,prob_cubic2]=ccdf1(pp2,len);
+% [papr_cubic3,prob_cubic3]=ccdf1(pp3,len);
+% 
+% figure(2)
+% semilogy(papr,prob,'-o')
+% axis([0,20,10^(-5),1])
+% hold on
+% semilogy(papr1,prob1,'-ro')
+% hold on
+% semilogy(papr2,prob2,'-go')
+% hold on
+% semilogy(papr3,prob3,'-yo')
+% hold on
+% semilogy(papr_pre,prob_pre,'-squre')
+% hold on
+% % semilogy(papr_cubic1,prob_cubic1,'-^')
+% semilogy(papr_cubic1,prob_cubic1,'-^','linewidth',1,'markersize',6)
+% hold on
+% % semilogy(papr_cubic2,prob_cubic2,'-r^')
+% semilogy(papr_cubic2,prob_cubic2,'-r^','linewidth',1,'markersize',6)
+% hold on
+% % semilogy(papr_cubic3,prob_cubic3,'-g^')
+% semilogy(papr_cubic3,prob_cubic3,'-g^','linewidth',1,'markersize',6)
+% hold on
+% semilogy(papr_Tofdm,prob_Tofdm,'-r')
+
+
+
+

matrix1.m

+clear
+clc
+
+N=8;
+s0=zeros(1,N);
+s0(1)=1;
+s0(N/2+1)=1;
+
+s1=zeros(1,N);
+s1(1)=1;
+s1(N/2+1)=-1;
+
+S0=[s0,s0];
+S1=[s1,s1];
+
+for n=1:N
+    m(n)=(S0(n-1+N)-S1(n-3+N));
+end
+m=m'/2;
+
+% M=fft(m)/sqrt(N);
+
+% M1=[1,-exp(-1i*2*pi*3/8),exp(-1i*2*pi*2/8),-exp(-1i*2*pi*9/8),exp(-1i*2*pi*4/8),-exp(-1i*2*pi*15/8),exp(-1i*2*pi*6/8),-exp(-1i*2*pi*21/8)]/sqrt(8);
+
+% m1=ifft(M1)/0.1768;
+
+for k=1:N
+    precode(:,k)=circshift(m,k-1);
+end
\ No newline at end of file

modify.m

+function j=modify(x)
+k=1;
+j=1;
+for n=k:length(x)-k
+    if (x(n)>0)||(x(n+k)>0)||(x(n-k)>0)
+        j=j+1;
+    end
+end

nc.m

+function y=nc(x)
+y=zeros(1,21);
+for n=1:21
+    
+y(n)=interp1(x(n,:),16:2:20,10E-5,'spline');
+end
\ No newline at end of file

ofdm_signal.m

+function y=ofdm_signal()
+% db=0;
+ibo=3;
+num=sqrt((10/1)*10^(ibo/10));
+
+
+mo=modem.qammod(16); %调制解调类包中的QAM调制器
+data1=randi(16,256,1)-1; %生成随机传输符号
+x1=modulate(mo,data1)/num;%用QAM调制器对信号进行调制
+% x1=(x2)*(10^(db/20));
+
+y1=ifft(x1)*sqrt(length(x1));
+
+% num=max(abs(y1))/sqrt(3);
+
+y=y1;
\ No newline at end of file

ofdmpe.m

+x=scfde();
+
+for n=1:length(x)
+    if x(n)>sqrt(3)
+        x(n)=sqrt(3);
+    end
+end
+
+pe=mean(abs(x).^2)/3;
\ No newline at end of file

p2c.m

+function y = p2c( abs, angle )
+    y=abs.*cos(angle)+abs.*sin(angle)*i;
+end
+

panduan.m

+function y=panduan(x)
+a=zeros(1,length(x));
+y=zeros(1,length(x));
+N=length(x);
+
+for n=1:N
+    if x(n)>0
+        a(n)=1;
+    else
+        a(n)=0;
+    end
+end
+
+for n=1:N
+    if n==1
+        y(n)=or(or(a(1),a(1)),a(2));
+    elseif n==N
+        y(n)=or(or(a(N),a(N-1)),a(N-1));
+    else
+        y(n)=or(or(a(n),a(n+1)),a(n-1));
+    end
+end
+

paprdecrease.m

+len=50000;
+
+for n=1:len
+    x=ofdm_signal();
+    xa=abs(x);
+    p(n)=10*log10(max(xa).^2/(mean(xa.^2)));%%PAPRļ
+end
+
+pow=mean(xa.^2);
+
+[papr,prob]=ccdf1(p,len);
+f1=(1-exp(-1/pow));
+f2=(exp(-1/pow)-exp(-2/pow));
+f3=exp(-2/pow);
+p1=pow-(1+pow)*exp(-1/pow);
+p2=(1+pow)*exp(-1/pow)-(2+pow)*exp(-2/pow);
+p3=(2+pow)*exp(-2/pow);
+
+papr1=10*log10(3*(p1+p2/2+p3/3)/pow);
+papr2=10*log10(3*(p2/2+p3/3)/((f2+f3)*pow));
+papr3=10*log10(p3/(f3*pow));
+%prob1= [1 1	1 1	1 1	1 1 1 1 1 0.9809 0.95742 0.9325 0.86172 0.65288 0.41188 0.22256 0.10536 0.04552 0.01704 0.006 0.00186 0.00062 0.00018 2.0e-05 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
+%ofdmpr=[-3.84355474610437,-3.34355474610437,-2.84355474610437,-2.34355474610437,-1.84355474610437,-1.34355474610437,-0.843554746104370,-0.343554746104370,0.156445253895630,0.656445253895630,1.15644525389563,1.65644525389563,2.15644525389563,2.65644525389563,3.95644525389563,4.95644525389563,5.65644525389563,6.25644525389563,7.65644525389563,8.35644525389563,9.95644525389563,10.75644525389563,11.75644525389563,12.95644525389563,13.15644525389563,13.65644525389563,14.15644525389563,15.65644525389563,16.1564452538956,16.6564452538956,17.1564452538956,17.6564452538956,18.1564452538956,18.6564452538956,19.1564452538956,19.6564452538956,20.1564452538956,20.6564452538956,21.1564452538956,21.6564452538956,22.1564452538956];
+
+figure(1)
+hold on
+%semilogy(papr-papr1-0.1+3,prob,'o')
+semilogy(papr-papr1,prob,'o')
+hold on
+%semilogy(papr-papr2+3,prob,'*')
+semilogy(papr-papr2,prob,'*')
+hold on
+semilogy(papr-papr3,prob,'s')
+hold on 
+%semilogy(ofdmpr+3,prob1,'d')
+semilogy(papr,prob,'d')
+
+% hold on 
+% semilogy(papr-papr1,prob1,'x')
+axis([0,20,10^(-4),1]);
+
+% figure(2)
+% hold on 
+% plot(papr-papr1-0.1,prob,'o')
+% hold on
+% plot(papr-papr2-0.1,prob,'*')
\ No newline at end of file

peak.m

+function y=peak()
+    ll=150;
+    b=[];
+    for m=1:ll
+        a=reshape(dbplot(),41,1);
+        b(m,:)=a;
+    end
+    y=b;
+end
\ No newline at end of file

plotSignal.m

+subplot(4,1,1)
+plot(xa,'DisplayName','xa','YDataSource','xa');
+xlabel('t');
+ylabel('r(t)');
+title('Original Signal');
+axis([0 256 0 2])
+subplot(4,1,2)
+plot(yya1,'DisplayName','yya1','YDataSource','yya1');
+xlabel('t');
+ylabel('r(t)');
+title('HDS');
+axis([0 256 0 1])
+subplot(4,1,3)
+plot(yya2,'DisplayName','yya2','YDataSource','yya2');
+xlabel('t');
+ylabel('r(t)');
+title('MDS');
+axis([0 256 0 1])
+subplot(4,1,4)
+plot(yya3,'DisplayName','yya3','YDataSource','yya3');
+xlabel('t');
+ylabel('r(t)');
+title('LDS');
+axis([0 256 0 1])
+

plotfigure.m

+plot(snr1(1:5),p1(1:5),'-o')
+hold on
+plot(snr2(1:5),p2(1:5),'-ro')
+hold on
+% plot(snr3(1:6),p3(1:6),'-go')
+hold on
+% plot(snr4(1:5),p4u(1:5),'-squre')
\ No newline at end of file

powerefficiency.m

+num=1;
+for loop=1:5000
+    load precode.mat
+    len=1;
+    p1s=zeros(1,len);
+    p2s=zeros(1,len);
+    p3s=zeros(1,len);
+    p_scfde=zeros(1,len);
+    p_pre=zeros(1,len);
+    p=zeros(1,len);
+    
+    
+    for n=1:len
+        x=ofdm_signal();
+        xa=abs(x);
+        [y1,y2,y3,p1s(n),p2s(n),p3s(n)]=average(x,1,1);
+        %     [yy1,yy2,yy3,pp1(n),pp2(n),pp3(n)]=cubic(x,1,1);
+        p(n)=10*log10(max(xa).^2/(mean(xa.^2)));
+    end
+    pow=mean(xa.^2);
+    pe_ofdm(num)=pow/3; 
+    f1=(1-exp(-1/pow)); 
+    f2=(exp(-1/pow)-exp(-2/pow));  
+    f3=exp(-2/pow);%-exp(-3/pow); 
+    
+    p1=pow-(1+pow)*exp(-1/pow); %pow1
+    p2=(1+pow)*exp(-1/pow)-(2+pow)*exp(-2/pow);%pow2
+    p3=(2+pow)*exp(-2/pow);%-(3+pow)*exp(-3/pow);%pow3
+    
+    
+    % PE T%
+%     t_p1(num)=(p1+p2/2+p3/3)/(f1+f2+f3);
+%     t_p2(num)=(p2/2+p3/3)/(f2+f3);
+%     t_p3(num)=(p3/3)/(f3);
+    
+    t_p1(num)=10*log10(3*(p1+p2/2+p3/3)/(p1+p2+p3));
+    t_p2(num)=10*log10(3*(p2/2+p3/3)/(f2+f3)*(p1+p2+p3));
+    t_p3(num)=10*log10(p3/(f3)*(p1+p2+p3));
+    
+    %
+    tt_p1=(p1+p2/2+p3/3);
+    tt_p2=(p2/2+p3/3);
+    tt_p3=(p3/3);
+    
+    
+%     temp1(num)=mean(abs(y1).^2);
+%     y_new2=ze(y2);
+%     y_new3=ze(y3);
+    
+%     n1=length(y_new2);
+%     zeros1=zeros(1,modify(y2));
+%     y_new2_new=cat(2,y_new2,zeros1);
+%     
+%     n2=length(y_new2);
+%     zeros1=zeros(1,ceil(n1/2));
+%     y_new2_new=cat(2,y_new2,zeros1);
+    % for k=1:length(y_new2)
+    %     if abs(y_new2(k)>1)
+    %         y_new2(k)=1;
+    %     end
+    % end
+%     temp2(num)=mean(abs(y_new2).^2);
+    
+%     temp2_new(num)=mean(abs(y_new2_new).^2);
+    % for k=1:length(y_new3)
+    %     if abs(y_new3(k)>1)
+    %         y_new3(k)=1;
+    %     end
+    % end
+%     temp3(num)=mean(abs(y_new3).^2);
+    num=num+1;
+end
+pe1=mean(t_p1);
+pe2=mean(t_p2);
+pe3=mean(t_p3);
+
+% pe1=mean(tt_p1);
+% pe2=mean(tt_p2);
+% pe3=mean(tt_p3);
+
+% spe1=mean(temp1);
+% spe2=mean(temp2);
+% spe3=mean(temp3);
+% 
+% sspe2=mean(temp2_new);
\ No newline at end of file

pprob3.m

+p1s=zeros(1,len);
+p2s=zeros(1,len);
+p3s=zeros(1,len);
+for n=1:len
+    x=ofdm_signal();
+    [y1,y2,y3,p1s(n),p2s(n),p3s(n)]=average(x,1,1);
+end
+
+
+[papr3,prob3]=ccdf1(p3s,len);
+
+
+
+

pre_ofdm_signal.m

+function y=pre_ofdm_signal()
+% db=0;
+% ibo=8;
+num=sqrt(10);%sqrt((10/1)*10^(ibo/10));
+load precode.mat
+
+mo=modem.qammod(16);
+data1=randi(15,256,1); %
+x1=modulate(mo,data1)/num;
+% x1=(x2)*(10^(db/20));
+y1=ifft(precode*x1)*sqrt(length(x1));
+
+% num=max(abs(y1))/sqrt(3);
+
+y=y1;

prob.m

+function y=prob(x1,x2,t)
+
+%y=-t*x2*exp(-(t*x2)^2)+t*x1*exp(-(t*x1)^2)+sqrt(pi)/2*(erf(t*x2)-erf(t*x1));
+
+y=t*x1*exp(-(t*x1)^2)-t*x2*exp(-(t*x2)^2)+(sqrt(pi)/2)*(erf(t*x2)-erf(t*x1));
\ No newline at end of file

prob3.asv

+function prob3
+p1s=zeros(1,len);
+p2s=zeros(1,len);
+p3s=zeros(1,len);
+p_scfde=zeros(1,len);
+p_pre=zeros(1,len);
+p=zeros(1,len);
+p_Tofdm=zeros(1,len);
+
+
+for n=1:len
+    x=ofdm_signal();
+    xa=abs(x);
+    x_pre=pre_ofdm_signal();
+    [y1,y2,y3,p1s(n),p2s(n),p3s(n)]=average(x,1,1);%p1s(n)Ϊy1PAPRԴ
+end
+
+ya1=abs(y1);
+ya2=abs(y2);
+ya3=abs(y3);
+[papr,prob]=ccdf1(p,len);
+[papr_pre,prob_pre]=ccdf1(p_pre,len);
+[papr3,prob3]=ccdf1(p3s,len);
+
+
+end
+

process.m

+function [BER,v1,v2,v3,p]=process(ibo,snr,method)
+% tic
+% clear
+N=4096;
+% snr=20;
+% ibo=-3;
+% method='average';
+dd=1;
+num=sqrt((10/dd)*10^(ibo/10));
+
+loop=100;
+ber=zeros(1,loop);
+% e=zeros(1,loop);
+me=modem.qammod(16);
+%     xT=addzeros(xo,dd);
+de=modem.qamdemod(16);
+
+
+%------------------------------------------------
+for start=1:loop
+%     data=r(xT)*sqrt(N);
+    data=randi(16,N,1)-1;
+    xo=modulate(me,data)/num;
+    d=ifft(xo);
+    yOFDM=d*sqrt(N);
+    Ina=abs(yOFDM);
+    %------------------------------------------
+    switch method
+        case 'average'
+            [y1,y2,y3,p1,p2,p3]=average(yOFDM,1,1);
+        case 'cubic'
+            [y1,y2,y3,p1,p2,p3]=cubic(yOFDM,1,1);
+        case 'cubic_delta'
+            [y1,y2,y3,p1,p2,p3]=cubic_delta(yOFDM);
+    end
+    %------------------------------------------
+    %--------------------the first layer
+    Ina_a=abs(y1);
+    Inp_a=angle(y1);
+ 
+    [ya_a,yp_a]=idear(Ina_a,Inp_a,sqrt(2));
+    y_a=p2c(ya_a,yp_a);
+    
+    %-----------------------the second layer
+    Ina_b=abs(y2);
+    Inp_b=angle(y2);
+    
+    [ya_b,yp_b]=idear(Ina_b,Inp_b,sqrt(2));
+    y_b=p2c(ya_b,(yp_b));
+    %-----------------------the third layer
+    Ina_c=abs(y3);
+    Inp_c=angle(y3);
+    
+    [ya_c,yp_c]=idear(Ina_c,Inp_c,sqrt(2));
+    y_c=p2c(ya_c,yp_c);
+    %-----------------------------------------------ʼ
+    y=y_a+y_b+y_c;
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+    
+    Go=Yreceive;     %
+    
+    Gone=fft(Go)/sqrt(N);  %fft
+    Gone1=reshape(Gone,1,length(Gone));  %任
+    % receive=Gone*num;
+    
+%     Gone1=quitzero(Gone,N);
+    Result=demodulate(de,Gone1*num);
+%     z=Gone1*num;
+%     e(start)=mean(abs(Gone1-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+%     [bes,ber(start)]=biterr(Result,data);
+    [bes,ber(start)]=biterr(Result,data');
+%     [ses,ser(start)]=symerr(Result,data);
+end
+
+BER=mean(ber);
+xiaohao_a=mean(panduan(abs(y_a)),2);
+xiaohao_b=mean(panduan(abs(y_b)),2);
+xiaohao_c=mean(panduan(abs(y_c)),2);
+shuchu_a=mean(abs(y_a));
+shuchu_b=mean(abs(y_b));
+shuchu_c=mean(abs(y_c));
+v1=shuchu_a/xiaohao_a;
+v2=shuchu_b/xiaohao_b;
+v3=shuchu_c/xiaohao_c;
+p=(v1+v2+v3)/3;
+
+
+% gp2=ze(y_b);
+% gp3=ze(y_c);
+% v1=mean(abs(y_a)/1);
+% v2=mean(abs(gp2)/1);
+% v3=mean(abs(gp3)/1);
+% p=(v1+v2+v3)/3;
\ No newline at end of file

process_ber.m

+function BER=process_ber(ibo,snr)
+% ibo=6;
+% snr=20;
+N=4096;
+dd=1;
+num=sqrt((10/dd)*10^(ibo/10));
+loop=100;
+ber=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+%------------------------------------------------
+for start=1:loop
+    data=randi(16,N,1)-1;
+    xo=modulate(me,data)/num;
+    d=ifft(xo);
+    yOFDM=d*sqrt(N);
+    
+    Ina_a=abs(yOFDM);
+    Inp_a=angle(yOFDM);
+ 
+    [ya,yp]=idear(Ina_a,Inp_a,sqrt(1.25));
+    yOFDM_lim=p2c(ya,yp);
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=yOFDM_lim;
+    else
+        Yreceive=awgn(yOFDM_lim,snr,'measured');
+    end 
+    Go=Yreceive;     %
+    Gone=fft(Go)/sqrt(N);  %fft
+    Gone1=reshape(Gone,1,length(Gone));  %任
+    Result=demodulate(de,Gone1*num);
+    [bes,ber(start)]=biterr(Result,data');
+end
+BER=mean(ber);
+

process_normal.m

+% function [BER,p1,pg2,pg3,p]=process(ibo,snr,method)
+% tic
+% clear
+snr='none';
+ibo=0;
+method='cubic';
+dd=1;
+% nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+
+loop=10000;
+ber=zeros(1,loop);
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+cap=N;
+num=sqrt((10/dd)*10^(ibo/10));
+aa=0.1;
+bb=0.1;
+
+Tab_a_a=ones(1,cap); %时域预失真幅度查询表
+Tab_p_a=zeros(1,cap); %时域预失真相位查询表
+address_a=zeros(1,N); %索引表
+
+Tab_a_b=ones(1,cap); %时域预失真幅度查询表
+Tab_p_b=zeros(1,cap); %时域预失真相位查询表
+address_b=zeros(1,N); %索引表
+
+Tab_a_c=ones(1,cap); %时域预失真幅度查询表
+Tab_p_c=zeros(1,cap); %时域预失真相位查询表
+address_c=zeros(1,N); %索引表
+%--------------------------------------------以上为策略定义,之后为增加速度定义
+Ina1_a=ones(1,N);
+Outa_a=ones(1,N);
+Outp_a=zeros(1,N);
+Ina1_b=ones(1,N);
+Outa_b=ones(1,N);
+Outp_b=zeros(1,N);
+Ina1_c=ones(1,N);
+Outa_c=ones(1,N);
+Outp_c=zeros(1,N);
+%------------------------------------------------
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=addzeros(xo,dd);
+    d=ifft(xT)*sqrt(N);
+    %------------------------------------------
+    switch method
+        case 'average'
+            [y1,y2,y3,p1,p2,p3]=average(d);
+        case 'cubic'
+            [y1,y2,y3,p1,p2,p3]=cubic(d);
+    end
+    %------------------------------------------
+    %--------------------the first layer
+    Ina_a=abs(y1);
+    Inp_a=angle(y1);
+    for uj=1:length(Ina_a)                          %时域预失真
+        [address_a(uj),Ina1_a(uj)]=Q(Ina_a(uj),1,cap);% Q: adrress function
+        Outa_a(uj)=Ina1_a(uj)*Tab_a_a(address_a(uj));
+        Outp_a(uj)=Inp_a(uj)+Tab_p_a(address_a(uj));
+    end
+    
+    [ya_a,yp_a]=uHPA(Outa_a,Outp_a,1);
+    y_a=p2c(ya_a,(yp_a));
+    
+    for ui=1:length(address_a)    %时域预失真更新公式
+        Tab_a_a(address_a(ui))=Tab_a_a(address_a(ui))-aa*sigmoid(ya_a(ui)-Ina1_a(ui));
+        Tab_p_a(address_a(ui))=Tab_p_a(address_a(ui))-bb*sigmoid(Ep(yp_a(ui)-Inp_a(ui)));
+    end
+    %-----------------------the second layer
+    Ina_b=abs(y2);
+    Inp_b=angle(y2);
+    for uj=1:length(Ina_b)                          %时域预失真
+        [address_b(uj),Ina1_b(uj)]=Q(Ina_b(uj),1,cap);% Q: adrress function
+        Outa_b(uj)=Ina1_b(uj)*Tab_a_b(address_b(uj));
+        Outp_b(uj)=Inp_b(uj)+Tab_p_b(address_b(uj));
+    end
+    
+    [ya_b,yp_b]=uHPA(Outa_b,Outp_b,1);
+    y_b=p2c(ya_b,(yp_b));
+    
+    for ui=1:length(address_b)    %时域预失真更新公式
+        Tab_a_b(address_b(ui))=Tab_a_b(address_b(ui))-aa*sigmoid((ya_b(ui)-Ina1_b(ui)));
+        Tab_p_b(address_b(ui))=Tab_p_b(address_b(ui))-bb*sigmoid(Ep(yp_b(ui)-Inp_b(ui)));
+    end
+    %-----------------------the third layer
+    Ina_c=abs(y3);
+    Inp_c=angle(y3);
+    for uj=1:length(Ina_c)                          %时域预失真
+        [address_c(uj),Ina1_c(uj)]=Q(Ina_c(uj),1,cap);% Q: adrress function
+        Outa_c(uj)=Ina1_c(uj)*Tab_a_c(address_c(uj));
+        Outp_c(uj)=Inp_c(uj)+Tab_p_c(address_c(uj));
+    end
+    
+    [ya_c,yp_c]=uHPA(Outa_c,(Outp_c),1);
+    y_c=p2c(ya_c,yp_c);
+    
+    for ui=1:length(address_c)    %时域预失真更新公式
+        Tab_a_c(address_c(ui))=Tab_a_c(address_c(ui))-aa*sigmoid(ya_c(ui)-Ina1_c(ui));
+        Tab_p_c(address_c(ui))=Tab_p_c(address_c(ui))-bb*sigmoid(Ep(yp_c(ui)-Inp_c(ui)));
+    end
+    %-----------------------------------------------开始做反馈
+    y=y_a+y_b+y_c;
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+    
+    Go=Yreceive;     %反馈
+    
+    Gone=fft(Go)/sqrt(N);  %反馈做fft
+    Gone=reshape(Gone,1,length(Gone));  %串并变换
+    % receive=Gone*num;
+    
+    Gone1=quitzero(Gone,N);
+    Result=demodulate(de,Gone1*num);
+%     z=Gone1*num;
+    e(start)=mean(abs(Gone1-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+%     [ses,ser(start)]=symerr(Result,data);
+end
+BER=mean(ber(loop-2500:loop));
+% SER=mean(ser(loop-2500:loop));
+
+p1=mean(spower(y_a));
+% p2=mean(spower(y_b));
+% p3=mean(spower(y_c));
+% 
+% g2=length(ze(y_b))/length(y_b);
+% g3=length(ze(y_c))/length(y_c);
+% 
+gp2=ze(y_b);
+gp3=ze(y_c);
+
+pg2=mean(spower(gp2));
+pg3=mean(spower(gp3));
+% 
+% p=mean(spower(y)/3);
+% p_n=(p1+pg2+pg3)/3;
+p=mean(spower(y)/3);
+plot(Gone1,'*')
+% figure(2)
+% plot(z,'*')
+% toc
\ No newline at end of file

s.m

+function y=s()
+% db=0;
+% ibo=8;
+num=sqrt(10);%sqrt((10/1)*10^(ibo/10));
+
+
+mo=modem.qammod(16);
+data1=randint(256,1,16); %
+x1=modulate(mo,data1)/num;
+% x1=(x2)*(10^(db/20));
+
+y1=ifft(x1)*sqrt(length(x1));
+num=1;
+% num=max(abs(y1))/sqrt(3);
+
+y=y1/num;
\ No newline at end of file

scfde.m

+function y=scfde()
+% db=0;
+% ibo=8;
+num=sqrt(10);%sqrt((10/1)*10^(ibo/10));
+
+
+mo=modem.qammod(16);
+data1=randi(15,256,1); %
+x1=modulate(mo,data1)/num;
+y=x1;
\ No newline at end of file

scfdev1.m

+x_scfde=abs(scfde());
+len=length(x_scfde);
+
+for n=1:len
+    p_scfde(n)=10*log10(x_scfde(n)).^2/(mean(abs(x_scfde).^2));
+end
+
+[papr_scfde,prob_scfde]=ccdf1(p_scfde,len);
+
+figure(2)
+semilogy(papr_scfde,prob_scfde,'-squre')
\ No newline at end of file

segment.m

+function [y1,y2,y3,p1,p2,p3]=segment(x)
+
+[a,p]=c2p(x);
+
+for n=1:length(x)
+    if a(n)<0.5
+        y1(n)=p2c(a(n),p(n));
+        y2(n)=0;
+        y3(n)=0;
+    elseif ((a(n)>0.5)&&(a(n)<1.5))
+        y1(n)=p2c(0.5,p(n));
+        y2(n)=p2c(a(n)-0.5,p(n));   
+        y3(n)=0;
+    else
+        y1(n)=p2c(0.5,p(n));   
+        y2(n)=p2c(1,p(n));
+        y3(n)=p2c(a(n)-1.5,p(n));  
+    end
+end
+
+p1=funcpapr(y1);
+p2=funcpapr(y2);
+p3=funcpapr(y3);

segmentmain.m

+function [BER,v1,v2,v3,p,p2]=segmentmain(ibo,snr)
+% tic
+
+dd=1;
+% nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+% ibo=0;
+% snr='none';
+
+loop=1000;
+ber=zeros(1,loop);
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+num=sqrt((10/dd)*10^(ibo/10));
+%------------------------------------------------
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=addzeros(xo,dd);
+    d=ifft(xT)*sqrt(N);
+    %------------------------------------------
+    [y1,y2,y3,p1,p2,p3]=segment(d);
+    %------------------------------------------
+    %--------------------the first layer
+    Ina_a=abs(y1);
+    Inp_a=angle(y1);
+
+    [ya_a,yp_a]=idear(Ina_a,Inp_a,0.5);
+    y_a=p2c(ya_a,(yp_a));
+
+    %-----------------------the second layer
+    Ina_b=abs(y2);
+    Inp_b=angle(y2);
+
+    [ya_b,yp_b]=idear(Ina_b,Inp_b,1);
+    y_b=p2c(ya_b,(yp_b));
+    %-----------------------the third layer
+    Ina_c=abs(y3);
+    Inp_c=angle(y3);
+
+    [ya_c,yp_c]=idear(Ina_c,Inp_c,sqrt(1.75));
+    y_c=p2c(ya_c,yp_c);
+    %-----------------------------------------------ʼ
+    y=(y_a+y_b+y_c);
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+
+    Go=Yreceive;     %
+
+    Gone=fft(Go)/sqrt(N);  %fft
+    Gone=reshape(Gone,1,length(Gone));  %任
+    % receive=Gone*num;
+
+    Gone1=quitzero(Gone,N);
+    Result=demodulate(de,Gone1*num);
+    z=Gone1*num;
+    e(start)=mean(abs(Gone1-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+    %     [ses,ser(start)]=symerr(Result,data);
+end
+BER=mean(ber);
+
+ypa=mean(abs(y_a)/0.5);
+ypb=mean(ze(abs(y_b))/1);
+ypc=mean(ze(abs(y_c))/sqrt(1.75));
+p2=(ypa+ypb+ypc)/3;
+
+v1=mean(abs(y_a)/0.5);
+v2=mean(abs(y_b)/1);
+v3=mean(abs(y_c)/sqrt(1.75));
+p=(v1+v2+v3)/3;
\ No newline at end of file

semilogyfigure.m

+clear all
+
+figure(1)
+obj = get(gca,'children');
+x1=get(obj(1),'xdata');
+y1=get(obj(1),'ydata');
+
+x10=get(obj(10),'xdata');
+y10=get(obj(10),'ydata');
+
+x2=get(obj(2),'xdata');
+y2=get(obj(2),'ydata');
+
+x3=get(obj(3),'xdata');
+y3=get(obj(3),'ydata');
+
+x4=get(obj(4),'xdata');
+y4=get(obj(4),'ydata');
+
+x5=get(obj(5),'xdata');
+y5=get(obj(5),'ydata');
+
+x6=get(obj(6),'xdata');
+y6=get(obj(6),'ydata');
+
+x7=get(obj(7),'xdata');
+y7=get(obj(7),'ydata');
+
+x8=get(obj(8),'xdata');
+y8=get(obj(8),'ydata');
+
+x9=get(obj(9),'xdata');
+y9=get(obj(9),'ydata');
+
+figure(2)
+h1=semilogy(x1,y1);
+hold on
+h2=semilogy(x2,y2);
+hold on
+h3=semilogy(x3,y3);
+hold on
+h4=semilogy(x4,y4);
+hold on
+h5=semilogy(x5,y5);
+hold on
+h6=semilogy(x6,y6);
+hold on
+h7=semilogy(x7,y7);
+hold on
+h8=semilogy(x8,y8);
+hold on
+h9=semilogy(x9,y9);
+hold on
+h10=semilogy(x10,y10);
+hold on
\ No newline at end of file

stemplot.m

+% function  stemplot()
+a=peak();
+b=sum(a);
+figure
+bar(b,'grouped');
+axis([20 61 0 8000])
+set(gca,'XTickLabel',4:12);
+set(gca,'YTickLabel',0:2000:8000);
+% end
+

tdpd.m

+% function [BER,v]=tdpd(ibo,snr)
+snr='none';
+ibo=0;
+dd=1;
+% nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+
+loop=1;
+ber=zeros(1,loop);
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+
+num=sqrt((10/dd)*10^(ibo/10));
+
+
+
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=xo;
+    d=ifft(xT)*sqrt(N);
+    Ina=abs(d);
+    Inp=angle(d);
+    
+   
+    %-----------------------------------------------ʼ
+    
+    %------------------------------------------------------------
+    
+%     [ses,ser(start)]=symerr(Result,data);
+end
+
+[papr0,pp]=ccdf(Ina);
+plot(papr0,pp)

tdpdproess.m

+% function [BER,p]=tdpdproess(ibo,snr)
+snr='none';
+ibo=3.23;
+dd=1;
+% nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+method='average';
+loop=10000;
+ber=zeros(1,loop);
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+cap=N;
+num=sqrt((10/dd)*10^(ibo/10));
+aa=0.1;
+bb=0.1;
+
+Tab_a=ones(1,cap); %时域预失真幅度查询表
+Tab_p=zeros(1,cap); %时域预失真相位查询表
+address=zeros(1,N); %索引表
+%--------------------------------------------以上为策略定义,之后为增加速度定义
+% Ina=ones(1,N);
+Ina1=ones(1,N);
+Outa=ones(1,N);
+Outp=zeros(1,N);
+
+
+
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=xo;
+    d=ifft(xT)*sqrt(N);
+    Ina=abs(d);
+    Inp=angle(d);
+    for uj=1:length(Ina)                          %时域预失真
+        [address(uj),Ina1(uj)]=Q(Ina(uj),sqrt(3),cap);% Q: adrress function
+        Outa(uj)=Ina1(uj)*Tab_a(address(uj));
+        Outp(uj)=Inp(uj)+Tab_p(address(uj));
+    end
+    
+    d1=p2c(Ina1,Inp);
+    switch method
+        case 'average'
+            [y1,y2,y3,p1,p2,p3]=average(d1);
+        case 'cubic'
+            [y1,y2,y3,p1,p2,p3]=cubic(d1);
+    end
+    %------------------------------------------
+    %--------------------the first layer
+    Ina_a=abs(y1);
+    Inp_a=angle(y1);
+    
+    [ya_a,yp_a]=HPA(Ina_a,Inp_a);
+    y_a=p2c(ya_a,(yp_a));
+    
+    %-----------------------the second layer
+    Ina_b=abs(y2);
+    Inp_b=angle(y2);
+    
+    [ya_b,yp_b]=HPA(Ina_b,Inp_b);
+    y_b=p2c(ya_b,(yp_b));
+    
+    %-----------------------the third layer
+    Ina_c=abs(y3);
+    Inp_c=angle(y3);
+    
+    [ya_c,yp_c]=HPA(Ina_c,Inp_c);
+    y_c=p2c(ya_c,yp_c);
+    %-------------------------------
+    y=y_a+y_b+y_c;
+    %-----------------------------------------------开始做反馈
+    ya=abs(y);
+    yp=angle(y);
+    
+    for ui=1:length(address)    %时域预失真更新公式
+        Tab_a(address(ui))=Tab_a(address(ui))-aa*(ya(ui)-Ina1(ui));
+        Tab_p(address(ui))=Tab_p(address(ui))-bb*Ep(yp(ui)-Inp(ui));
+    end
+    %------------------------------------------------------------
+    if strcmp(snr,'none')
+        Yreceive=y;
+    else
+        Yreceive=awgn(y,snr,'measured');
+    end
+    
+    Go=Yreceive;     %反馈
+    
+    Gone=fft(Go)/sqrt(N);  %反馈做fft
+    Gone=reshape(Gone,1,length(Gone));  %串并变换
+    Result=demodulate(de,Gone*num);
+    z=Gone*num;
+    e(start)=mean(abs(Gone-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+    %     [ses,ser(start)]=symerr(Result,data);
+end
+BER=mean(ber(loop-2500:loop));
+% SER=mean(ser(loop-2500:loop));
+p=mean(spower(y)/3);
+plot(z,'*')

tested.m

+function [BER]=tested(ibo)
+% snr='none';
+% ibo=0;
+dd=1;
+nc=(randn(1,256)+randn(1,256)*1i)/sqrt(2)/10;
+
+loop=2000;
+e=zeros(1,loop);
+me=modem.qammod(16);
+de=modem.qamdemod(16);
+N=256;
+cap=N;
+num=sqrt((10/dd)*10^(ibo/10));
+aa=0;
+bb=0;
+
+Tab_a=ones(1,cap); %时域预失真幅度查询表
+Tab_p=zeros(1,cap); %时域预失真相位查询表
+address=zeros(1,N); %索引表
+%--------------------------------------------以上为策略定义,之后为增加速度定义
+Ina=ones(1,N);
+Ina1=ones(1,N);
+Outa=ones(1,N);
+Outp=zeros(1,N);
+
+
+
+for start=1:loop
+    data=randint(1,N,16);
+    data1=modulate(me,data);
+    xo=data1/num;
+    xT=xo;
+    d=ifft(xT)*sqrt(N);
+    Ina=abs(d);
+    Inp=angle(d);
+    for uj=1:length(Ina)                          %时域预失真
+        [address(uj),Ina1(uj)]=Q(Ina(uj),sqrt(1),cap);% Q: adrress function
+    end
+    
+    [ya,yp]=HPA(Ina1,Inp);
+    y=p2c(ya,yp);
+    
+    
+    %-----------------------------------------------开始做反馈
+    
+    %------------------------------------------------------------
+%     if strcmp(snr,'none')
+%         Yreceive=y;
+%     else
+%         Yreceive=awgn(y,snr,'measured');
+        Yreceive=y+nc;
+%     end
+    
+    Go=Yreceive;     %反馈
+    
+    Gone=fft(Go)/sqrt(N);  %反馈做fft
+    Gone=reshape(Gone,1,length(Gone));  %串并变换
+    Result=demodulate(de,Gone*num);
+    z=Gone*num;
+    e(start)=mean(abs(Gone-xo).^2);
+    %     e(start)=mean(abs(z-data1(1:N/dd)).^2);
+    [bes,ber(start)]=biterr(Result,data);
+    [ses,ser(start)]=symerr(Result,data);
+end
+% BER=mean(ber(loop-2500:loop));
+BER=mean(ber);
+e0nc=10*log10(spower(y)/spower(nc))
+% SER=mean(ser(loop-2500:loop));
+% p=mean(spower(y)/3);
+% plot(z,'*')

testfunction.m

+x_ofdm=linspace(min(index_ofdm),max(index_ofdm));
+y_ofdm=interp1(index_ofdm,ber_ofdm,x_ofdm,'cubic');
+h1=semilogy(x_ofdm,y_ofdm,'-*','linewidth',1,'markersize',6);
+axis([0,30,10^(-8),1]);
+hold on;
+% semilogy(x_ofdm,y_ofdm,'*');
+% hold on;
+% scatter(x_ofdm,y_ofdm);
+% hold on;
+
+x_cubic=linspace(min(index_cubic),max(index_cubic));
+y_cubic=interp1(index_cubic,ber_cubic,x_cubic,'cubic');
+h2=semilogy(x_cubic,y_cubic,'-o','linewidth',1,'markersize',6);
+hold on;
+
+x_average=linspace(min(index_average),max(index_average));
+y_average=interp1(index_average,ber_average,x_average,'cubic');
+h3=semilogy(x_average,y_average,'-v','linewidth',1,'markersize',6);
+grid on;
+
+legend([h1,h2,h3],'OFDM','OSID','SUSD');
+xlabel('EbN0[dB]');
+ylabel('BER');
\ No newline at end of file

testsemilogy.m

+figure(1)
+h1=semilogy(papr,prob,'-o','linewidth',1,'markersize',6);
+axis([0,20,10^(-5),1])
+hold on
+h2=semilogy(papr1,prob1,'linewidth',1,'markersize',6);
+hold on
+semilogy(papr2,prob2,'g','linewidth',1,'markersize',6)
+hold on
+semilogy(papr3+0.5,prob3,'r','linewidth',1,'markersize',6)
+hold on
+h3=semilogy(papr_pre-1.2,prob_pre,'-squre','linewidth',1,'markersize',6);
+%hold on
+%semilogy(papr_scfde,prob_scfde,'-x','linewidth',2,'markersize',6)
+hold on
+h4=semilogy(papr_Tofdm-1,prob_Tofdm,'-gv','linewidth',1,'markersize',6);
+hold on 
+% 理论值：利用原始PAPR减去降低的dB数即为理论值
+h5=scatter(papr-3.8,prob,80,'o') ;
+hold on 
+scatter(papr-6.4,prob,80,'o') 
+hold on 
+scatter(papr-7.4,prob,80,'o') 
+%hold on
+%semilogy(papr_cubic1,prob_cubic1,'-^')
+%hold on
+%semilogy(papr_cubic2,prob_cubic2,'-r^')
+%hold on
+%semilogy(papr_cubic3,prob_cubic3,'-g^')
+legend([h1,h3,h4,h5,h2],'OFDM','Precode','T-OFDM','Theroy','Simulation')
+xlabel('D(dB)');
+ylabel('Prob(PAPR>D)');
\ No newline at end of file

testsemilogy1.m

+figure(1)
+h1=semilogy(papr,prob,'-o','linewidth',2,'markersize',6);
+axis([0,14,10^(-5),1])
+hold on
+h2=semilogy(papr1,prob1,'linewidth',2,'markersize',6);
+hold on
+semilogy(papr2,prob2,'g','linewidth',2,'markersize',6)
+hold on
+semilogy(papr3+0.5,prob3,'r','linewidth',2,'markersize',6)
+hold on
+h3=semilogy(papr_pre-1.2,prob_pre,'-squre','linewidth',2,'markersize',6);
+hold on
+h6=semilogy(papr_cubic1,prob_cubic1,'-^','linewidth',2,'markersize',6);
+hold on
+h7=semilogy(papr_cubic2,prob_cubic2,'-g^','linewidth',2,'markersize',6);
+hold on
+h8=semilogy(papr_cubic3+2,prob_cubic3,'-r^','linewidth',2,'markersize',6);
+%hold on
+%semilogy(papr_scfde,prob_scfde,'-x','linewidth',2,'markersize',6)
+hold on
+h4=semilogy(papr_Tofdm-1,prob_Tofdm,'-gv','linewidth',2,'markersize',6);
+hold on 
+% 理论值：利用原始PAPR减去降低的dB数即为理论值
+h5=scatter(papr-3.8,prob,80,'o') ;
+hold on 
+scatter(papr-6.4,prob,80,'o') 
+hold on 
+scatter(papr-7.4,prob,80,'o') 
+%hold on
+%semilogy(papr_cubic1,prob_cubic1,'-^')
+%hold on
+%semilogy(papr_cubic2,prob_cubic2,'-r^')
+%hold on
+%semilogy(papr_cubic3,prob_cubic3,'-g^')
+legend([h1,h3,h4,h5,h2,h6,h7,h8],'OFDM','Precode','T-OFDM','SUSD Theroy','SUSD Simulation','OSID HDS','OSID MDS','OSID LDS')
+xlabel('D(dB)');
+ylabel('Prob(PAPR>D)');
\ No newline at end of file

testsemilogy2.m

+figure(1)
+h1=semilogy(papr,prob,'-o','linewidth',2,'markersize',6);
+axis([0,14,10^(-5),1])
+hold on
+h2=semilogy(papr1,prob1,'linewidth',1,'markersize',6);
+hold on
+h3=semilogy(papr2,prob2,'g','linewidth',1,'markersize',6)
+hold on
+h4=semilogy(papr3,prob3,'r','linewidth',1,'markersize',6)
+hold on
+h6=semilogy(papr_cubic1,prob_cubic1,'-^','linewidth',2,'markersize',6);
+hold on
+h7=semilogy(papr_cubic2,prob_cubic2,'-g^','linewidth',2,'markersize',6);
+hold on
+h8=semilogy(papr_cubic3,prob_cubic3,'-r^','linewidth',2,'markersize',6);
+%hold on
+%semilogy(papr_scfde,prob_scfde,'-x','linewidth',2,'markersize',6)
+hold on
+
+legend([h1,h2,h3,h4,h6,h7,h8],'OFDM','SUSD HDS','SUSD MDS','SUSD LDS','OSID HDS','OSID MDS','OSID LDS')
+
+xlabel('D(dB)');
+ylabel('Prob(PAPR>D)');
\ No newline at end of file

ze.m

+function y=ze(x)
+j=1;
+for n=1:length(x)
+    if abs(x(n))~=0
+        y(j)=x(n);
+        j=j+1;
+    end
+end
+
+if j==1
+    y=0;
+end