diff --git a/Classes/AWG.m b/Classes/AWG.m index ac1d637..2553933 100644 --- a/Classes/AWG.m +++ b/Classes/AWG.m @@ -63,7 +63,7 @@ classdef AWG % obj.fdac = 256e9; obj.dac_max = 0.5; obj.dac_min = -.5; - obj.f_cutoff = 80e9; + obj.f_cutoff = 800e9; end @@ -75,9 +75,12 @@ classdef AWG function signalclass_out = process(obj,signalclass_in) + + len_in = length(signalclass_in.signal); + % actual processing of the signal (steps 1. - 3.) signalclass_in.signal = obj.process_(signalclass_in.signal); - + % 4. Apply LPF on the signal if obj.lpf_active lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); @@ -93,6 +96,12 @@ classdef AWG % write to output signalclass_out = signalclass_in; + len_out = length(signalclass_out.signal); + + if len_out ~= len_in * obj.kover + warning("AWG: Output length maybe not correct.") + end + end function elec_out = process_(obj,data_in) @@ -129,7 +138,7 @@ classdef AWG elec_out = repmat(elec_out,obj.repetitions,obj.kover); elec_out = reshape(elec_out',[],1); - % 3. Add skew + % 3. Add skew (not working so far) if obj.skew_active elec_out = obj.skew(elec_out); end diff --git a/Classes/Scope.m b/Classes/Scope.m index 66f3d4c..23b1786 100644 --- a/Classes/Scope.m +++ b/Classes/Scope.m @@ -93,7 +93,8 @@ classdef Scope % sample signal % TODO: implement and test the delays. Also look for delay of % lpf filter - yout = obj.sampleSignal(xin); + %yout = obj.sampleSignal(xin); + yout = resample(xin,obj.fadc,obj.fsimu); % quantize signal yout = obj.quantize(yout); @@ -148,6 +149,8 @@ classdef Scope % Sample the signal using linear interpolation sampout = xin(floor(inx_re)) .* (1 - val_re) + xin(ceil(inx_re)) .* val_re ; + sampout2 = resample(xin,obj.fadc,obj.fsimu); + end function quantout = quantize(obj,xin) diff --git a/Functions/applyPulseShaping.m b/Functions/applyPulseShaping.m new file mode 100644 index 0000000..3b76620 --- /dev/null +++ b/Functions/applyPulseShaping.m @@ -0,0 +1,34 @@ +function yout = applyPulseShaping(xin,fsym,fdac) + +if ~rem(fdac,fsym) + %ist ein Vielfaches + sps = fdac / fsym; + up = sps; + dn = 1; +else + %ist kein Vielfaches + up = fdac / gcd(fdac, fsym); + dn = fsym / gcd(fdac, fsym); + sps= up; +end + +%Bau das Filter (hier rrc) +racos_len = 2048; +alpha = 0.1; +h = rcosdesign(alpha,racos_len,sps); + +%Apply Filter using Matlab build in fctn. +yout = upfirdn(xin,h,up,dn); + +%cut signal, which is longer due to fir filter +st = round(up/dn*racos_len/2); %we need to cut y_out +en = round(st + (length(xin)*up/dn) -1); + +yout = yout(st:en); + +%Check output integrity +if round(up/dn * length(xin)) ~= length(yout) + warning('Check signal length after pulse shaping'); +end + +end \ No newline at end of file diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m new file mode 100644 index 0000000..fe8d5e6 --- /dev/null +++ b/Functions/calc_ber.m @@ -0,0 +1,28 @@ +function [bits,errors,BER] = calc_ber(data_in,data_ref,skip) + +data_ref=logical(data_ref); +data_in = logical(data_in); +bits = 0; +errors=0; + +data_ref_overlap=zeros(size(data_ref,1),skip+(length(data_ref)-size(data_in,2))); +data_ref_pointer=0; +% Determine BER + +bits = bits+size(data_in,2)-skip; +try + errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end,:),2 ); +catch + %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') + errors = NaN; +end +try + errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end-1,:),2 ); +end +try + errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end-2,:),2 ); +end + +BER = sum(errors)/sum(bits); + +end \ No newline at end of file diff --git a/construct_classes.m b/construct_classes.m new file mode 100644 index 0000000..99a8091 --- /dev/null +++ b/construct_classes.m @@ -0,0 +1,2 @@ + + diff --git a/loop_simulation.m b/loop_simulation.m new file mode 100644 index 0000000..99a8091 --- /dev/null +++ b/loop_simulation.m @@ -0,0 +1,2 @@ + + diff --git a/main_simulation.m b/main_simulation.m deleted file mode 100644 index 8719fa2..0000000 --- a/main_simulation.m +++ /dev/null @@ -1,183 +0,0 @@ -%clear all - -O = 17; %order of prbs -N = 2^(O-1); %length of prbs -[~,seed] = prbs(O,1); %initialize first seed of prbs - -% Modulation -M = 4; %PAM-M -bitpattern = zeros(N,log2(M)); - -% Symbol Rate -fsym = 92e9; -% DAC Rate -fdac = 120e9; -% Simulation oversampling rate "k"; -kover = 16; -% ADC Rate -fadc = 256e9; -% Simulation frequency in "analog domain" -fsimu = kover * fdac ; - -%CONSTRUCTION -pam_mapper = PAMmapper(M,0); - -awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',80e9,'lpf_type',filtertypes.gaussian); - - -fil_tx = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); - -u_pi = 3.5; -vbias = (0.5*u_pi)-u_pi; -extmodlaser = EML("mode",emlmodes.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",0); - -amp = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power"); - -fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",17,"lambda0",1550); - -%noise loading -zeroDB_attenuator = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0); - -edfa_ = Amplifier("amp_mode","edfa_increase_nase","nase_mode","generate_ase","amplification_db",0,"noifig",33,"gain_mode","output_power"); - -rop_amplifier = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power","nase_mode","pass_ase"); - - -phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); - -fil_diode = Filter('filtdegree',1,"f_cutoff",70e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); - -scp = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",6,"quantbuffer",0.1); - -eq = EQ("K",2,"plottrain",0,"plotfinal",1,"training_length",2048,"Ne",[25,5,5],"Nb",[2,0,0],"training_loops",5,"DCmu",0.005,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"thres",[0.005 0.004 0.0005 ]); - -%SIMULATE - -% INFORMATION SIGNAL -for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); -end - -bits = Informationsignal(bitpattern); - -mod_out = pam_mapper.map(bits); -reference = mod_out; -mod_out.signal = applyPulseShaping(mod_out.signal,fsym,fdac); - -% ELECTRICAL DOMAIN -X = awg.process(mod_out); - -% X = fil_tx.process(X); -%X = fil_tx.process(X); - -X = X.normalize; - -% OPTICAL DOMAIN -X = extmodlaser.process(X); - -X = amp.process(X); - -X = fib.process(X); - -% X = zeroDB_attenuator.process(X); -% X = edfa_.process(X); -% X = rop_amplifier.process(X); - -X = phdiode.process(X); - -X = fil_diode.process(X); - -% ELECTRICAL DOMAIN -X = scp.process(X); - -X = X.resample("fs_out",2*fsym,"fs_in",fadc); - -% INFORMATION SIGNAL -X = X.normalize; -X = eq.process(X,reference); -rx_series = X.signal; - -X = pam_mapper.demap(X); - - -% BER -[bits,errors,BER] = calc_ber(X.signal(:,10000:end-20),bitpattern(10000:end-20,:)',0); - -disp(['BER: ', sprintf('%2E',BER)]); - - -figure() -hold on -plot(reference.signal); -plot(rx_series); -hold off - -disp(X.logbook); - -function yout = applyPulseShaping(xin,fsym,fdac) - -if ~rem(fdac,fsym) - %ist ein Vielfaches - sps = fdac / fsym; - up = sps; - dn = 1; -else - %ist kein Vielfaches - up = fdac / gcd(fdac, fsym); - dn = fsym / gcd(fdac, fsym); - sps= up; -end - -%Bau das Filter (hier rrc) -racos_len = 2048; -alpha = 0.1; -h = rcosdesign(alpha,racos_len,sps); - -%Apply Filter using Matlab build in fctn. -yout = upfirdn(xin,h,up,dn); - -%cut signal, which is longer due to fir filter -st = round(up/dn*racos_len/2); %we need to cut y_out -en = round(st + (length(xin)*up/dn) -1); - -yout = yout(st:en); - -%Check output integrity -if round(up/dn * length(xin)) ~= length(yout) - warning('Check signal length after pulse shaping'); -end - -end - -function [bits,errors,BER] = calc_ber(data_in,data_ref,skip) - -data_ref=logical(data_ref); -data_in = logical(data_in); -bits = 0; -errors=0; - -data_ref_overlap=zeros(size(data_ref,1),skip+(length(data_ref)-size(data_in,2))); -data_ref_pointer=0; -% Determine BER - -bits = bits+size(data_in,2)-skip; -try - errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end,:),2 ); -catch - %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') - errors = NaN; -end -try - errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end-1,:),2 ); -end -try - errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end-2,:),2 ); -end - -BER = sum(errors)/sum(bits); - -end - diff --git a/process_simulation.m b/process_simulation.m new file mode 100644 index 0000000..5fa3831 --- /dev/null +++ b/process_simulation.m @@ -0,0 +1,37 @@ + + + + +figure(22) +clf +subplot(2,1,1) +hold on +plot(reference.signal(10000:end-20)); +plot(rx_series(10000:end-20)); +hold off +subplot(2,1,2) +hold on +plot(reference.signal(4150:4175)); +plot(rx_series(4150:4175)); +hold off + +figure() +hold on +stem(X.signal(1,:)) +stem(bitpattern(:,1)') +hold off + +a = fftshift(xcorr(X.signal(1,:),circshift(bitpattern(:,1)',0))); + +% BER +if (length(X.signal) ~= length(bitpattern')) + warning("TX and RX bitstreams have different length...") +end + +[bits,errors,BER] = calc_ber(X.signal(:,10000:end-20),bitpattern(10000:end-19,:)',0); + +disp(X.logbook); +disp(['BER: ', sprintf('%2E',BER)]); + + + diff --git a/setup_simulation.m b/setup_simulation.m new file mode 100644 index 0000000..34b9c18 --- /dev/null +++ b/setup_simulation.m @@ -0,0 +1,192 @@ +clear all + +O = 18; %order of prbs +N = 2^(O-1); %length of prbs +[~,seed] = prbs(O,1); %initialize first seed of prbs + +% Modulation +M = 4; %PAM-M +bitpattern = zeros(N,log2(M)); + +% Symbol Rate +fsym = 112e9; +% DAC Rate +fdac = 160e9; +% Simulation oversampling rate "k"; +kover = 16; +% ADC Rate +fadc = 256e9; +% Simulation frequency in "analog domain" +fsimu = kover * fdac ; + +%% CONSTRUCT ALL CLASSES + +pam_mapper = PAMmapper(M,0); + +awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',0,'f_cutoff',80e9,'lpf_type',filtertypes.gaussian); + +fil_50 = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); + +u_pi = 3.5; +vbias = (0.5*u_pi)-u_pi; +extmodlaser = EML("mode",emlmodes.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",0); + +amp = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power"); + +fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",17,"lambda0",1550); + +% %noise loading +% zeroDB_attenuator = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0); +% +% edfa_ = Amplifier("amp_mode","edfa_increase_nase","nase_mode","generate_ase","amplification_db",0,"noifig",33,"gain_mode","output_power"); +% +% rop_amplifier = Amplifier("amp_mode","ideal_no_noise","amplification_db",-5,"gain_mode","output_power","nase_mode","pass_ase"); + +% X = zeroDB_attenuator.process(X); +% X = edfa_.process(X); +% X = rop_amplifier.process(X); + +phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + +fil_diode = Filter('filtdegree',1,"f_cutoff",120e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); + +scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",6,"quantbuffer",0.1); + +eq = EQ("K",2,"plottrain",0,"plotfinal",1,... + "training_length",4096*2,"training_loops",5,... + "Ne",[50,9,9],"Nb",[0,0,0],... + "DCmu",0.005,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0.000,... + "dd_loops",2,"epsilon",[10 100 1000 ],"M",4,... + "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0); + + +%% PROCESS +% INFORMATION SIGNAL +for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); +end + +bits = Informationsignal(bitpattern); + +mod_out = pam_mapper.map(bits); + +reference = mod_out; + +mod_out.signal = applyPulseShaping(mod_out.signal,fsym,fdac); + +% ELECTRICAL DOMAIN +X = awg.process(mod_out); + +X = X.normalize; + +% OPTICAL DOMAIN +X = extmodlaser.process(X); + +X = amp.process(X); + +X = fib.process(X); + +X = phdiode.process(X); + +X = fil_diode.process(X); + +% ELECTRICAL DOMAIN +% X = scp.process(X); + +X = X.resample("fs_out",2*fsym,"fs_in",fsimu); + +% INFORMATION SIGNAL +X = X.normalize; +eq_out = eq.process(X,reference); +rx_series = eq_out.signal; + +demap_out = pam_mapper.demap(eq_out); + +[bits,errors,BER] = calc_ber(demap_out.signal(:,10000:end-20),bitpattern(10000:end-19,:)',0); + + +%% +figure(22) +clf +subplot(3,1,1) +hold on +plot(reference.signal(10000:end-20)); +plot(rx_series(10000:end-20)); +hold off +subplot(3,1,2) +hold on +plot(reference.signal(4150:4175)); +plot(rx_series(4150:4175)); +hold off +subplot(3,1,3) +hold on +stem(demap_out.signal(1,4150:4175)) +stem(bitpattern(4150:4175,1)') +hold off +legend + +xax = 1:eq_out.length; +col = cbrewer2('Set1',4); + +figure(12) +sgtitle('Laser Linewidth = 10 MHz; SIR = 24 dB ; $N_{1,2}$ = 100') +subplot(1,4,1:2) +scatter(xax,eq_out.signal,4,'.','MarkerEdgeColor',col(1,:),'DisplayName','No MPI Mitigation'); +xlim([1, xax(end)]); +ylim([-2 2]); +xlabel('Sampling Index') +ylabel('Amplitude') +legend + +subplot(1,4,3:4) +scatter(1:X.length,X.signal,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','A1'); +xlim([1, xax(end)]); +ylim([-2 2]); +xlabel('Sampling Index') +ylabel('Amplitude') +legend +% +% subplot(1,4,3) +% scatter(xax,a2_out,4,'.','MarkerEdgeColor',col(3,:),'DisplayName','A2'); +% xlim([1, xax(end)]); +% ylim([-2 2]); +% xlabel('Sampling Index') +% ylabel('Amplitude') +% legend +% +% subplot(1,4,4) +% scatter(xax,a1_a2_out,4,'.','MarkerEdgeColor',col(4,:),'DisplayName','A1+A2'); +% xlim([1, xax(end)]); +% ylim([-2 2]); +% xlabel('Sampling Index') +% ylabel('Amplitude') +% legend + +% set(gcf,'Units','centimeters') +% set(gcf,'Position',[2 2 25 10]) + +a = fftshift(xcorr(X.signal(1,:),circshift(bitpattern(:,1)',0))); + +% BER +if (length(X.signal) ~= length(bitpattern')) + warning("TX and RX bitstreams have different length...") +end + + + +disp(demap_out.logbook); +disp(['BER: ', sprintf('%2E',BER)]); + + + + + + + + + + +