diff --git a/Classes/02_optical/dp_fiber_lib/CNLSE.m b/Classes/02_optical/dp_fiber_lib/CNLSE.m new file mode 100644 index 0000000..d2773dc --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/CNLSE.m @@ -0,0 +1,135 @@ + +function [opt_out_struct,state] = CNLSE(opt_in_struct,state) + + % init transfer functions h.X and h.Y + h = struct('X',0,'Y',0); + state.common_beta=struct('X',0,'Y',0); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % pre calculations + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % calculate transfer function and rotate coordines for both + % polarizations + + for n=1:2 + % get current polarization name and contrary one + curPol = state.polNames{n}; + + % extend linear transfer function depending on beta values for the + % current polarization + for n_beta = 1:length(state.beta.(curPol)) +% h.(curPol) = h.(curPol) - 1j*state.beta.(curPol)(n_beta)*(state.omega).^(n_beta-1)/factorial(n_beta-1); +% if n_beta ~= 2 + state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); +% end + end + + opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; + end + + state.h=h; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Splitstep method + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% state.SS_dzs = zeros(1,state.max_nonlin_its); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Split Step Method + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % get nonlinear step size + [state.dz] = getNLstepsize(state,opt_out_struct); + + state.n_step = 0; + state.z_prop = 0; + state.test_dz = []; + state.powers = []; + + while state.z_prop < state.L + + + + + if state.z_prop + state.dz > state.L + state.dz = state.L - state.z_prop; + end + + + % lin conv + % opt_out_struct = [ opt_out_struct 0 0 0 0 0 ]; + % opt_out_struct + %%%%%%%%%%%%% + % STEP + % update step number + + state.n_step=state.n_step+1; + state.dzs(state.n_step)=state.dz; + + % half linear step + [opt_out_struct,state] = lin_step(state,opt_out_struct,state.dz/2); + + + % complete nonlinear step + [opt_out_struct,state] = nl_step(state,opt_out_struct,state.dz); + + % half linear step + [opt_out_struct,state] = lin_step(state,opt_out_struct,state.dz/2); + + %%%%%%%%%%%%% + % prepare next STEP + + % overlap(n_step+1,:) = opt_out_struct(M+1:end); + % opt_out_struct = opt_out_struct(1:M); + + + % get nonlinear step size + [state.dz] = getNLstepsize(state,opt_out_struct); + + end + + +% figure(88);clf;subplot(2,1,1);stem(state.test_plates);subplot(2,1,1); hold all;stem(-1000*state.test_plate_numbers);subplot(2,1,2);stem(state.dzs) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Post Calculations + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% opt_out_struct.X.envelope = ( cos(state.psi)*cos(state.chi) + 1j*sin(state.psi)*sin(state.chi))*opt_out_struct.X + ... +% (-sin(state.psi)*cos(state.chi) - 1j*cos(state.psi)*sin(state.chi))*opt_out_struct.Y; +% + buffer.X.envelope = opt_out_struct.X; + buffer.X.type = opt_in_struct.X.type; + buffer.X.wavelength = opt_in_struct.X.wavelength; + if isfield(buffer.X,'Nase') + buffer.X.Nase = opt_in_struct.X.Nase; + else + buffer.X.Nase = 0; + end + + opt_out_struct.X =[]; + opt_out_struct.X.envelope = buffer.X.envelope; + opt_out_struct.X.type = buffer.X.type; + opt_out_struct.X.wavelength = buffer.X.wavelength; + opt_out_struct.X.Nase = buffer.X.Nase; + +% +% opt_out_struct.Y.envelope = ( sin(state.psi)*cos(state.chi) - 1j*cos(state.psi)*sin(state.chi))*opt_out_struct.X.envelope + ... +% ( cos(state.psi)*cos(state.chi) - 1j*sin(state.psi)*sin(state.chi))*opt_out_struct.Y; +% + buffer.Y.envelope = opt_out_struct.Y; + buffer.Y.type = opt_in_struct.Y.type; + buffer.Y.wavelength = opt_in_struct.Y.wavelength; + if isfield(buffer.Y,'Nase') + buffer.Y.Nase = opt_in_struct.Y.Nase; + else + buffer.Y.Nase = 0; + end + + opt_out_struct.Y =[]; + opt_out_struct.Y.envelope = buffer.Y.envelope; + opt_out_struct.Y.type = buffer.Y.type; + opt_out_struct.Y.wavelength = buffer.Y.wavelength; + opt_out_struct.Y.Nase = buffer.Y.Nase; + +% figure(100+loop);plot([real(opt_out_struct.X.envelope);real(opt_out_struct.Y.envelope)].'); +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m new file mode 100644 index 0000000..cfa7294 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m @@ -0,0 +1,120 @@ + +function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) + + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % pre calculations + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + state.common_beta=struct('X',0,'Y',0); + + for n=1:2 + % get current polarization name and contrary one + curPol = state.polNames{n}; + + % extend linear transfer function depending on beta values for the current polarization + % Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41] + for n_beta = 1:length(state.beta.(curPol)) + state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); + end + + %opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; + end + + beta_const = state.beta.('X')(1); + beta_1 = state.beta.('X')(2); + beta_2 = state.beta.('X')(3); + beta_3 = state.beta.('X')(4); + deltaomega = state.omega; + beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3; + +% opt_in_x = gpuArray(opt_in_x); +% opt_in_y = gpuArray(opt_in_y); + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Split Step Method + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% [opt_out_x,opt_out_y] = split_step_loop(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... +% state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... +% state.chi,state.manakov,state.beat_len); + +% [opt_out_x,opt_out_y] = split_step_loop_mex(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... +% state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... +% state.chi,state.manakov,state.beat_len); + + % get nonlinear step size + [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); + %[state.dz] = getNLstepsize_original(state,opt_out_struct); + + state.n_step = 0; + state.z_prop = 0; + state.test_dz = []; + state.powers = []; + + tic + + while state.z_prop < state.L + + % reduce step length (dz) if we are to overshoot the fiber length + % (L) in the next step + if state.z_prop + state.dz > state.L + state.dz = state.L - state.z_prop; + end + + % update step number (n) + state.n_step=state.n_step+1; + + % append current step length to logbook (dzs) + state.dzs(state.n_step)=state.dz; + + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... + state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... + state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... + = lin_step(... + opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + state.dz/2,state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... + state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... + state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); + + % complete nonlinear step + + [opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y); + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... + state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... + state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... + = lin_step... + (opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + state.dz/2,state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... + state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... + state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); + + % get nonlinear step size + [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); + %[state.dz] = getNLstepsize_original(state,opt_out_struct); + + + + end + + toc + + + + + + opt_out_x = (opt_in_x); + opt_out_y = (opt_in_y); + +% opt_out_x = gather(opt_in_x); +% opt_out_y = gather(opt_in_y); + +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/getNLstepsize.m b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m new file mode 100644 index 0000000..a70edb9 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m @@ -0,0 +1,25 @@ + +function [rDZ] = getNLstepsize(ux,uy,gamma,dzmin,dzmax,dphimax,alpha_lin) + + + maxPow = max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2)); + + Leff = dphimax/maxPow; + alpha_lin = max([alpha_lin.X alpha_lin.Y]); + nl_att_len_ratio = alpha_lin*Leff; + + if nl_att_len_ratio >= 1 + rDZ = dzmax; + else + if alpha_lin == 0 + step = Leff; + else + %effective length? + step = -1/alpha_lin*log(1-nl_att_len_ratio); + end + + rDZ = min([step dzmax]); + rDZ = max([rDZ dzmin]); + end + +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/getNLstepsize_original.m b/Classes/02_optical/dp_fiber_lib/getNLstepsize_original.m new file mode 100644 index 0000000..93a3d63 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/getNLstepsize_original.m @@ -0,0 +1,27 @@ + +function [rDZ] = getNLstepsize_original(state,aOpt) + + ux = aOpt.X; + uy = aOpt.Y; + + maxPow = max(state.gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2)); + + Leff = state.SS_dphimax/maxPow; + alpha_lin = max([state.alpha_lin.X state.alpha_lin.Y]); + nl_att_len_ratio = alpha_lin*Leff; + + if nl_att_len_ratio >= 1 + rDZ = state.SS_dzmax; + else + if alpha_lin == 0 + step = Leff; + else + %effective length? + step = -1/alpha_lin*log(1-nl_att_len_ratio); + end + + rDZ = min([step state.SS_dzmax]); + rDZ = max([rDZ state.SS_dzmin]); + end + +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/lin_step.m b/Classes/02_optical/dp_fiber_lib/lin_step.m new file mode 100644 index 0000000..eea8a1e --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/lin_step.m @@ -0,0 +1,125 @@ + +%function [rOpt,state] = lin_step(state,aOpt,aStepSize) + +function [rOpt_x,rOpt_y,z_prop,lin_z_test,... + corr_length,n_plates_done,missing_dz,n_step,test_plates,... + test_plate_numbers,brf,common_beta_x,common_beta_y,alpha_lin_x,alpha_lin_y]... + = lin_step(... + opt_x,opt_y,z_prop,lin_z_test,aStepSize,corr_length,... + n_plates_done,missing_dz,n_step,test_plates,test_plate_numbers,... + brf,common_beta_x,common_beta_y,alpha_lin_x,alpha_lin_y) + +%%%%%% 1) Update and Check Distances etc. %%%%%% + +% update propgated distance z_prop +z_prop = z_prop + aStepSize; + +% calculate the number of plates needed for the so far propagated fiber length +n_plates = ceil(z_prop/corr_length); + +% subtract the number of plates which were already processed +n_plates_left = n_plates - n_plates_done; + +% compute last plate size ( if it fits, it should be 0) +if missing_dz > aStepSize + + last_plate = aStepSize; + missing_dz = missing_dz-aStepSize; + plate_sizes = last_plate; + plate_numbers = n_plates; + +else + + last_plate = aStepSize - missing_dz - (n_plates_left-1)*corr_length; + + if missing_dz == 0 + missing_dz = []; + end + + %build vector of plate lengths with missing plate part from prev. + %iterartion , then some normal plates and finally a fraction of a plate + %to fit into the step length + plate_sizes = [missing_dz corr_length*ones(1,n_plates_left-1) last_plate]; + + if n_plates_done == 0 + plate_numbers =[(n_plates_done+1):(n_plates-1) n_plates]; + else + plate_numbers = [n_plates_done (n_plates_done+1):(n_plates-1) n_plates]; % not wrking yet + end + + %remember for next step + missing_dz = corr_length - last_plate; + +end + +plate_steps = repmat(n_step,1,length(plate_sizes)); + +% +%figure;stem(plate_sizes); +test_plates = [test_plates,plate_sizes]; + +test_plate_numbers = [test_plate_numbers, plate_numbers]; + +%%%%%% 2) Apply Waveplate Model %%%%%% + +% transfer optical envelope to frequency domain for effective convolution with transfer function h +opt_x=fft(opt_x); +opt_y=fft(opt_y); + +% db1 = gpuArray(brf.db1); +% db0 = gpuArray(brf.db0); +% common_beta_x = gpuArray(common_beta_x); +% common_beta_y = gpuArray(common_beta_y); + +db1 = (brf.db1); +db0 = (brf.db0); +common_beta_x = (common_beta_x); +common_beta_y = (common_beta_y); + +% process every waveplate with given sizes in plate_sizes +for n=1:length(plate_sizes) + dz = plate_sizes(n); + + % figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150))); + % MOV1=[MOV1 getframe(87)]; + + % extract rotation matrix from pre calculated matrices + matR = brf.matR{plate_numbers(n)}; + + % transform to eigenvalue of of fiber segment + tOpt.X = conj(matR(1,1))*opt_x + conj(matR(2,1))*opt_y; + tOpt.Y = conj(matR(1,2))*opt_x + conj(matR(2,2))*opt_y; + + % calculate statistical delta beta for pmd + delta_beta = 0.5*(db1+db0(n))/corr_length; + % build transfer function with delta beta + %common.beta = beta1+beta2*omega^2 + + %accumulate delta beta for log... + brf.simdgd = brf.simdgd + (db1(length(db1)/2+1)+db0(n))/corr_length; + + h.X = exp(-1j*(common_beta_x-delta_beta)*dz); + h.Y = exp(-1j*(common_beta_y+delta_beta)*dz); + % delta_beta has to be added to the transfer function + + % process with transfer function + tOpt.X = h.X.*tOpt.X ; + tOpt.Y = h.Y.*tOpt.Y ; + + % rotate back + opt_x = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y; + opt_y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y; + +end + +lin_z_test = lin_z_test + sum(plate_sizes,2); + +%update the number of processed plates so far +n_plates_done = n_plates_done + n_plates_left; + +% attanuate the signal each linear state with alpha +% ( 0.2dB = 4.6052e-05 ) +rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); % /2 not sure why (have to find it in formulas) +rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); % but not relevant for now + +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/lin_step_original.m b/Classes/02_optical/dp_fiber_lib/lin_step_original.m new file mode 100644 index 0000000..dac5b32 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/lin_step_original.m @@ -0,0 +1,106 @@ + +function [rOpt,state] = lin_step_original(state,aOpt,aStepSize) + + + % update propgated distance z_prop + state.z_prop = state.z_prop + aStepSize; + +% if state.synchronous_plates % not waveplate model (just rotation with dz) +% state.plate_sizes = aStepSize; +% +% else + % calculate the number of plates needed for the so far propagated fiber + % length + state.n_plates = ceil(state.z_prop/state.corr_length); + + % subtract the number of plates which were already be processed + state.n_plates_left = state.n_plates - state.n_plates_done; + + % compute last plate size ( if it fits, it should be 0) + if state.missing_dz > aStepSize + + state.last_plate = aStepSize; + state.missing_dz = state.missing_dz-aStepSize; + state.plate_sizes = state.last_plate; + state.plate_numbers = state.n_plates; + + else + + state.last_plate = aStepSize - state.missing_dz - (state.n_plates_left-1)*state.corr_length; + + if state.missing_dz == 0 + state.missing_dz = []; + end + + state.plate_sizes = [state.missing_dz state.corr_length*ones(1,state.n_plates_left-1) state.last_plate]; + + if state.n_plates_done == 0 + state.plate_numbers =[(state.n_plates_done+1):(state.n_plates-1) state.n_plates]; + else + state.plate_numbers = [state.n_plates_done (state.n_plates_done+1):(state.n_plates-1) state.n_plates]; % not wrking yet + end + + state.missing_dz = state.corr_length - state.last_plate; + end + + state.plate_steps = repmat(state.n_step,1,length(state.plate_sizes)); + + % figure;stem(state.plate_sizes); + state.test_plates = [state.test_plates,state.plate_sizes]; + state.test_plate_numbers = [state.test_plate_numbers, state.plate_numbers]; +% end + + % transfer optical envelope to frequency domain for effective + % convolution with transfer function h + aOpt.X=fft(aOpt.X); + aOpt.Y=fft(aOpt.Y); + + % process every waveplate with given sizes in state.plate_sizes + for n=1:length(state.plate_sizes) + dz = state.plate_sizes(n); + +% figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150))); +% state.MOV1=[state.MOV1 getframe(87)]; + + % extract rotation matrix from pre calculated matrices + matR = state.brf.matR{state.plate_numbers(n)}; + + % transform to eigenvalue of of fiber segment + tOpt.X = conj(matR(1,1))*aOpt.X + conj(matR(2,1))*aOpt.Y; + tOpt.Y = conj(matR(1,2))*aOpt.X + conj(matR(2,2))*aOpt.Y; + + % calculate statistical delta beta for pmd + delta_beta = 0.5*(state.brf.db1+state.brf.db0(n))/state.corr_length; + % db1 = sqrt(3*pi/8)*(para.dgd/para.fa)/state.wave_plates.*state.omega; +% delta_beta = 0.5*(state.brf.db0(n))/state.corr_length; + + % build transfer function with delta beta + % common.beta = beta1+beta2*omega^2 + % delta beta + h.X = exp(-1j*(state.common_beta.X-delta_beta)*dz); + h.Y = exp(-1j*(state.common_beta.Y+delta_beta)*dz); + % delta_beta has to be added to the transfer function + + % process with transfer function + tOpt.X = h.X.*tOpt.X ; + tOpt.Y = h.Y.*tOpt.Y ; + + % rotate back + aOpt.X = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y; + aOpt.Y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y; + + end + + state.lin_z_test = state.lin_z_test + sum(state.plate_sizes,2); + + %update the number of processed plates so far + state.n_plates_done = state.n_plates_done + state.n_plates_left; + + + % attanuate the signal each linear state with alpha + % ( 0.2dB = 4.6052e-05 ) + rOpt.X=ifft(exp(-state.alpha_lin.X*aStepSize/2).*aOpt.X); % /2 not sure why (have to find it in formulas) + rOpt.Y=ifft(exp(-state.alpha_lin.Y*aStepSize/2).*aOpt.Y); % but not relevant for now + + +end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/nl_step.m b/Classes/02_optical/dp_fiber_lib/nl_step.m new file mode 100644 index 0000000..800c502 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/nl_step.m @@ -0,0 +1,39 @@ + +%function [rOpt,state] = nl_step(state,aOpt,aDz) + +function [rOpt_x,rOpt_y] = nl_step(opt_x,opt_y, dz, gamma, chi, use_manakov, beatlength, alpha_lin_x, alpha_lin_y) + + if ~use_manakov % CNLSE + + rOpt_x = opt_x .* exp( (-1j*(1/3)*gamma*dz).* ... + ( (2 + cos(2*chi)^2)*(abs(opt_x).^2) + ... + (2+2*sin(2*chi)^2)*(abs(opt_y).^2) ) ); + + rOpt_y = opt_y .* exp( (-1j*(1/3)*gamma*dz).* ... + ( (2 + cos(2*chi)^2)*(abs(opt_y).^2) + ... + (2+2*sin(2*chi)^2)*(abs(opt_x).^2) ) ); + +% A_x = opt_x; +% A_y = opt_y; +% +% rOpt_x = 1i* gamma * (abs(A_x).^2 + (2/3 .* abs(A_y).^2) ) .* A_x + ((1i * gamma / 3) * conj(A_x).*(A_y.^2) * exp(-2i * dz * 2*pi / beatlength )); +% rOpt_y = 1i* gamma * (abs(A_y).^2 + (2/3 .* abs(A_x).^2) ) .* A_y + ((1i * gamma / 3) * conj(A_y).*(A_x.^2) * exp(-2i * dz * 2*pi / beatlength )); + + + else + % estimate effective length of dz (ref?) + if (alpha_lin_x == 0) && (alpha_lin_y == 0) + Leff = dz; + else + Leff = (1-exp(-alpha_lin_x*dz))/alpha_lin_x; + end + + %compute power + power = real(opt_x).^2+imag(opt_x).^2+real(opt_y).^2+imag(opt_y).^2; +% power= abs(opt_x).^2+abs(opt_y).^2; +% powers = [powers;power]; + Hnl = exp( -1j*8/9*gamma*power*Leff); + rOpt_x = opt_x .* Hnl; + rOpt_y = opt_y .* Hnl; + end +end diff --git a/Classes/02_optical/dp_fiber_lib/nl_step_original.m b/Classes/02_optical/dp_fiber_lib/nl_step_original.m new file mode 100644 index 0000000..4fb9593 --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/nl_step_original.m @@ -0,0 +1,30 @@ + +function [rOpt,state] = nl_step(state,aOpt,aDz) + + if ~state.manakov % CNLSE + + rOpt.X = aOpt.X .* exp( (-1j*(1/3)*state.gamma*aDz).* ... + ( (2 + cos(2*state.chi)^2)*(abs(aOpt.X).^2) + ... + (2+2*sin(2*state.chi)^2)*(abs(aOpt.Y).^2) ) ); + + rOpt.Y = aOpt.Y .* exp( (-1j*(1/3)*state.gamma*aDz).* ... + ( (2 + cos(2*state.chi)^2)*(abs(aOpt.Y).^2) + ... + (2+2*sin(2*state.chi)^2)*(abs(aOpt.X).^2) ) ); + + else + % estimate effective length of dz (ref?) + if (state.alpha_lin.X == 0) && (state.alpha_lin.Y == 0) + Leff = aDz; + else + Leff = (1-exp(-state.alpha_lin.X*aDz))/state.alpha_lin.X; + end + + %compute power + power = real(aOpt.X).^2+imag(aOpt.X).^2+real(aOpt.Y).^2+imag(aOpt.Y).^2; +% power= abs(aOpt.X).^2+abs(aOpt.Y).^2; +% state.powers = [state.powers;power]; + Hnl = exp( -1j*8/9*state.gamma*power*Leff); + rOpt.X = aOpt.X .* Hnl; + rOpt.Y = aOpt.Y .* Hnl; + end +end diff --git a/Classes/02_optical/dp_fiber_lib/split_step_loop.m b/Classes/02_optical/dp_fiber_lib/split_step_loop.m new file mode 100644 index 0000000..e74f5bf --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/split_step_loop.m @@ -0,0 +1,93 @@ +function [opt_x,opt_y] = split_step_loop(L,opt_x,opt_y,gamma,SS_dzmin,SS_dzmax,SS_dphimax,alpha_lin,... + lin_z_test,corr_length,n_plates_done,missing_dz,brf,common_beta,... + chi,manakov,beat_len) + +%SPLIT_STEP_LOOP Summary of this function goes here +% Detailed explanation goes here + %Optical Input +% opt_x; +% opt_y; +% +% %required for loop condition +% z_prop = 0; +% L; +% +% %required for NLstepsize +% gamma; +% SS_dzmin; +% SS_dzmax; +% SS_dphimax; +% alpha_lin; +% +% %required for lin_step +% z_prop; +% lin_z_test; +% corr_length; +% n_plates_done; +% missing_dz; +% n_step = 0; +% brf; +% common_beta.X; +% common_beta.Y; +% alpha_lin.X; +% alpha_lin.X; +% +% %required fr nonlin step +% chi; +% manakov; +% beat_len ; +% alpha_lin.X; +% alpha_lin.Y; + + + % get nonlinear step size + [dz] = getNLstepsize(opt_x,opt_y,gamma,SS_dzmin,SS_dzmax,SS_dphimax,alpha_lin); + + n_step = 0; + z_prop = 0; + + while z_prop < L + + % reduce step length (dz) if we are to overshoot the fiber length + % (L) in the next step + if z_prop + dz > L + dz = L - z_prop; + end + + % update step number (n) + n_step=n_step+1; + + % half linear step + [opt_x,opt_y,z_prop,lin_z_test,... + corr_length,n_plates_done,missing_dz,n_step,... + brf,common_beta.X,... + common_beta.Y,alpha_lin.X,alpha_lin.X]... + = lin_step(... + opt_x,opt_y,z_prop,lin_z_test,... + dz/2,corr_length,n_plates_done,missing_dz,n_step,... + brf,common_beta.X,... + common_beta.Y,alpha_lin.X,alpha_lin.X); + + % complete nonlinear step + + [opt_x,opt_y] = nl_step(opt_x,opt_y, dz, gamma, chi, manakov, beat_len ,alpha_lin.X, alpha_lin.Y); + + % half linear step + [opt_x,opt_y,z_prop,lin_z_test,... + corr_length,n_plates_done,missing_dz,n_step,... + brf,common_beta.X,... + common_beta.Y,alpha_lin.X,alpha_lin.X]... + = lin_step... + (opt_x,opt_y,z_prop,lin_z_test,... + dz/2,corr_length,n_plates_done,missing_dz,n_step,... + brf,common_beta.X,... + common_beta.Y,alpha_lin.X,alpha_lin.X); + + % get nonlinear step size + [dz] = getNLstepsize(opt_x,opt_y,gamma,SS_dzmin,SS_dzmax,SS_dphimax,alpha_lin); + + end + + +end + diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m index 14db6dd..37a39fa 100644 --- a/projects/WDM/WDM_model.m +++ b/projects/WDM/WDM_model.m @@ -21,8 +21,8 @@ laser_linewidth = 0e6; % EQ SETTINGS vnle_order1 = 50; -vnle_order2 = 0; -vnle_order3 = 0; +vnle_order2 = 3; +vnle_order3 = 3; vnle_order=[vnle_order1,vnle_order2,vnle_order3]; dfe_order = [0 0 0]; len_tr = 4096*2; @@ -34,15 +34,15 @@ mu_dc = 0.005; mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; mu_dfe = 0.0004; - -rcalpha = 0.05; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - +%DB Stuff db_precode = 0; db_encode = 0; duob_mode = db_mode.no_db; apply_pulsef = 0; +rcalpha = 0.05; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + N = numel(wavelengthplan); f_plan = physconst('lightspeed')./(wavelengthplan.*1e-9); @@ -67,6 +67,7 @@ num_realiz = 1; gmi_vnle_bitwise = NaN(length(wavelengthplan),length(rop),num_realiz); snr_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); ber_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); +output = cell(length(wavelengthplan),length(rop),num_realiz); for realiz = 1:num_realiz @@ -166,30 +167,46 @@ for realiz = 1:num_realiz Rx_sig = Scpe_cell{1}; Rx_sig = Rx_sig.normalize("mode","rms"); - % FFE or VNLE - eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); - - [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols{l}); - showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); - [mi_gomez] = calc_air(eq_signal_sd, Symbols{l}, "skip_front", 100, "skip_end", 100); - [gmi_vnle_bitwise(l,ri,realiz)] = calc_ngmi(eq_signal_sd,Symbols{l}); - % [gmi_bitwise_2] = calc_gmi_bitwise(eq_signal_sd,Symbols{l}); - snr_vnle(l,ri,realiz) = calc_snr(Symbols{l}, eq_signal_sd-Symbols{l}); - - % eq_signal_sd.plot("displayname",'bla','fignum',199); - % eq_signal_sd.eye(fsym,M,"fignum",103837); - - eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); - rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); - [~,tot_err,ber_vnle(l,ri,realiz),a] = calc_ber(rx_bits.signal,Tx_bits{l}.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - burst_vnle = count_error_bursts(a, 10)./tot_err; - - % showLevelConfusionMatrix(eq_signal_hd,Symbols{l},"M",M,"fignum",200,"displayname",'bla'); - % showLevelScatter(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','f_sym',fsym,'fignum',201); - % show2Dconstellation(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','fignum',2241); - - fprintf('CH %d :BER VNLE: %.2e \n',l,ber_vnle(l,ri,realiz)); - fprintf('CH %d :NGMI VNLE: %.2f \n',l,gmi_vnle_bitwise(l,ri,realiz)./m); + + + + ffe_order = [50, 0, 0]; + eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + ffe_results = ffe(eq_ffe,M,Rx_sig,Symbols{l},Tx_bits{l},... + "precode_mode",duob_mode,... + 'showAnalysis',0,... + "postFFE",[],... + "eth_style_symbol_mapping",0); + + output{l,ri,realiz} = ffe_results; + + + + % % FFE or VNLE + % eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + % + % [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols{l}); + % showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); + % [mi_gomez] = calc_air(eq_signal_sd, Symbols{l}, "skip_front", 100, "skip_end", 100); + % [gmi_vnle_bitwise(l,ri,realiz)] = calc_ngmi(eq_signal_sd,Symbols{l}); + % % [gmi_bitwise_2] = calc_gmi_bitwise(eq_signal_sd,Symbols{l}); + % snr_vnle(l,ri,realiz) = calc_snr(Symbols{l}, eq_signal_sd-Symbols{l}); + % + % % eq_signal_sd.plot("displayname",'bla','fignum',199); + % % eq_signal_sd.eye(fsym,M,"fignum",103837); + % + % eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); + % rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); + % [~,tot_err,ber_vnle(l,ri,realiz),a] = calc_ber(rx_bits.signal,Tx_bits{l}.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % burst_vnle = count_error_bursts(a, 10)./tot_err; + % + % % showLevelConfusionMatrix(eq_signal_hd,Symbols{l},"M",M,"fignum",200,"displayname",'bla'); + % % showLevelScatter(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','f_sym',fsym,'fignum',201); + % % show2Dconstellation(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','fignum',2241); + % + % fprintf('CH %d :BER VNLE: %.2e \n',l,ber_vnle(l,ri,realiz)); + % fprintf('CH %d :NGMI VNLE: %.2f \n',l,gmi_vnle_bitwise(l,ri,realiz)./m); end diff --git a/projects/WDM/WDM_settings.m b/projects/WDM/WDM_settings.m index 49ddaf4..a0bfe85 100644 --- a/projects/WDM/WDM_settings.m +++ b/projects/WDM/WDM_settings.m @@ -1,3 +1,24 @@ + + + + + + + + + + + + + + + + + + + + + % Add the imdd_simulation framework to the path if ispc addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));