光纤通信matlab-simulink仿真.rar

% This file solves the nonlinear Schrodinger equation for % pulse propagation in an optical fiber using the split-step % Fourier method described in: % % Agrawal, Govind. Nonlinear Fiber Optics, 2nd ed. Academic % Press, 1995, Chapter 2 % % The following effects are included in the model: group % velocity dispersion (GVD), GVD-slope / third-order % dispersion, loss, and self-phase modulation (n2). The core % routine for implementing the split-step propagation is % called ssprop_matlabfunction_modified.m clear all close all % DQPSK - 2bits per symbol --> Baudrate = 1/2 Bitrate % CONSTANTS c = 299792458; % speed of light (m/s) % NUMERICAL PARAMETERS P0 = 0.01 % peak power (W) % Pulsewidth of the Bit --> BitRate FWHM = 25 % pulse width FWHM (ps) % Define the change %halfwidth = FWHM/1.6651 % for Gaussian pulse halfwidth = FWHM % for square pulse % DQPSK - 2bits per symbol --> Baudrate = 1/2 Bitrate BitRate = 1/halfwidth; % GHz BaudRate = BitRate/2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%************************ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % for DPSK Vpi=5; halfVpi = Vpi/2; twoVpi=Vpi*2; samplingtime = FWHM/1e3; decisiontime = (FWHM/samplingtime)/2; %index of the output element in the matrix taken for detection simtime = 10; tout = 0; % threshold =0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % nt = 2^8; % number of points in FFT PRBSlength = 2^7; % Make sure : FFT time window (=nt*dt) = PRBSlength * FWHM... % FFTlength nt = PRBSlength/block * numbersamples/bit = PRBSlength * (FWHM/dt) % num_samplesperbit = FWHM/dt should be about 8 - 16 samples/bit num_samplesperbit = 32; % should be 2^n dt = FWHM/num_samplesperbit; dt_baudrate = (1/(BaudRate))/num_samplesperbit ; % sampling time(ps); % time step (ps) nt = PRBSlength*num_samplesperbit; % FFT length % dt = 1.5; dz = 0.5; % distance stepsize (km) nz = 200; % number of z-steps maxiter = 10; % max # of iterations tol = 1e-5; % error tolerance % OPTICAL PARAMETERS nonlinearthreshold = 0.005; % 5mW -- % Nonlinear Threshold Peak Power lambda = 1550; % wavelength (nm) optical_carrier = c/(lambda*1e-9); %dBperkm = 0.2; % loss (dB/km) alpha_indB = 0.2; % loss (dB/km) D = 17; % GVD (ps/nm.km); if anomalous dispersion(for compensation),D is negative beta3 = 0.3; % GVD slope (ps^3/km) ng = 1.46; % group index n2 = 2.6e-20; % nonlinear index (m^2/W) Aeff = 47; % effective area (um^2) % CALCULATED QUANTITIES T = nt*dt; % FFT window size (ps) -Agrawal: should be about 10-20 times of the pulse width alpha_loss = log(10)*alpha_indB/10; % alpha (1/km) beta2 = 1000*D*lambda^2/(2*pi*c); % beta2 (ps^2/km); %---------------------------------------------------- % beta 3 can be calculated from the Slope Dispersion (S) as follows:] % Slope Dispersion % S = 0.092; % ps/(nm^2.km) % beta31 = (S - (4*pi*c./lambda.^3))./(2*pi*c./lambda.^2) %---------------------------------------------------- gamma = 2e24*pi*n2/(lambda*Aeff); % nonlinearity coef (km^-1.W^-1) t = ((1:nt)'-(nt+1)/2)*dt; % vector of t values (ps) t1 = [(-nt/2+1:0)]'*dt; % vector of t values (ps) t2 = [(1:nt/2)]'*dt; % vector of t values (ps) w = 2*pi*[(0:nt/2-1),(-nt/2:-1)]'/T; % vector of w values (rad/ps) v = 1000*[(0:nt/2-1),(-nt/2:-1)]'/T; % vector of v values (GHz) vs = fftshift(v); % swap halves for plotting v_tmp = 1000*[(-nt/2:nt/2-1)]'/T; % STARTING FIELD % P0 = 0.001 % peak power (W) % FWHM = 20 % pulse width FWHM (ps) %halfwidth = FWHM/1.6651 % for Gaussian pulse %For square wave input, the FWHM = Half Width %halfwidth = FWHM; L = nz*dz Lnl = 1/(P0*gamma) % nonlinear length (km) Ld = halfwidth^2/abs(beta2) % dispersion length (km) N = sqrt(abs(Ld./Lnl)) % governing the which one is dominating: dispersion or Non-linearities ratio_LandLd = L/Ld % if L << Ld --> NO Dispersion Effect ratio_LandLnl = L/Lnl % if L << Lnl --> NO Nonlinear Effect % Monitor the broadening of the pulse with relative the Dispersion Length % Calculate the expected pulsewidth of the output pulse % Eq 3.2.10 in Agrawal "Nonlinear Fiber Optics" 2001 pp67 FWHM_new = FWHM*sqrt(1 + (L/Ld)^2) % N<<1 --> GVD ; N >>1 ---> SPM Leff = (1 - exp(-alpha_loss*L))/alpha_loss expected_normPout = exp(-alpha_loss*2*L) NlnPhaseshiftmax = gamma*P0*Leff betap = [0 0 beta2 beta3]'; % Constants for ASE of EDFA % PSD of ASE: N(at carrier freq) = 2*h*fc*nsp*(G-1) with nsp = Noise % Figure/2 (assume saturated gain) %**************** Standdard Constant ******************************** h = 6.626068e-34; %Plank's Constant %****************************************** %************************** %************************** % u0 = sqrt(P0).*prbs_18May(PRBSlength,FWHM,dt); %************************** %************************** % %PROPAGATE % tmp = cputime; % tic % fftw('planner','exhaustive'); % %PROPAGATE % tmp = cputime; % tic %fftw('planner','exhaustive'); % % Nonlinear Threshold Peak Power % P_non_thres = 0.005; % 5mW % % FORMAT of function u1 % %function u1 = ssprop_modified(u0,dt,dz,nz,alpha_indB,betap,gamma,P_non_thres,maxiter,tol); % u = ssprop_modified(u0,dt,dz,nz,alpha_indB,betap,gamma,P_non_thres,maxiter); % u = ssprop_modified(u,dt,dz,20,alpha_indB,[0 0 -beta2*100/20 -beta3]',gamma,P_non_thres,maxiter); % %u = ssprop_modified(u,dt,dz,nz,alpha_indB,-betap,gamma,P_non_thres,maxiter); % % for i = 1: num_span % % u_tmp = ssprop(u,dt,dz,nz,alpha_indB,betap,gamma,maxiter); % % u = u_tmp + wgn(size(u_tmp,1),1,1e-5,'linear','complex'); % % end % fprintf(1,'ssprop completed (%.2f s)\n',cputime-tmp); % toc % % % PLOT OUTPUT % % figure(1); % iv = find(abs(t)<370); % plot (t(iv),abs(u0(iv)).^2/P0,t(iv),abs(u(iv)).^2/P0); % %plot (t(iv),abs(u(iv)).^2/P0); % grid on; % xlabel ('t (ps)'); % ylabel ('|u(z,t)|^2 (W)'); % title ('Initial and Final Pulse Shapes'); % % % figure(2); % iv = find(abs(vs)<280); % S0 = fftshift(abs(dt*fft(u0)/sqrt(2*pi)).^2); % S = fftshift(abs(dt*fft(u)/sqrt(2*pi)).^2); % plot (vs(iv),S(iv)); % grid on; % xlabel ('\nu (GHz)'); % ylabel ('|U(z,\nu)|^2'); % title ('Initial and Final Pulse Spectra');