split1.m

利用分步傅立叶积分法

% Nonlinear simulation model.
% only self phase modulation with dispersion is considered here
% ----------------------------------------------------------------------------
  nB = 64;			      % number of bits in the pattern
  npb = 64;			      % number of points per bit
  Tsample = npb*nB;		% total number of samples in time domain
  B = 4*2.48832e9;		% bit rate (# of bits/second)
  T = 1/B;			      % bit period
  Tstep = T/npb;		   % step size in time domain = period/number of points per bit
  Fsample = Tsample;    % Sampling rate in freq. domain
  FLoop = 1:Fsample;    % generate an array of frequency sample
  Fstep = B/nB;         % freq. sampling interval
  
 nf=Tsample;			   %how many sample points used in FFT
 fre=(-nf/2:nf/2-1)/Tstep/nf;	  %frequency scale
 time=(1:nf)*Tstep;		        %time scale
 omg=2*pi*fre;
 
%load data with chirp at the output of the external modulator.
a=rand(1,64)>0.5;
b=transpose(a)*ones(1,64);
cw=transpose(b);
d=reshape(cw,1,4096);
pppln=d';
%load pppln.prn     %signal optical intensity
fffln=ones(size(pppln))*0;  %no chirp considere dpp
%fiber dispersion parameters
%------------------------------------------------------------------
lam=1.553e-6;			   %signal wavelength (um)
s0=85;				      %fiber dispersion slope (s/m^3) 0.075ps/km/nm^2)
c=3e8;				      %lighe velocity
lam0=1320*1e-9;			%zero dispersion wavelength (m)
dlam=lam^2/c.*fre;		%wavelength shift due to signal baseband.
wl=lam+dlam;			   %actual wavelength range
D(1:nf)=s0/4*(wl(1:nf)-(lam0^4 ./(wl(1:nf).^3)));
				            %wavelength dependent dispersion factor
Ds=s0/4*(lam-(lam0^4 /(lam.^3)));
Ds0=Ds;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pwr=pppln;
phi=(fffln-mean(fffln))*0;
pave=mean(pppln);			%signal average power

ppr=sqrt(pwr).*cos(phi);    %re-construct signal (real part)
ppi=sqrt(pwr).*sin(phi);    %re-construct signal (imaginary part)

% Nonlinear effects of the fiber
%----------------------------------------------------------------------------
n2=2.35e-20; 			%Intensity dependent refractive index (m^2/W)
%Aeff=5.5e-11;			%effective fiber cross section (DSF)
Aeff=8e-11;			    %effective fiber cross section (SMF)
loss=0.25;			   %fiber loss/10 (dB/km)
alpha1=1e-3*2.303*(loss/10);		%fiber loss in Neper/m.

csp=ppr+i*ppi;	       % normalized complex optical field at input
ppr=(csp+conj(csp))/2;
ppi=imag((csp-conj(csp))/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cascade fiber spans

al=80;             	%fiber length in kilometer
l=al*1e3;				%fiber length in meter
pwr=(abs(ppr+i*ppi)).^2;
ppo=10;				                 %signal average power in dBm
power=10 .^(ppo/10)*1e-3/mean(pwr);	  %signal peak power in Watt

%calculate section lengths dl(mm)
%-----------------------------------------------------------------------------
mm=100;				%number of sections in the fiber
dl=l/mm;

%-----------------------------------------------------------------------------

   Input=sqrt(power)*(ppr+i*ppi);	% complex optical field at input
   csp=Input;


   for k=1:mm
	Hf=exp(-i*lam^2/(4*pi*c)*(2*pi*fre).^2.*D*dl);
	Hf=[Hf(nf/2+1:nf) Hf(1:nf/2)];
	sp=fft(csp)/(B*npb);
	sp=sp.*Hf';
	csp=ifft(sp)*B*npb;
	csb=abs(csp.*conj(csp));
	gggm=csb*(1-exp(-alpha1*dl))/alpha1*n2*2*pi/Aeff/lam;
	csp=csp.*exp(-i*gggm)*sqrt(exp(-alpha1*dl));
   end
   
   sp=fftshift(fft(csp)); %signal in frequency domain
   
   % A rise-cosine filter at the receiver
   fc=7.5e9; %filter 3-dB bandwidth in Hz.
   HHh=(1+exp(j*fre/fc*pi/2))/2;
   HHh(1:4096/2+1-98)=zeros(size(HHh(1:4096/2+1-98)));
   HHh(4096/2+1+98:4096)=zeros(size(HHh(4096/2+1+98:4096)));
   
   ssp=sp.*(HHh.^2)';
   csp=ifft(ssp);
   
   pwr=csp.*conj(csp);
   pwr=pwr/mean(pwr)/2;
   ppr=(csp+conj(csp))/2;
   ppi=imag((csp-conj(csp))/2);

csp=ppr+i*ppi;	% complex optical field at input
ppr=(csp+conj(csp))/2;
ppi=imag((csp-conj(csp))/2);
pwr=csp.*conj(csp);
pwr=pwr/mean(pwr)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pwo=pwr;
pwopp=sqrt(abs(pwo.*conj(pwo)));
pwopp=pwopp./mean(pwopp)/2;
Deye = eyediagram(pwopp,npb,nB);
axis([0 125 0 3])
  title(['Signal Eye Diagram, Pin=',num2str(ppo)]);




return