klse822.asv

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

     H2f=(exp(-i*4*pi*f*T)-exp(-i*pi/4)*exp(i*detaPHI1))/2;
     %Hf=exp(-i*2*pi^2*b2*Long*f.^2);
     Hf=exp(i*2*pi*f*(-0.5*detaTao)-i*2*pi^2*b2*Long*f.^2);%%% Hf+ 
     Hof=G*exp(-(f./(0.65*Bo)).^4);
     H1(l+N*R/2+1)=Hf*Hof*H1f;
     H2(l+N*R/2+1)=Hf*Hof*H2f;
 end
 Aw=fft(A);H1=fftshift(H1);
 Aw1=Aw.*H1;H1=ifftshift(H1);
 A1=ifft(Aw1);H2=fftshift(H2);
 Aw2=Aw.*H2;H2=ifftshift(H2);
 A2=ifft(Aw2);
 IA=A1.*conj(A1)-A2.*conj(A2);%%%光探测器输出的差分电流
 Aw=fft(IA);
 [b,a]=besself(5,2*pi*Be);
 w=((-N*R/2):(N*R/2-1))*2*pi/(2*N*T);
 Hes=freqs(b,a,w);
 Hes=fftshift(Hes);
 vtr_d_k=ifft(abs(Hes).*Aw);%%%完成对电流信号的电滤波
  
 %%求矩阵mtx_v%%%
 term=zeros(1,N*R);%%中间变量，用于暂时存储一个数组
 for ii=1:2*M
     f=(ii-M-1)/T0;
     Hof=G*exp(-(f./(0.65*Bo)).^4);
     H1f=Hof*(exp(-i*4*pi*f*T)+exp(-i*pi/4)*exp(i*detaPHI1))/2;
     H2f=Hof*(exp(-i*4*pi*f*T)-exp(-i*pi/4)*exp(i*detaPHI1))/2;
     
     term=Aw1*conj(H1f)-Aw2*conj(H2f);
     w=(-N*R/2:(N*R/2-1))*2*pi/(2*N*T)-(ii-M-1)*2*pi/T0;
     %w=([0:N*R/2-1 -N*R/2:-1])*2*pi/(N*T)-(ii-M-1)*2*pi/T0;
     Hens=freqs(b,a,w);
     mtx_v(ii,:)=ifft(term.*Hens);    
 end
 
 %%求矩阵mtx_A%%
[mtx_A]=MXA_up(M,T0,T,detaPHI1,Be,Bo,G);

%%%%%%求vtr_z and vtr_b%%%
[vtr_U1,vtr_eig]=eig(mtx_A);
[vtr_U,value]=MXU(vtr_U1,vtr_eig,mtx_A,M);
vtr_lammda=real(value);
mtx_b=vtr_U'*mtx_v;
sigma=sqrt(N0/T0/2);
%w_t=normrnd(0,sigma,1,2*M)+i*normrnd(0,sigma,1,2*M);
for ii=1:N*R
    w_t=normrnd(0,sigma,1,2*M)+i*normrnd(0,sigma,1,2*M);
    vtr_b=mtx_b(:,ii);
    vtr_z=vtr_U'*(w_t.');
    n_k(ii)=sum( vtr_lammda.*abs(vtr_z+vtr_b./vtr_lammda).^2);
    nu_k(ii)=sum( vtr_b.*conj(vtr_b)./vtr_lammda);
end
y_k=real(vtr_d_k+n_k-nu_k);                    %%%考虑PMD，求出y_k的正交分量





%eyediagram(y_k,2*R,2,R/2);
%figure(1);
%eyediagram(vtr_d_k,2*R,2,R/2);
%figure(2);
%Ana_Res=zeros(1,6);
%Ana_Res= SampleValueAnalyzer(y_k, N, R, KK1);
%OptTk=Ana_Res(1);%%%最佳判决时刻

%Decth=Ana_Res(2);%%%最佳门限，直接设为0

dth=0;%Decth;
xi_k=real( dth-vtr_d_k+nu_k);
%eyediagram(xi_k,R*2,2,R/2);
%figure(4)

%%%%%%%saddle point method 计算误码率
BER_up=zeros(1,N);
for ii=1:N
    ii
    jj=(ii-1)*R+16;
    xi=xi_k(jj);
    vtr_b=mtx_b(:,jj);
    vtr_alpha=(vtr_b.*conj(vtr_b))./vtr_lammda;
    vtr_beta=2*vtr_lammda*sigma^2;
    BER_up(ii)=SaddlePointAppro(vtr_alpha,vtr_beta,xi,1-a1(ii));
end

H1=zeros(0,N*R);
H2=zeros(0,N*R);
 for l=-N*R/2:(N*R/2-1)
     f=l/(2*N*T);%**********************************************************T change to 2*T
     H1f=(exp(-i*4*pi*f*T)+exp(i*pi/4)*exp(i*detaPHI2))/2;
     H2f=(exp(-i*4*pi*f*T)-exp(i*pi/4)*exp(i*detaPHI2))/2;
     Hf=exp(i*2*pi*f*(-0.5*detaTao)-i*2*pi^2*b2*Long*f.^2);%%% Hf+ 
     Hof=G*exp(-(f./(0.65*Bo)).^4);
     H1(l+N*R/2+1)=Hf*Hof*H1f;
     H2(l+N*R/2+1)=Hf*Hof*H2f;
 end
 Aw=fft(A);H1=fftshift(H1);
 Aw1=Aw.*H1;H1=ifftshift(H1);
 A1=ifft(Aw1);H2=fftshift(H2);
 Aw2=Aw.*H2;H2=ifftshift(H2);
 A2=ifft(Aw2);
 IA=A1.*conj(A1)-A2.*conj(A2);%%%光探测器输出的差分电流
 Aw=fft(IA);
 [b,a]=besself(5,2*pi*Be);
 w=((-N*R/2):(N*R/2-1))*2*pi/(2*N*T);
 Hes=freqs(b,a,w);
 Hes=fftshift(Hes);
 vtr_d_k=ifft(abs(Hes).*Aw);%%%完成对电流信号的电滤波
  
 %%求矩阵mtx_v%%%
 term=zeros(1,N*R);%%中间变量，用于暂时存储一个数组
 for ii=1:2*M
     f=(ii-M-1)/T0;
     Hof=G*exp(-(f./(0.65*Bo)).^4);
     H1f=Hof*(exp(-i*4*pi*f*T)+exp(i*pi/4)*exp(i*detaPHI2))/2;
     H2f=Hof*(exp(-i*4*pi*f*T)-exp(i*pi/4)*exp(i*detaPHI2))/2;
     
     term=Aw1*conj(H1f)-Aw2*conj(H2f);
     w=(-N*R/2:(N*R/2-1))*2*pi/(2*N*T)-(ii-M-1)*2*pi/T0;
     %w=([0:N*R/2-1 -N*R/2:-1])*2*pi/(N*T)-(ii-M-1)*2*pi/T0;
     Hens=freqs(b,a,w);
     mtx_v(ii,:)=ifft(term.*Hens);    
 end
 
 %%求矩阵mtx_A%%
[mtx_A]=MXA_down(M,T0,T,detaPHI2,Be,Bo,G);

%%%%%%求vtr_z and vtr_b%%%
[vtr_U1,vtr_eig]=eig(mtx_A);
[vtr_U,value]=MXU(vtr_U1,vtr_eig,mtx_A,M);
vtr_lammda=real(value);
mtx_b=vtr_U'*mtx_v;
sigma=sqrt(N0/T0/2);
%w_t=normrnd(0,sigma,1,2*M)+i*normrnd(0,sigma,1,2*M);
for ii=1:N*R
    w_t=normrnd(0,sigma,1,2*M)+i*normrnd(0,sigma,1,2*M);
    vtr_b=mtx_b(:,ii);
    vtr_z=vtr_U'*(w_t.');
    n_k(ii)=sum( vtr_lammda.*abs(vtr_z+vtr_b./vtr_lammda).^2);
    nu_k(ii)=sum( vtr_b.*conj(vtr_b)./vtr_lammda);
end
y_k=real(vtr_d_k+n_k-nu_k);              %%%

%eyediagram(y_k,2*R,2,R/2);
%figure(1);
%eyediagram(vtr_d_k,2*R,2,R/2);
%figure(2);
%Ana_Res=zeros(1,6);
%Ana_Res= SampleValueAnalyzer(y_k, N, R, b1);
%OptTk=Ana_Res(1);%%%最佳判决时刻

%Decth=Ana_Res(2);%%%最佳门限，直接设为0

dth=0;%Decth;
xi_k=real( dth-vtr_d_k+nu_k);
%eyediagram(xi_k,R*2,2,R/2);
%figure(4)

%%%%%%%saddle point method 计算误码率
BER_down=zeros(1,N);
for ii=1:N
    ii
    jj=(ii-1)*R+16;
    xi=xi_k(jj);
    vtr_b=mtx_b(:,jj);
    vtr_alpha=(vtr_b.*conj(vtr_b))./vtr_lammda;
    vtr_beta=2*vtr_lammda.*sigma^2;
    BER_down(ii)=SaddlePointAppro(vtr_alpha,vtr_beta,xi,1-b1(ii));
end



BER=0.5*(BER_up+BER_down);
pp=0;
NN=N;
for ii=1:N
    
    if BER(ii)>0.48||isinf(BER(ii))||isnan(BER(ii))
        BER(ii)=0;
        NN=NN-1;
    end
    
    pp=pp+BER(ii);
end
P_4RZ(oo+1)=pp/NN;
end
xx=11:15;
%plot(xx,log10(P_RZ),'b*');
%title('2-RZ-DPSK误码率性噪比仿真结果');
%xlabel('Eb/N0 (dB)');
%ylabel('log10(BER)');
p=polyfit(xx,log10(P_4RZ),5);
xx1=8:0.01:14;
plot(xx,log10(P_4RZ),'r*',xx1,polyval(p,xx1),'k-');
title('4-RZ-DPSK误码率性噪比仿真结果');
xlabel('Eb/N0 (dB)');
ylabel('log10(BER)');
grid;
save BER P_4RZ;
y=spline(log10(P_4RZ),xx,-9);

%              [m,n] = size(ber);
%              Minber = sum(ber')/n;                                           %% find the minimum ber and save it in file
%              Minber = Minber';
              
%              totalber = [ KK;ber];
%              save BER\TotalBER.dat  totalber   -ascii
%              save BER\Minber.dat  Minber    -ascii

             

% %***************************************************************************** 
% %数据存储
% fid1=fopen('Fiber_BER.dat','a');                         
%'a'表示创建一个新文件或删除已存在的文件内容，并进行数据读写操作。
% fprintf(fid1,'%e  %e\n',P_pin_av_preamplifier,BER); 
% status=fclose('all');              
% %************************************************************************