diff --git a/Classes/00_signals/Electricalsignal.m b/Classes/00_signals/Electricalsignal.m index 61b872c..2c067f2 100644 --- a/Classes/00_signals/Electricalsignal.m +++ b/Classes/00_signals/Electricalsignal.m @@ -33,11 +33,6 @@ classdef Electricalsignal < Signal end - function obj = normalize(obj) - - obj.signal = obj.signal/sqrt(mean(abs(obj.signal),"all")); - - end end diff --git a/Classes/00_signals/Informationsignal.m b/Classes/00_signals/Informationsignal.m index b3c7940..1e32d05 100644 --- a/Classes/00_signals/Informationsignal.m +++ b/Classes/00_signals/Informationsignal.m @@ -21,11 +21,6 @@ classdef Informationsignal < Signal end - function obj = normalize(obj) - - obj.signal = obj.signal/sqrt(mean(abs(obj.signal),"all")); - - end end end diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index 1fa5d94..d63c452 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -33,12 +33,6 @@ classdef Opticalsignal < Signal end - function obj = normalize(obj) - - obj.signal = obj.signal/sqrt(mean(abs(obj.signal).^2)); - - end - function pow = power(obj) pow = mean(abs(obj.signal.^2)) ; diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index e5b551c..9e9fc7e 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -106,8 +106,13 @@ classdef Signal end + %% Add signals from one signal to another, the first object will sustain + function Sum = plus(X,Y) + Sum = X; %first input object will sustain + Sum.signal = X.signal + Y.signal; + end - %% Display length + %% Display length function return_length = length(obj) %METHOD1 Summary of this method goes here % Detailed explanation goes here @@ -135,7 +140,7 @@ classdef Signal end - %% Resample Signal + %% Resample Signal function obj = resample(obj,options) arguments @@ -203,7 +208,7 @@ classdef Signal psd = psd/length(Fsignal); %smoothing - psd = smooth(psd,50); + psd = smooth(psd,1000); psd_plot = 20*log10(psd); @@ -255,6 +260,57 @@ classdef Signal end + %% + function obj = normalize(obj,options) + + arguments + obj Signal + options.mode normalization_mode = normalization_mode.rms + end + + switch options.mode + case normalization_mode.rms + obj.signal = obj.signal/sqrt(mean(abs(obj.signal).^2,"all")); + case normalization_mode.oneone + obj.signal = obj.signal/max(abs(obj.signal)); + end + + end + + %% + function obj = delay(obj,options) + + arguments + obj Opticalsignal + options.delay_meter double = 0 + end + + delay_t = options.delay_meter /(physconst("LightSpeed")/1.4677); + + delay_n = round(delay_t .* obj.fs); + + % build "long" hann window to fade the signal in and out + % -> prevent hard step in the signal! + hann_wind = hann(200); + ones_wind = ones(size(obj.signal)); + ones_wind(1:100) = hann_wind(1:100); + ones_wind(end-100:end) = hann_wind(end-100:end); + + % subtract average + mu = mean(obj.signal,"all"); + + obj.signal = obj.signal - mu; + + %apply hann + obj.signal = obj.signal .* ones_wind; + + %add average again + obj.signal = obj.signal + mu; + + % finally circshift the signal + obj.signal=circshift(obj.signal,delay_n); + + end end end diff --git a/Classes/01_transmit/AWG.m b/Classes/01_transmit/AWG.m index c1bf453..cb8f558 100644 --- a/Classes/01_transmit/AWG.m +++ b/Classes/01_transmit/AWG.m @@ -32,7 +32,7 @@ classdef AWG % Detailed explanation goes here arguments - options.preset = [] ; + options.preset = 'none'; options.kover = 16; options.repetitions = 1; options.normalize = 1; diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 462b61e..89bdb54 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -50,6 +50,8 @@ classdef PAMmapper pam_sig=2*pam_sig-3; end + pam_sig = pam_sig .* 1/sqrt(5); + case 3 % 8-ASK: @@ -86,8 +88,6 @@ classdef PAMmapper %28.03.2023 - Silas Oett. - Extracted from digi_demod.m % - obj.thresholds = 0; - switch log2(obj.M) case 1 @@ -106,6 +106,7 @@ classdef PAMmapper elseif obj.unipolar==1 thres=[0.5,1.5,2.5]; end + thres = thres .* 1/sqrt(5); case 3 @@ -174,6 +175,8 @@ classdef PAMmapper 1-comp_real(:,:,4)+comp_real(:,:,12)]; end + data_out = data_out'; + end diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index c04ea8e..aae5c49 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -41,7 +41,7 @@ classdef Pulseformer signalclass_in.signal = obj.process_(signalclass_in.signal); % append to logbook - lbdesc = ['Applied Pulseshaping']; + lbdesc = 'Applied Pulseshaping'; signalclass_in = signalclass_in.logbookentry(lbdesc); % write to output @@ -83,23 +83,56 @@ classdef Pulseformer if obj.pulseform == pulseform.rrc %Bau das Filter (hier rrc) - racos_len = obj.pulselength ; + racos_len = obj.pulselength*2 ; alpha = obj.rrcalpha; h = rcosdesign(alpha,racos_len,sps); + h = h./ max(h); end - %Apply Filter using Matlab build in fctn. - data_out = upfirdn(data_in,h,up,dn); + % Apply filter the long way (from move_it) + % block length in samples + + data_in = data_in'; + blen = length(data_in)*sps; + + % oversample symbol sequence + + symbolov=zeros(size(data_in,1),blen); + symbolov(:,1:sps:blen-sps+1)=data_in; - %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(data_in)*up/dn) -1); + H=fft(h,blen); + + % Convolution of Bit sequence with impulse response - data_out = data_out(st:en); + % cyclic convolution + data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; + + data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); + + if rem(obj.fdac,obj.fsym) + data_out = data_out(1:dn:end); %! + end + + +% %Apply Filter using Matlab build in fctn. +% data_out = upfirdn(data_in,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(data_in)*up/dn) -1); +% +% data_out = data_out(st:en); + + %scaling?! see pulsef module line 696 + scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part + data_out = data_out./scale; + + data_out = data_out'; %Check output integrity if round(up/dn * length(data_in)) ~= length(data_out) - warning('Check signal length after pulse shaping'); + %warning('Check signal length after pulse shaping'); + disp('Check signal length after pulse shaping'); end diff --git a/Classes/02_etc/Amplifier.m b/Classes/02_etc/Amplifier.m index c85a680..de3970c 100644 --- a/Classes/02_etc/Amplifier.m +++ b/Classes/02_etc/Amplifier.m @@ -149,6 +149,7 @@ classdef Amplifier end function nase_numeric = generateAseNoise(~, nase, fs, dimension) + rng(2023); nase_numeric = (randn(dimension) + 1i*randn(dimension))*sqrt(nase/2*fs) ; end diff --git a/Classes/02_etc/Filter.m b/Classes/02_etc/Filter.m index e328e68..2a0e8b7 100644 --- a/Classes/02_etc/Filter.m +++ b/Classes/02_etc/Filter.m @@ -43,14 +43,13 @@ classdef Filter function signalclass_out = process(obj,signalclass_in) % actual processing of the signal - signalclass_in.signal = obj.process_(signalclass_in.signal); + signalclass_out = signalclass_in; + signalclass_out.signal = obj.process_(signalclass_in.signal); % append to logbook filterdesc = [num2str(obj.filtdegree),'. order ',char(obj.filterType),' filter with f_cutoff at ', num2str(obj.f_cutoff*1e-9), ' GHz.']; - signalclass_in = signalclass_in.logbookentry(filterdesc); + signalclass_out = signalclass_out.logbookentry(filterdesc); - % write to output - signalclass_out = signalclass_in; end diff --git a/Classes/02_optical/Fiber.m b/Classes/02_optical/Fiber.m index bda66a6..b903d8c 100644 --- a/Classes/02_optical/Fiber.m +++ b/Classes/02_optical/Fiber.m @@ -23,7 +23,7 @@ classdef Fiber %FIBER Construct an instance of this class % Detailed explanation goes here arguments - options.fsimu + options.fsimu options.fiber_length = 0 options.alpha = 0.2 options.D = 17 @@ -42,22 +42,20 @@ classdef Fiber obj.gamma = options.gamma; obj.dphimax = options.dphimax; - obj.b2 = -obj.D*obj.lambda0^2/(2*pi*Constant.LightSpeed); - obj.b3 = ((obj.lambda0.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope); - obj.alpha_lin = obj.alpha/10*log(10)/1000; + end function signalclass_out = process(obj,signalclass_in) - % actual processing of the signal (steps 1. - 3.) - signalclass_in.signal = obj.process_(signalclass_in.signal); + % actual processing of the signal (steps 1. - 3.) + signalclass_in.signal = obj.process_(signalclass_in.signal); - % append to logbook - lbdesc = 'Fiber '; + % append to logbook + lbdesc = 'Fiber '; signalclass_in = signalclass_in.logbookentry(lbdesc); - % write to output + % write to output signalclass_out = signalclass_in; end @@ -66,6 +64,10 @@ classdef Fiber %METHOD1 Summary of this method goes here % Detailed explanation goes here + obj.b2 = -obj.D*obj.lambda0^2/(2*pi*Constant.LightSpeed); + obj.b3 = ((obj.lambda0.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope); + obj.alpha_lin = obj.alpha/10*log(10)/1000; + N = length(opt_in); faxis = linspace(-obj.fsimu/2,obj.fsimu/2,N+1); faxis = ifftshift(faxis(:,1:end-1)); @@ -109,7 +111,7 @@ classdef Fiber z_prop = z_prop + dz; maxPow = obj.gamma*max(abs(yout).^2); - Leff = obj.dphimax/maxPow; + Leff = obj.dphimax/maxPow; dz_new = Leff; diff --git a/Classes/03_receive/Scope.m b/Classes/03_receive/Scope.m index 23b1786..d7125f9 100644 --- a/Classes/03_receive/Scope.m +++ b/Classes/03_receive/Scope.m @@ -20,6 +20,8 @@ classdef Scope filtertype lpf_bw + block_dc + %during construction %during process @@ -47,24 +49,17 @@ classdef Scope options.filtertype = filtertypes.bessel_bilin; options.lpf_bw = 120e9; + + options.block_dc = 1; end - obj.fsimu = options.fsimu; - obj.fadc = options.fadc; - obj.adcresolution = options.adcresolution; - obj.quantbuffer = options.quantbuffer; + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end - obj.rand_samplingdelay = options.rand_samplingdelay; % use a randomized sample delay INSTEAD of samplingdelay - obj.samplingdelay = options.samplingdelay; %specifiy a sampling delay - - obj.freq_offset = options.freq_offset; %offset of the sampler - obj.samp_jitter = options.samp_jitter; %include jitter in [s] - - obj.fixed_delay = options.fixed_delay; %fix the delay of the filter or use minimal delay for kausal system - obj.delay = options.delay; %specify a fixed delay of the filter - - obj.filtertype = options.filtertype; - obj.lpf_bw = options.lpf_bw; end @@ -91,14 +86,17 @@ classdef Scope % Detailed explanation goes here % sample signal - % TODO: implement and test the delays. Also look for delay of - % lpf filter - %yout = obj.sampleSignal(xin); + % TODO: implement and test the delays. + + % resample yout = resample(xin,obj.fadc,obj.fsimu); % quantize signal yout = obj.quantize(yout); + if obj.block_dc + yout = yout-mean(yout,1); + end end diff --git a/Datatypes/normalization_mode.m b/Datatypes/normalization_mode.m new file mode 100644 index 0000000..768dd8f --- /dev/null +++ b/Datatypes/normalization_mode.m @@ -0,0 +1,8 @@ +classdef normalization_mode < int32 + + enumeration + rms (1) + oneone (2) + end + +end \ No newline at end of file diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index fe8d5e6..387ea9d 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -1,7 +1,7 @@ function [bits,errors,BER] = calc_ber(data_in,data_ref,skip) -data_ref=logical(data_ref); -data_in = logical(data_in); +data_ref=logical(data_ref)'; +data_in = logical(data_in)'; bits = 0; errors=0; @@ -10,18 +10,14 @@ 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); diff --git a/setup_simulation.m b/setup_simulation.m index b095486..125c2a3 100644 --- a/setup_simulation.m +++ b/setup_simulation.m @@ -1,6 +1,6 @@ -rng(2023); +rng(2020); %% Set Simulation Variables @@ -32,48 +32,50 @@ fsimu = kover * fdac ; digimod = PAMmapper(M,0); -pulsef = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.05); +pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.05); -awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); +awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); -lp_laser = Filter('filtdegree',1,"f_cutoff",60e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp); +lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); u_pi = 4.6; vbias = (0.5*u_pi)-u_pi; -extmodlaser = EML("mode",eml_mode.im_cosinus,"power",15,"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"); +extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1e6); fib = Fiber("fsimu",fdac*kover,"fiber_length",1,"alpha",0.2,"D",17,"lambda0",1550,"gamma",0); -optatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0); +reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-20); +reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",1,"alpha",0.2,"D",17,"lambda0",1550,"gamma",0); + + +opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0); 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); +lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"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); + "adcresolution",6,"quantbuffer",0.1,'block_dc',1); eq = EQ("K",2,"plottrain",0,"plotfinal",0,... "training_length",4096,"training_loops",5,... - "Ne",[30,0,0],"Nb",[0,0,0],... - "DCmu",0.005,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.002,"FFEmu",0.00,... + "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.00,... "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,"DB_aim",0); + "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",1,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); %% PROCESS -output_pow = 0 -for lp1 = 1:length(output_pow) +fiblen = 0; +for lp1 = 1:length(fiblen) % change parameters - optatten.amplification_db = output_pow(lp1); + fib.fiber_length = fiblen(lp1); % PRBS Generation @@ -93,25 +95,38 @@ for lp1 = 1:length(output_pow) reference = mod_out; % shape shape - mod_out = pulsef.process(mod_out); - test = applyPulseShaping(reference.signal,fsym,fdac); + X = pulseform.process(mod_out); + test = applyPulseShaping(reference.signal,fsym,fdac); % AWG -> ELECTRICAL DOMAIN - X = awg.process(mod_out); + awg_out = awg.process(X); + X = lp_laser.process(awg_out); X = lp_laser.process(X); - X = X.normalize; + X = X.normalize("mode","oneone"); + X.signal = X.signal .* 1.3800; % Laser; Modulation -> OPTICAL DOMAIN X = extmodlaser.process(X); - %X.spectrum(fsimu,"displayname",'laser out','figurename','after laser'); - - %X = amp.process(X); % Fiber Propagation X = fib.process(X); - X = optatten.process(X); + + %% Reflect with attenuation + R = reflectionpoint.process(X); + + % Propagate + R = reflectionprop.process(R); + + % Delay + R = R.delay("delay_meter",200); + + % Add together + X = X+R; + + %% + X = opticatten.process(X); X = edfaamp.process(X); % Photo Diode -> ELECTRICAL DOMAIN @@ -124,8 +139,8 @@ for lp1 = 1:length(output_pow) % Resample to Symbol Rate X = X.resample("fs_out",2*fsym,"fs_in",fadc); - % INFORMATION SIGNAL - X = X.normalize; + % Normalize + X = X.normalize("mode","rms"); % Equalizer eq_out = eq.process(X,reference); @@ -135,7 +150,8 @@ for lp1 = 1:length(output_pow) demap_out = digimod.demap(eq_out); % BER - [bits,errors,BER] = calc_ber(demap_out.signal(:,1:end-1),bitpattern(1:end,:)',0); + [bits,errors,BER] = calc_ber(demap_out.signal(1:end-2,:),bitpattern(1:end,:),0); + disp(['BER: ', sprintf('%2E',BER), ' ERRORS: ' ,num2str(sum(errors))]); bercurve(lp1) = BER; @@ -145,33 +161,13 @@ end - - - - - - - - - - - - - - %% Generate some Plots -if 0 +if 1 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); @@ -181,7 +177,7 @@ if 0 hold off subplot(2,1,2) hold on - stem(demap_out.signal(1,4150:4175),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) + stem(demap_out.signal(4150:4175,1),'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 @@ -193,7 +189,7 @@ if 0 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'); + scatter(1:4:X.length,X.signal(1:4:end),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Before EQ'); xlim([1, xax(end)]); %ylim([-2 2]); xlabel('Sampling Index')