📄 klse822.txt
字号:
function y=klse822()
clear;
clc;
close all;
%*******************相关参数定义************************%
for oo=0:6
SNRdB=8+oo;
nsp=3; %理想放大器
hv=6.626*10^(-34)*1.94*10^(14); %h*f,为光子能量，f为波长是1550纳米的频率
N0=nsp*hv; %噪声功率谱密度
Eb=N0*10^(SNRdB/10); %信号平均功率，即每个码元的能量
N=2^5; %2DPSK设N=2^5,4DPSK设N=2^10
Rb=10^10; %10Gbit/s
T=1/Rb; %
Bo=1.1*Rb; %适用于NRZ－2－DPSK
Be=1.2*Rb; %同上
Miu=1.4;
Tao=2;
T0=Miu*(1/Bo+1/Be);
L=floor(Tao*N*Bo*T);
M=floor(Tao*Bo*T0);
%*****************色散相关参数定义******************************************
%load BER1 DD;
%D=DD*10^(-24); %s/nm.nm
D=0*10^(-24);
wavelength=1550; %nm
c=3e005; %km/s
b2=-D*wavelength^2/(2*pi*c); %s^2/km
%b2=2.0400e-023;%s^2/Km,b2=20.4ps^2/Km
Long=50;%Km
%****************phase error***********************************************
%load BER1 detaPHI;
detaPHI=0;
%%***************PMD相关参数************************************************
%load BER1 deta
deta=0;
Dp=1; %ps/sqrt(Km),PMD参数。
detaTao=deta*T;
%detaTao=0;%Dp*sqrt(Long)*10^(-12); %差分群时延DGD。
eta=0.5; %使两个PSP对称，使PMD最大
%***************************传输信号参数定义********************************
Bitrate=10*1e9; %信号传输速率定义
Samplerate=32; %每比特周期采样率定义
code_order=5; %伪随机二进制码阶数定义
Period=1/Bitrate;
Bitlength=2^code_order; % 2^code_order-1
V=Bitrate; %信号传输速率（10Gbit/s或者40Gbit/s)
T=Period; %比特周期（10Gbit/s系统为100ps;40Gbit/s系统为25ps)
N=Bitlength; %码长(127)
R=Samplerate; %每比特周期采样率（32）
T_delta=T/R; %取样时间间隔 %光纤放大器噪声指数
G=1; %光纤放大器固定增益（dB）
%************产生debrujin序列***********
fbconnection=[1 0 0 1 0];
KK1=de_sequence(fbconnection);
%**************差分编码******************
b=zeros(1,length(KK1));
b(1)=xor(KK1(1),0);
for ii=2:length(KK1)
b(ii)=xor(KK1(ii),b(ii-1));
end
KK=b;
%周期采样
jj=0;
for n=1:N
for T_Sample=-T/2+T_delta/2:T_delta:T/2
jj=jj+1;
UU(jj)=KK(n);
end
end
% figure(1)
% plot(UU)
% hold on
% title('伪随机码序列输出')
%Bessel滤波用以产生NRZ码
U_fft=fft(UU);
[b,a]=besself(10,2*pi*V*4); %10阶Bessel函数，截止频率=(Bitrate*4)GHz
w=0:2*pi/(N*T):((2*pi*R/T)-2*pi/(N*T));
h=freqs(b,a,w);
for kk=(N*R/2+1):N*R
h(kk)=conj(h(N*R+2-kk));
end
U=ifft(abs(h).*U_fft);
U=abs(U);
%%%%%%%调试
%U=UU; %直接使用矩形NRZ脉冲
for ii=1:N*R
if U(ii)==0
U(ii)=eps;
end
end
%********************************U为输出的PRBS
A_dpsk=sqrt(Eb/(2*T));
jj=0;
for nn=1:N
for tt=1:R
jj=jj+1;
fai(jj)=U(jj)*pi;
A(jj)=A_dpsk*exp(i*fai(jj));%%对相对码的绝对相移调制即是DPSK调制
end
end
% figure(2)
% plot(angle(A))
% title('激光器输出信号的瞬态相位')
%
% figure(3)
% plot(abs(fftshift(fft(A))).^2)
% title('激光器输出功率谱')
% figure(4);
%t = (1:N*R)/(N*T);
%B=A.*conj(A);
%plot(t,B)
%title('Optical Pulse Train')
%eyediagram(A, 2*R, 2, R/2);
% figure(1);
% eyescat(A,V/2,R*V,R/2); %调制器输出波形眼图
%title('光发射机输出信号眼图')
%****************************************DPSK****************************
%****************************yk+*****************************************
H1=zeros(0,N*R);
H2=zeros(0,N*R);
for l=-N*R/2:(N*R/2-1)
f=l/(N*T);
H1f=(exp(-i*2*pi*f*T)+exp(i*detaPHI))/2;
H2f=(exp(-i*2*pi*f*T)-exp(i*detaPHI))/2;
Hf=sqrt(eta)*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/(N*T);
Hes=freqs(b,a,w);
Hes=fftshift(Hes);
vtr_d_k1=ifft(abs(Hes).*Aw);%%
%*****************************yk-****************************************
H1=zeros(0,N*R);
H2=zeros(0,N*R);
for l=-N*R/2:(N*R/2-1)
f=l/(N*T);
H1f=(exp(-i*2*pi*f*T)+exp(i*detaPHI))/2;
H2f=(exp(-i*2*pi*f*T)-exp(i*detaPHI))/2;
Hf=sqrt(1-eta)*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/(N*T);
Hes=freqs(b,a,w);
Hes=fftshift(Hes);
vtr_d_k2=ifft(abs(Hes).*Aw);%%
vtr_d_k=vtr_d_k1+vtr_d_k2;
%%求矩阵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=((exp(-i*2*pi*f*T)+exp(i*detaPHI))/2)*Hof;
H2f=((exp(-i*2*pi*f*T)-exp(i*detaPHI))/2)*Hof;
term=Aw1*conj(H1f)-Aw2*conj(H2f);
w=(-N*R/2:(N*R/2-1))*2*pi/(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(M,T0,T,detaPHI,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, 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=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(ii)=SaddlePointAppro(vtr_alpha,vtr_beta,xi,1-KK1(ii));%****1-KK1(ii)
end
pp=0;
NN=N;
for ii=1:N
if BER(ii)>1||isinf(BER(ii))||isnan(BER(ii))
BER(ii)=0;
NN=NN-1;
end
pp=pp+BER(ii);
end
P_2dpsk_NRZ(oo+1)=pp/NN;
end
xx=8:14;
p=polyfit(xx,log10(P_2dpsk_NRZ),5);
xx1=8:0.01:14;
plot(xx,log10(P_2dpsk_NRZ),'r*',xx1,polyval(p,xx1),'b-');
title('2-NRZ-DPSK误码率性噪比仿真结果');
xlabel('Eb/N0 (dB)');
ylabel('log10(BER)');
grid;
save BER P_2dpsk_NRZ;
y=spline(log10(P_2dpsk_NRZ),xx,-9);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -