saddlepointappro.asv

一种新的光纤通信调制技术

function y = SaddlePointAppro(vtr_alpha,vtr_beta , xi_k, bitvalue) 
%% 由n_k的 MGF 求 P(n_k<xi_k, bitvalue作参考，得误码率)
	eta_nk = sum( vtr_alpha + vtr_beta);
	sigma_nk2 =   sum( vtr_beta.*(2*vtr_alpha + vtr_beta) );
   
	% 插入迭代, 求出 鞍点
	Posi_Bmax = max(vtr_beta); Nega_Bmax = -1*min(vtr_beta);
    1/Posi_Bmax ;
    -1/Nega_Bmax;
	if bitvalue == 0
		s0 = ( xi_k - eta_nk + sqrt( (xi_k - eta_nk).^2 + 4*sigma_nk2) )/(2*sigma_nk2);
        if s0 < -1/Nega_Bmax || s0 > 1/Posi_Bmax
            s0 = 0.99/Posi_Bmax;
        end
	else
        s0 = ( xi_k - eta_nk - sqrt( (xi_k - eta_nk).^2 + 4*sigma_nk2) )/(2*sigma_nk2);
        if s0 < -1/Nega_Bmax || s0 > 1/Posi_Bmax
            s0 = -0.99/Nega_Bmax;
        end
	end
    
	s0 = real( s0 ); s_n = s0; s_error = 1; n=1;
    while s_error > 1e-10                                        %% newton method 
        y1=Phi_nk_one(vtr_alpha, vtr_beta, s_n, xi_k);
		y2=Phi_nk_two(vtr_alpha, vtr_beta, s_n);
        if y2==0
            y2=eps;
        end
        y3=Phi_nk_three(vtr_alpha, vtr_beta, s_n);
        s_n1 = s_n - y1/y2; 
        q = y1*y3/y2.^2;
        s_error = abs( (s_n1-s_n)*q/(1-q) );
        s_n = s_n1;
        n = n+1;         % 计数器
    end
    u0=s_n         %%%%%%有值输出                                                  % -1.064e3
    if u0>1.0e10
        u0=1.0e10;
    else if u0<-
                    %%  For Test    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
					%     for n=-10:10
					%         x =u0*(1+n/100);
					%         s(n+11)  = x;
					%         y0(n+11) = Phi_nk( vtr_alpha, vtr_beta, x , xi_k);
					%         for m = -20:20
					%             y=  x+ i*u0*(m/100);
					%             map(n+11, m+21) = Phi_nk( vtr_alpha, vtr_beta, y , xi_k);
					%         end
					%     end
					%     figure (6)
					%     plot(s, real(y0) );
					%     figure(7)
					%     surf(abs(map));
                    
    phi_nk_two_u0 = Phi_nk_two(vtr_alpha, vtr_beta, u0); % 9.68e-7
    if phi_nk_two_u0==0
        phi_nk_two_u0=eps;
    end
    phi_nk_three_u0 = Phi_nk_three(vtr_alpha, vtr_beta, u0);            % 1.9e-9
    Kai = phi_nk_three_u0/phi_nk_two_u0/3 ;        % 169
       if (u0>0) && (Kai<0)
           Kai=0;
       end
       if isnan(Kai)||isinf(Kai)
           Kai=0;
       end
       Kai
    delta_v = 1/sqrt(2*phi_nk_two_u0);                                  % 718
    
    v = 0;
    f0 = InteFun(vtr_alpha, vtr_beta, u0, Kai, v, xi_k);
    Prob = 0.5*f0; p_error =1;
    while p_error>1e-8                   %%%  该循环用于搜索积分上边界
        v = v+delta_v;
        fn = InteFun(vtr_alpha, vtr_beta, u0, Kai, v, xi_k);
        Prob = Prob + fn;
        if Prob==0
            Prob=eps;
        end
        p_error =abs( fn/Prob );
    end 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  用于调试  %%%%%%%%%%%%%%%%%% 
%                     for n=1:10000
%                         x = v*n/10000;
%                         s(n)  = x;
%                         yt(n) = InteFun(vtr_alpha, vtr_beta, u0, Kai, x, xi_k);
%                     end
%                     figure (8)
%                     plot(s, yt );  
%                     
%                     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    积分部分trapzd    %%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%  integration program: qsimp(...), Numerical Recipes in C++ (2ndEd.), pp 144  %%%%%%
     Prob = 0 ; 
    JMAX = 15;          
    EPS = 1e-5;              
    ost = 0;                
    os = 0;
    for jj=1:JMAX
        if jj ==1                                                          %%% trapzd(... ),  pp 141, program begin,  %%%%%%
            fv = InteFun(vtr_alpha, vtr_beta, u0, Kai, v, xi_k);
            TrapzdS = 0.5*v*(f0 + fv);
        else
            it = 2^(jj-2);           
            tnm = it;
            del = v/tnm;           
            x = 0.5*del;            
            sss = 0; 
            for ss = 1:it
                sss = sss+ InteFun(vtr_alpha, vtr_beta, u0, Kai, x, xi_k);
                x = x + del;
            end 
            TrapzdS = 0.5 * (TrapzdS + del * sss ) ;
        end                                                                 %%% trapzd(... ),  pp 141, program end,  %%%%%%          
        st = TrapzdS;            
        s = (4*st - ost)/3.0;
        if jj>5
            if (abs(s-os) < EPS*abs(os) || (s == 0 & os == 0) )
                break;
            end
        end 
        os = s;        ost = st;
    end
    Prob =  abs( s/pi ) ; 
    %prob_formu = exp( Phi_nk(vtr_alpha, vtr_beta, u0, xi_k) )/sqrt( 2*pi*Phi_nk_two(vtr_alpha, vtr_beta, u0) )   %%做对比%%
    y = Prob;