From 8b2129ac037737db5cb13e3b5840d648bcd86857 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Thu, 25 May 2023 16:29:27 +0200 Subject: [PATCH] =?UTF-8?q?WIP=20B=C3=BCro?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Classes/AWG.m | 9 +- Classes/EML.m | 4 +- Classes/Fiber.m | 7 +- Classes/Filter.m | 35 ++++- Classes/Opticalsignal.m | 1 + Classes/Signal.m | 108 ++++++++++++- Datatypes/{emlmodes.m => eml_mode.m} | 2 +- Functions/applyPulseShaping.m | 4 +- setup_simulation.m | 220 ++++++++++++++------------- 9 files changed, 270 insertions(+), 120 deletions(-) rename Datatypes/{emlmodes.m => eml_mode.m} (78%) diff --git a/Classes/AWG.m b/Classes/AWG.m index abead51..c1bf453 100644 --- a/Classes/AWG.m +++ b/Classes/AWG.m @@ -4,6 +4,7 @@ classdef AWG properties(Access=public) + preset kover %oversampling factor e.g. 16 repetitions %repeat the signal to generate a longer sequence? fdac %needed @@ -64,13 +65,8 @@ classdef AWG obj.dac_max = 0.5; obj.dac_min = -.5; obj.f_cutoff = 80e9; - end - - - - end function signalclass_out = process(obj,signalclass_in) @@ -91,7 +87,8 @@ classdef AWG signalclass_in = Electricalsignal(signalclass_in,"fs",obj.fdac*obj.kover,"logbook",signalclass_in.logbook); % append to logbook - signalclass_in = signalclass_in.logbookentry(); + lbdesc = ['AWG preset ', obj.preset, 'k_over:',num2str(obj.kover),'. f_dac:',num2str(obj.fdac*1e-9),'GHz. Resolution:',num2str(obj.bit_resolution),' bits.']; + signalclass_in = signalclass_in.logbookentry(lbdesc); % write to output signalclass_out = signalclass_in; diff --git a/Classes/EML.m b/Classes/EML.m index e3e533d..9eda123 100644 --- a/Classes/EML.m +++ b/Classes/EML.m @@ -32,7 +32,7 @@ classdef EML %EML Construct an instance of this class % Detailed explanation goes here arguments - options.mode; + options.mode eml_mode = eml_mode.im_cosinus ; options.fsimu; options.lambda; options.power; @@ -70,7 +70,7 @@ classdef EML signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0); % append to logbook - lbdesc = ['EML ']; + lbdesc = [num2str(obj.lambda),' nm Laser with ',num2str(obj.power),' dBm P_out. Linew.=',num2str(obj.linewidth*1e-6),' MHz. Modulation mode: ',char(obj.mode) ]; signalclass_in = signalclass_in.logbookentry(lbdesc); % write to output diff --git a/Classes/Fiber.m b/Classes/Fiber.m index 9f0c3c4..bda66a6 100644 --- a/Classes/Fiber.m +++ b/Classes/Fiber.m @@ -73,8 +73,11 @@ classdef Fiber obj.linstep = -obj.alpha_lin/2 - 2*1j*pi^2*obj.b2*faxis.^2 - 4/3*1j*pi^3*obj.b3*faxis.^3; - opt_out = obj.NLSE(opt_in); - + if obj.gamma ~= 0 + opt_out = obj.NLSE(opt_in); + else + opt_out = ifft(fft(opt_in).*exp(obj.linstep*obj.fiber_length)); % only one linear step + end %attenuate nase end diff --git a/Classes/Filter.m b/Classes/Filter.m index c5fdb2a..f6164c6 100644 --- a/Classes/Filter.m +++ b/Classes/Filter.m @@ -151,14 +151,38 @@ classdef Filter faxis=linspace(-obj.fsamp/2,obj.fsamp/2,blocklen+1)';%generates arow vector faxis of blocklen+1 points linearly spaced between and including -para.fs/2 and para.fs/2 faxis=ifftshift(faxis(1:end-1)); - + H=exp(-((faxis)/(obj.f_cutoff*2)).^(2*obj.filtdegree)*log(2)*2^(2*obj.filtdegree-1)); + +% figure() +% hold on +% xline(obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--'); +% xline(-obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--'); +% plot(faxis*1e-9,20*log10(abs(H)),'LineWidth',3); +% ax = gca; +% ylim([-6 0]) +% grid on +% xlabel('Freq in GHz') +% ylabel('Magnitude (dB)') + end % Build Filter from coefficients if filterType ~= 10 [H,w] = freqz(B, A, obj.signal_length,'whole'); - + +% figure() +% hold on +% plot(w/(2*pi)*obj.fsamp*1e-9,20*log10(abs(H)),'LineWidth',3); +% ax = gca; +% ylim([-6 0]) +% grid on +% xlabel('Freq in GHz') +% ylabel('Magnitude (dB)') + + + % freqz(B, A) + %hfvt = fvtool(B,A); end end @@ -167,7 +191,12 @@ classdef Filter function show(obj) obj.signal_length = 1024; [H,w] = obj.buildFilter(obj.filterType); - plot(w,20*log10(abs(H))) + plot(w/pi,20*log10(abs(H))); + ax = gca; + ax.XTick = 0:.5:2; + grid on + xlabel('Normalized Frequency (\times\pi rad/sample)') + ylabel('Magnitude (dB)') %freqz(H); end end diff --git a/Classes/Opticalsignal.m b/Classes/Opticalsignal.m index 8b58ea9..1fa5d94 100644 --- a/Classes/Opticalsignal.m +++ b/Classes/Opticalsignal.m @@ -42,6 +42,7 @@ classdef Opticalsignal < Signal function pow = power(obj) pow = mean(abs(obj.signal.^2)) ; + % pow = pow2db(pow)+30; end end diff --git a/Classes/Signal.m b/Classes/Signal.m index 1240422..7d76852 100644 --- a/Classes/Signal.m +++ b/Classes/Signal.m @@ -24,6 +24,7 @@ classdef Signal end + %% CONVERT TO INFORMATIONSIGNAL function [i_sig, varargout] = Informationsignal(obj) if isa(obj,'Electricalsignal') @@ -40,6 +41,7 @@ classdef Signal end + %% CONVERT TO Electricalsignal function [e_sig, varargout] = Electricalsignal(obj,options) arguments @@ -70,6 +72,7 @@ classdef Signal end + %% CONVERT TO Opticalsignal function o_sig = Opticalsignal(obj, options) arguments @@ -91,9 +94,9 @@ classdef Signal if isa(obj,'Electricalsignal') %convert to optical - + o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase); - + elseif isa(obj,'Informationsignal') @@ -104,12 +107,14 @@ classdef Signal end + %% Display length function return_length = length(obj) %METHOD1 Summary of this method goes here % Detailed explanation goes here return_length = length(obj.signal); end + %% Write Logbook Entry function obj = logbookentry(obj,varargin) if nargin > 1 @@ -130,6 +135,7 @@ classdef Signal end + %% Resample Signal function obj = resample(obj,options) arguments @@ -146,6 +152,104 @@ classdef Signal end + %% + function spectrum(obj,fsamp,options) + + arguments + obj + fsamp + options.figurename = []; + options.displayname = []; + end + + %Get figure if there is already a spectrum plot -> I want to add the new + %spectum "onto" the existing plot to have a better comparison + if isempty(options.figurename) + fig = findall(groot, 'Type', 'figure', 'Name', 'power density'); + if isvalid(fig) + fig = get(fig); + hold on + else + figure('name','power density'); + end + else + fig = findall(groot, 'Type', 'figure', 'Name', options.figurename); + if isvalid(fig) + fig = get(fig); + hold on + else + figure('name',options.figurename); + end + + end + + + %compute FFT of input + Fsignal = fft(obj.signal); + + %POWER spectral density (todo: toggle?) + psd = Fsignal.*conj(Fsignal); + + %Use only magnitude of FFT (which was complex) + psd = abs(psd); + + %Shift the spectrum to yield + psd = fftshift(psd); + + %divide by N + psd = psd/length(Fsignal); + + %smoothing + psd = smooth(psd,100); + + psd_plot = 20*log10(psd); + + + testParseval = 1; + if testParseval == 1 + E_FreqDomain = sum(psd); + %test parseval + E_TimeDomain = sum(abs(Fsignal.^2)); + + if isequal(round(E_FreqDomain,1),round(E_TimeDomain,1)) + %disp('Parseval is right!'); + else + disp('Parseval theorem is not right...'); + end + end + + if fsamp <= 1e+100 + %Frequency Axis + freq_vec = linspace(-fsamp/2,fsamp/2,length(psd)); + freq_vec = reshape(freq_vec,size(psd_plot)); + + + if ~isempty(options.displayname) + plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname); + else + plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5); + end + + xlabel('Frequency [GHz]') + else + %Wavelength Axis + freq_vec = physconst('LightSpeed')*linspace(-fsamp/2,fsamp/2,length(psd))./((physconst('LightSpeed')/1550e-9)^2); + if ~isempty(options.displayname) + plot(freq_vec*1e9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname); + else + plot(freq_vec*1e9,psd_plot,'Linewidth',0.5); + end + + xlabel('Wavelength [nm]') + end + + + ylabel('Magnitude [dB]') + legend + grid minor; + + end + end end diff --git a/Datatypes/emlmodes.m b/Datatypes/eml_mode.m similarity index 78% rename from Datatypes/emlmodes.m rename to Datatypes/eml_mode.m index 9486815..5be8357 100644 --- a/Datatypes/emlmodes.m +++ b/Datatypes/eml_mode.m @@ -1,4 +1,4 @@ -classdef emlmodes < int32 +classdef eml_mode < int32 enumeration im_linear (1) diff --git a/Functions/applyPulseShaping.m b/Functions/applyPulseShaping.m index 3b76620..9a43655 100644 --- a/Functions/applyPulseShaping.m +++ b/Functions/applyPulseShaping.m @@ -13,8 +13,8 @@ else end %Bau das Filter (hier rrc) -racos_len = 2048; -alpha = 0.1; +racos_len = 32; +alpha = 0.05; h = rcosdesign(alpha,racos_len,sps); %Apply Filter using Matlab build in fctn. diff --git a/setup_simulation.m b/setup_simulation.m index b8a8366..3085f72 100644 --- a/setup_simulation.m +++ b/setup_simulation.m @@ -1,4 +1,6 @@ -clear all + + +%% Set Simulation Variables O = 18; %order of prbs N = 2^(O-1); %length of prbs @@ -9,7 +11,8 @@ M = 4; %PAM-M bitpattern = zeros(N,log2(M)); % Symbol Rate -fsym = 80e9; +fsym = 112e9; + % DAC Rate fdac = 120e9; @@ -22,144 +25,157 @@ fadc = 256e9; % Simulation frequency in "analog domain" fsimu = kover * fdac ; + + %% CONSTRUCT ALL CLASSES -pam_mapper = PAMmapper(M,0); +digimod = PAMmapper(M,0); -awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',80e9,'lpf_type',filtertypes.gaussian); +awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); -fil_50 = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); +lp_laser = Filter('filtdegree',1,"f_cutoff",60e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); - -u_pi = 3.5; +u_pi = 4.6; 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",10000); +extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",10000); -amp = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power"); +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); +fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",17,"lambda0",1550,"gamma",0); -phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); +optatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",5); -fil_diode = Filter('filtdegree',1,"f_cutoff",120e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); +edfaamp = Amplifier("amp_mode","edfa_increase_nase","gain_mode","output_power","amplification_db",0,"nase_mode","generate_ase","noifig",5); + +phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + +lp_diode = Filter('filtdegree',1,"f_cutoff",80e9,"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,"training_loops",5,... + "Ne",[50,0,0],"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",1,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); -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,"training_loops",5,... - "Ne",[50,0,0],"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",2,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); %% PROCESS -% INFORMATION SIGNAL + +% PRBS Generation for i = 1:log2(M) [bitpattern(:,i),seed] = prbs(O,N,seed); end -bits = Informationsignal(bitpattern); +% Build Inf. signal class +bits = Informationsignal(bitpattern); -mod_out = pam_mapper.map(bits); +% Digi Mod +mod_out = digimod.map(bits); -reference = mod_out; + +reference = mod_out; mod_out.signal = applyPulseShaping(mod_out.signal,fsym,fdac); -% ELECTRICAL DOMAIN -X = awg.process(mod_out); -X = fil_50.process(X); -X = X.normalize; -% OPTICAL DOMAIN -X = extmodlaser.process(X); +% AWG -> ELECTRICAL DOMAIN +X = awg.process(mod_out); -X = amp.process(X); +X = lp_laser.process(X); -X = fib.process(X); +X = X.normalize; -X = phdiode.process(X); +% Laser; Modulation -> OPTICAL DOMAIN +X = extmodlaser.process(X); +X.spectrum(fsimu,"displayname",'laser out','figurename','after laser'); -X = fil_diode.process(X); +%X = amp.process(X); -% ELECTRICAL DOMAIN -X = scp.process(X); +% Fiber Propagation +X = fib.process(X); +X = optatten.process(X); +X = edfaamp.process(X); -X = X.resample("fs_out",2*fsym,"fs_in",fadc); +% Photo Diode -> ELECTRICAL DOMAIN +X = phdiode.process(X); +X = lp_diode.process(X); +X.spectrum(fsimu,"displayname",'diode out','figurename','after diode'); +% Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) +X = scp.process(X); + +% Resample to Symbol Rate +X = X.resample("fs_out",2*fsym,"fs_in",fadc); % INFORMATION SIGNAL -X = X.normalize; -eq_out = eq.process(X,reference); -rx_series = eq_out.signal; +X = X.normalize; -demap_out = pam_mapper.demap(eq_out); +% Equalizer +eq_out = eq.process(X,reference); +rx_series = eq_out.signal; -[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(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:4) -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 - -a = fftshift(xcorr(X.signal(1,:),circshift(bitpattern(:,1)',0))); +% Digi Demod +demap_out = digimod.demap(eq_out); % BER -if (length(X.signal) ~= length(bitpattern')) - warning("TX and RX bitstreams have different length...") -end - - - -disp(demap_out.logbook); +[bits,errors,BER] = calc_ber(demap_out.signal(:,1:end-1),bitpattern(1:end,:)',0); disp(['BER: ', sprintf('%2E',BER)]); +% disp(demap_out.logbook); - - - - - - - +%% Generate some Plots +if 0 + col = cbrewer2('Paired',8); + + figure(22) + clf +% subplot(3,1,1) +% hold on +% plot(reference.signal(10000:end-20),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); +% plot(rx_series(10000:end-20),'DisplayName','Rx','Color',col(6,:)); +% title('Modulated Sequence Tx - Rx'); +% legend +% hold off + subplot(2,1,1) + hold on + plot(reference.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); + plot(rx_series(4150:4175),'DisplayName','Rx after EQ','Color',col(6,:),'LineWidth',1); + title('Modulated Sequence Zoom'); + legend + hold off + subplot(2,1,2) + hold on + stem(demap_out.signal(1,4150:4175),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) + stem(bitpattern(4150:4175,1)','DisplayName','Rx','Color',col(6,:),'LineStyle','--','LineWidth',2) + title('Bitpattern Tx - Rx'); + legend + hold off + + + xax = 1:eq_out.length; + + figure(12) + sgtitle('') + subplot(1,4,1:2) + scatter(1:4:X.length,X.signal(1:4:end),4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Before EQ'); + xlim([1, xax(end)]); + %ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + + subplot(1,4,3:4) + scatter(xax,eq_out.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ'); + xlim([1, xax(end)]); + %ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend +end