ppp.m

来自「一种新的光纤通信调制技术」· M 代码 · 共 100 行

P=zeros(1,N);
for k=1:N
    if de(k)==1
       s=u01(k);%
   
       phai2=0;
       for I=1:2*M+1
           phai2=phai2+(2*alpha(I)*beta(I)+(beta(I)^2)*(1-beta(I)*s))/(1-beta(I)*s)^3;
       end
       phai2=phai2+1/(s^2);

       phai3=0;
       for I=1:2*M+1
            phai3=phai3+(6*alpha(I)*beta(I)^2+2*(beta(I)^3)*(1-beta(I)*s))/(1-beta(I)*s)^4;
       end
       phai3=phai3-2/(s^3);

       kin=0;%phai3/(3*phai2);
       
%
       vin=0;
       s=u01(k)+0.5*kin*vin^2+j*vin;%
            phai=0;
            for I=1:2*M+1
                phai=phai+alpha(I)*s/(1-beta(I)*s)-log(1-beta(I)*s);
            end
            phai=phai-log(s)-s*ek(k);
       fv0=real(exp(phai)*(1-j*kin*vin));
%求fv0



        detav=1/sqrt(2*phai2);
        sum0=0;
        n=1;
        fv=1;
        while(abs(fv)>1.0000e-003)
              vin=n*detav;
              s=u01(k)+0.5*kin*vin^2+j*vin;
            phai=0;
            for I=1:2*M+1
                phai=phai+alpha(I)*s/(1-beta(I)*s)-log(1-beta(I)*s);
            end
            phai=phai-log(s)-s*ek(k);
              fv=real(exp(phai)*(1-j*kin*vin));
              sum0=sum0+fv;
              n=n+1;
        end
        sum0=sum0+0.5*fv0;
        P(k)=sum0*detav/pi;
    else s=u02(k);
        phai2=0;
        for I=1:2*M+1
            phai2=phai2+(2*alpha(I)*beta(I)+(beta(I)^2)*(1-beta(I)*s))/(1-beta(I)*s)^3;
        end
        phai2=phai2+1/(s^2);

        phai3=0;
        for I=1:2*M+1
            phai3=phai3+(6*alpha(I)*beta(I)^2+2*(beta(I)^3)*(1-beta(I)*s))/(1-beta(I)*s)^4;
        end
        phai3=phai3-2/(s^3);

        kin=phai3/(3*phai2);
%
        vin=0;
        s=u02(k)+0.5*kin*vin^2+j*vin;%
        phai=0;
        for I=1:2*M+1
            phai=phai+alpha(I)*s/(1-beta(I)*s)-log(1-beta(I)*s);
        end
        phai=phai-log(s)-s*ek(k);
        fv0=real(exp(phai)*(1-j*kin*vin));
      %求fv0



       detav=1/sqrt(2*phai2);
       sum0=0;
       n=1;
       fv=1;
       while(abs(fv)>1.0000e-003)
            vin=n*detav;
            s=u02(k)+0.5*kin*vin^2+j*vin;
            phai=0;
            for I=1:2*M+1
                phai=phai+alpha(I)*s/(1-beta(I)*s)-log(1-beta(I)*s);
            end
            phai=phai-log(s)-s*ek(k);
            fv=real(exp(phai)*(1-j*kin*vin));
            sum0=sum0+fv;
            n=n+1;
       end
       sum0=sum0+0.5*fv0;
       P(k)=-sum0*detav/pi;
    end
end