📄 klse822.m
字号:
function y=klse()
%***************************************************
%输出电信号误码率计算
%***************************************************
clear;
clc;
close all;
%****************
% 系统参数定义
%****************
%***************************************************************************************************************************
%输入光功率（平均值）（峰值为：2*平均值）和用于误码测试的可调光衰减器
P_input_average=1e-3; %单位（mW）
att=20; %单位（-dB）
%*****************************************************************************************************************************
%*****************************************************************************************************************************
%传输信号参数定义
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; %取样时间间隔
F_sample=linspace(-R*V/2,+R*V/2-V/N,N*R); %取样频率分量序列(2*pi*F_sample才为角频率序列)
%*******************************************************************************************************************************
%*******************************************************************************************************************************
%激光器和调制器参数设置
P_max=P_input_average; %激光器发射的峰值功率（mW）
P_ext=30; %调制器消光比（dB）
P_min=P_max*10^(-P_ext/10); %调制器0码输出光功率
chirp=0; %调制器啁啾系数
laser_width=10e6; %激光器线宽（MHz）
%********************************************************************************************************************************
%********************************************************************************************************************************
%% 光信号部分 %光发射机模型（包括编码、采样、滤波和调制）
%************************************************************************
%伪随机二进制序列（PRBS）编码
for j=1:code_order
DD(j)=round(rand(1)); %四舍五入取整
end
for ii=1:N
KK(ii)=mod((DD(code_order)+DD(code_order-4)),2);
DDD=DD;
for jj=1:code_order
if jj==1
DD(jj)=KK(ii);
else
DD(jj)=DDD(jj-1);
end
end
end
%周期采样
j=0;
for n=1:N
for T_Sample=-T/2+T_delta/2:T_delta:T/2
j=j+1;
UU(j)=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);
%激光器发射功率、线宽、啁啾和调制器的消光比
%激光器输出：A(t)=sqrt(P(t))*exp(j*fai(t))
%光源的相位噪声fai(t)包括两个部分：1.光源线宽造成的，对应于频率起伏；2.光源调制造成的（光功率变化）
delta_t=T/R; %取样时间间隔
average_fai=sqrt(2*pi*delta_t*laser_width); %取样时间间隔之间相位变化的均方根
j=1;
fai(1)=0;
A(1)=sqrt(U(1)*P_max+P_min);
for nn=1:N %对于所有码
for tt=1:R %对于一个周期
j=j+1;
if j<N*R
fai(j)=fai(j-1)+randn*average_fai-chirp/4*((U(j+1)-U(j-1))/U(j)); %光源相位噪声fai(t)的两个部分
A(j)=sqrt(U(j)*P_max+P_min)*exp(fai(j)*i);
elseif j==N*R
fai(j)=fai(j-1)+randn*average_fai-chirp/4*((U(1)-U(j-1))/U(j)); %光源相位噪声fai(t)的两个部分
A(j)=sqrt(U(j)*P_max+P_min)*exp(fai(j)*i);
end
end
end
% figure(2)
% plot(angle(A))
% title('激光器输出信号的瞬态相位')
%
% figure(3)
% plot(abs(fftshift(fft(A))).^2)
% title('激光器输出功率谱')
t = (1:N*R)/(N*T);
B=A.*conj(A);
%plot(t,B)
%title('Optical Pulse Train')
% eyediagram(B, 2*R, 2, R/2)
% figure(5)
% eyescat(B,V/2,R*V,R/2); %调制器输出波形眼图
% title('光发射机输出信号眼图')
%*******************************************************************************************************************************
%光纤放大器（EDFA）和光接收机相关参数定义
NF=5; %光纤放大器噪声指数
G=23.5; %光纤放大器固定增益（dB）
c=2.9979*1e8;
lamda = 1550*1e-9; % 波长 m
q=1.6e-19; %电子电荷电量
hp=6.63e-34; %Plank常数
Vs=c/lamda; %光频率 1e14 Hz
Bo=4*10e9; %光滤波器带宽
Re=1; %光探测器的响应度
ST=400e-24; %% 热噪声频谱密度
I_dark=1e-9; %% dark current (1nA)
Band_rate=0.75; %% ratio of receiver bandwidth with signal rate;
Be=Band_rate*V; %% electrical filter bandwidth, receiver bandwidth:(NRZ: 0.75B; RZ: 0.5B)
N0 =10.^(NF/10) * hp * Vs * (10.^(G/10) - 1) ; %/ (10.^(G/10));
%N0 = 0;
mu = 1.38;
T0 = mu*(1/Bo + 1/Be);
M = 2.^( ceil( log2( T0/T_delta )));
M=19;
T0 = M*T_delta;
mu = T0/(1/Bo + 1/Be);
%%%%%%%%%%%%%%% Opt. Filter
[b,a]=besself(10,pi*Bo); %% 10th order Bessel function, with cutoff freq. Bo
w=(-N*R/2:(N*R/2-1))*2*pi/(N*T);
wn = (-M:(M-1))*2*pi/T0;
Hos=freqs(b,a,w); %% filtering function for signal
Hon=freqs(b,a,wn); %% filtering function for noise
vtr_H = diag( abs(Hon) );
Aw = fft(A); %% optical signal filtering
Hos = fftshift(Hos);
A = abs(Hos).*Aw; %% How to Calculate the signal-noise ratio ? 信噪比如何计算
wd = (12.5e9); %% method one
ntmp = ceil( wd *N*T/2 );
A1 = A(1:ntmp); A2= A( (N*R-ntmp+1):(N*R));
Eb = ( sum(A1.*conj(A1)) + sum(A2.*conj(A2)) ) / 2/ntmp;
PASE= N0*wd;
SNR_Freq = 10*log10(Eb/PASE);
A = ifft(A);
Hos = ifftshift(Hos);
Eb = sum(A.*conj(A))/(N*R); %% method two
PASE= N0*Bo;
SNR2_time = 10*log10(Eb/PASE);
Signal = A; %% optical signal back up 备份
%%%%%%%%%%%%%%%%%%% Receiver %%%%%%%%%%%%%%%%%%%%%%%%%%%%
IA = Re * A.*conj(A); %% square law detection
[b,a]=besself(10,2*pi*Be); %% Bessel Electrical Filter, 10-th order, Cutoff freq = Be
w=((-N*R/2):(N*R/2-1))*2*pi/(N*T);
wn = (-M:(M-1))*2*pi/T0;
Hes=freqs(b,a,w);
Hen=freqs(b,a,wn);
w=1:2*M;
I = ones(1,2*M);
wq = ( (I.')*w - (w.')*I )*2*pi/T0;
vtr_Q = freqs(b,a,wq); %% Q, Hertimitian symmetric matrix
mtx_A = vtr_H' * vtr_Q * vtr_H;
[vtr_U, vtr_eig] = eig(mtx_A); %% precision value: 10^-16. mtx_A*vtr_U - vtr_U * vtr_eig, or, mtx_A - vtr_U * vtr_eig * vtr_U';
vtr_U*vtr_U;
vtr_lammda=real(diag(vtr_eig));
Aw = fft(IA);
Hes= fftshift(Hes);
vtr_d_k = ifft(abs(Hes).*Aw); %% electrical filter for electrical signal
Hes = ifftshift(Hes);
mtx_v=zeros(2*M,N*R); %% vi
S_l = fft(Signal) ;
for ii=1:2*M
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(S_l.*Hens);
end
mtx_b = vtr_U'*vtr_H'* mtx_v; %% vector b
sigma = sqrt( N0/T0/2 );
for ii=1:N*R %% every sample point of signal
w_t= normrnd( 0, sigma, 1, M+M) + i * normrnd(0, sigma, 1, M+M);
%w_f = fftshift( fft(w_t) );
vtr_z = vtr_U' * w_t.';
vtr_b = mtx_b(:,ii);
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 ); %% d_k + n_k - nu_k
%% eyediagram(y_k,Samplerate*2,2,Samplerate/2);
%% find the point with maximum Q value, next function return value: [OptTk, DecTh, mu1,mu0,sigma1,sigma0]
Ana_Res= SampleValueAnalyzer(y_k, N, R, KK);
OptTk = Ana_Res(1); DecTh = Ana_Res(2); SPStat=Ana_Res(3:6);
dDecTh = ( Ana_Res(3) - Ana_Res(4) - Ana_Res(5) -Ana_Res(6) ) ;
ber=[];
%for kk = 1:5:21
dth = DecTh + (11-11)*dDecTh/50 ; %% decision threshold
xi_k = real( dth*ones(size(A)) - vtr_d_k + nu_k );
eyediagram(xi_k,Samplerate*2,2,Samplerate/2);
for ii=1:N %% P(e_k|{a_n}) = P(y_t_k<gamma_th) = P(n_k<zeta_k)
jj = (ii-1)*R+OptTk;
xi = xi_k(jj); %% gamma_th - d_k + nu_k
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, KK(ii)) ;
end
ber = [ber; BER]; %% mean value of BER
%end
[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');
% %*****************************************************************************
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -