diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index ec7c55e..279a913 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -51,7 +51,7 @@ classdef PAMsource warning('No Pulseformer given. Proceeding with RRC and alpha 0.05'); options.pulseformer = Pulseformer("fsym",obj.fsym,"fdac",obj.fs_out,"pulse","rrc","pulselength",16,"rrcalpha",0.05); end - + end @@ -67,12 +67,12 @@ classdef PAMsource if obj.useprbs for i = 1:log2(obj.M) [bitpattern(:,i),seed] = prbs(O,N,seed); - end + end else s = RandStream('twister','Seed',obj.randkey); for i = 1:log2(obj.M) bitpattern(:,i) = randi(s,[0 1], N, 1); - end + end end if obj.M == 6 @@ -102,7 +102,12 @@ classdef PAMsource %%%%% Hard clip digital signal to PAM range before DAC %%%%%% if obj.applyclipping - digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor); + try + digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor); + catch + digi_sig.signal(digi_sig.signal > sym_max * obj.clipfactor) = sym_max * obj.clipfactor; + digi_sig.signal(digi_sig.signal < sym_min * obj.clipfactor) = sym_min * obj.clipfactor; + end end % append to logbook diff --git a/Classes/03_receive/Photodiode.m b/Classes/03_receive/Photodiode.m index ecf9cac..2af5fb5 100644 --- a/Classes/03_receive/Photodiode.m +++ b/Classes/03_receive/Photodiode.m @@ -72,6 +72,7 @@ classdef Photodiode % Thermal Noise therm_current_psd = (2 * k * T / R ) ; %squared + therm_current_psd = 1.8e-11^2; Bw = obj.fsimu; therm_noise_pow = therm_current_psd * Bw; %squared diff --git a/Classes/04_DSP/FFE_DCremoval.m b/Classes/04_DSP/FFE_DCremoval.m new file mode 100644 index 0000000..7478bff --- /dev/null +++ b/Classes/04_DSP/FFE_DCremoval.m @@ -0,0 +1,179 @@ +classdef FFE_DCremoval < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + mu_dc + + constellation + + decide + end + + methods + function obj = FFE_DCremoval(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.mu_dc = 0.05; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + + arguments + obj + x + d + mio + epochs + N + training + showviz + end + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + % subplot(2,2,3) + % dplot = x(1:1+500); + % b = scatter(1:numel(dplot),dplot,5,'x'); + % xline(1) + % xline(obj.order) + % ylim([-3 3]) + % xlim([0 500]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + err = 0; + e_dc = 0; + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + U = x(obj.order+sample-1:-1:sample); + + y(symbol,1) = e_dc + obj.e.' * U; % Calculating output of LMS __ * | + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + % b.YData = x(symbol:symbol+500); + c.YData = obj.e; + drawnow; + end + + e_dc = e_dc - obj.mu_dc * err(symbol); + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + end + + end +end + diff --git a/Classes/04_DSP/FFE_FFDCAVG.m b/Classes/04_DSP/FFE_FFDCAVG.m new file mode 100644 index 0000000..1272e01 --- /dev/null +++ b/Classes/04_DSP/FFE_FFDCAVG.m @@ -0,0 +1,193 @@ +classdef FFE_FFDCAVG < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + mu_buff + + constellation + + decide + end + + methods + function obj = FFE_FFDCAVG(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.mu_buff = 0; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + + arguments + obj + x + d + mio + epochs + N + training + showviz + end + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + % subplot(2,2,3) + % dplot = x(1:1+500); + % b = scatter(1:numel(dplot),dplot,5,'x'); + % xline(1) + % xline(obj.order) + % ylim([-3 3]) + % xlim([0 500]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + + + for epoch = 1 : epochs + symbol = 0; + err_buffer = zeros(numel(obj.constellation),50); + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + U = x(obj.order+sample-1:-1:sample); + + y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * | + + if training + [~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + err_buffer(symbol_idx,1) = err(symbol); + err_buffer(symbol_idx,:) = circshift(err_buffer(symbol_idx,:),1); + + y(symbol) = y(symbol) - obj.mu_buff*mean(err_buffer(symbol_idx,:)); + + if training + [~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + % b.YData = x(symbol:symbol+500); + c.YData = obj.e; + drawnow; + end + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + end + + end +end + diff --git a/Classes/04_DSP/FFE_adaptive_decision.m b/Classes/04_DSP/FFE_adaptive_decision.m new file mode 100644 index 0000000..892305a --- /dev/null +++ b/Classes/04_DSP/FFE_adaptive_decision.m @@ -0,0 +1,149 @@ +classdef FFE_adaptive_decision < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + buffer_length + + constellation + + decide + end + + methods + function obj = FFE_adaptive_decision(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.buffer_length = 100; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + + arguments + obj + x + d + mio + epochs + N + training + showviz + end + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + for epoch = 1 : epochs + symbol = 0; + % y_buffer = zeros(numel(obj.constellation),500); + y_buffer = repmat(obj.constellation,1,obj.buffer_length); + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + U = x(obj.order+sample-1:-1:sample); + + y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * | + + if training + [~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = d(symbol); + else + y_buffer(y_buffer==0) = NaN; + adap_constellation = mean(y_buffer,2,"omitnan"); + [~,symbol_idx] = min(abs(y(symbol) - adap_constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + y_buffer(symbol_idx,1) = y(symbol); + y_buffer(symbol_idx,:) = circshift(y_buffer(symbol_idx,:),1); + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS + end + + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + end + + end +end + diff --git a/projects/MPI_August/EQ_development.m b/projects/MPI_August/EQ_development.m new file mode 100644 index 0000000..3cf1c3d --- /dev/null +++ b/projects/MPI_August/EQ_development.m @@ -0,0 +1,42 @@ + + +Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig"); +Scope_sig = Scope_sig.Scpe_sig; +Symbols = load("projects/MPI_August/symbols.mat","Symbols"); +Symbols = Symbols.Symbols; +Bits = load("projects/MPI_August/bits.mat","Bits"); +Bits = Bits.Bits; + +mdc = [1e-3,1e-2,2e-2,3e-2,4e-2,5e-2,1e-1]; +for m = 1:numel(mdc) + + Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_dc",mdc(m)); + [EQ_sig] = Eq.process(Scope_sig,Symbols); + + %%%%% DEMAP %%%%%% + Rx_bits = PAMmapper(4,0).demap(EQ_sig); + + % BER + [~,errors_bm,ber(m),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +end + +figure() +scatter(mdc,ber,20); +yline(3.8e-3,'DisplayName','HD-FEC'); +xlabel('Mu DC'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. Mu DC'); +set(gca,'yscale','log'); +set(gca,'xscale','log'); +grid on; +legend + + +%%%% Look at Pam levels %%%%% +if 1 + a = PAMmapper(4,0).separate_pamlevels(EQ_sig); + figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.'); +end + + + diff --git a/projects/MPI_August/EQ_development_ffe_adap_decision.m b/projects/MPI_August/EQ_development_ffe_adap_decision.m new file mode 100644 index 0000000..fc8fd2b --- /dev/null +++ b/projects/MPI_August/EQ_development_ffe_adap_decision.m @@ -0,0 +1,56 @@ + + + +Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig"); +Scope_sig = Scope_sig.Scpe_sig; +Symbols = load("projects/MPI_August/symbols.mat","Symbols"); +Symbols = Symbols.Symbols; +Bits = load("projects/MPI_August/bits.mat","Bits"); +Bits = Bits.Bits; + + +lenbuff = [80]; + +for m = 1:numel(lenbuff) + + Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"buffer_length",lenbuff(m)); + + [EQ_sig] = Eq.process(Scope_sig,Symbols); + + %%%%% DEMAP %%%%%% + Rx_bits = PAMmapper(4,0).demap(EQ_sig); + + % BER + [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + disp(['BER: ',sprintf('%.1E',ber)]); + +end + +% fs = 2*112e9; +% N = 2^(nextpow2(length(Eq.error(1,:)))-6); +% [p_lin,w] = pwelch(Eq.error(1,:),hanning(N),N/2,N,fs,"centered","power","mean"); +% p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power" + +% figure(123); % If figure does not exist, create new figure +% hold on +% plot(w.*1e-9,p_dbm,'DisplayName',['bla'],'LineWidth',1); +% xlabel("Frequency in GHz"); +% %ylabel("Power/frequency (dB/Hz)"); +% ylabel("Power (dBm)"); +% xlim([-fs/2 fs/2].*1e-9) +% edgetick = 2^(nextpow2(fs*1e-9)); +% xticks([-edgetick:16:edgetick]); +% xlim([-244, 244]) +% ylim([-120,-0]); +% yticks([-200:10:10]); +% legend + + +%%%% Look at Pam levels %%%%% +if 1 + a = PAMmapper(4,0).separate_pamlevels(EQ_sig); + figure(15);hold on;scatter(1:EQ_sig.length,a,1,'.'); +end + + diff --git a/projects/MPI_August/EQ_development_ffe_dc_avg.m b/projects/MPI_August/EQ_development_ffe_dc_avg.m new file mode 100644 index 0000000..7de32fd --- /dev/null +++ b/projects/MPI_August/EQ_development_ffe_dc_avg.m @@ -0,0 +1,56 @@ + + + +Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig"); +Scope_sig = Scope_sig.Scpe_sig; +Symbols = load("projects/MPI_August/symbols.mat","Symbols"); +Symbols = Symbols.Symbols; +Bits = load("projects/MPI_August/bits.mat","Bits"); +Bits = Bits.Bits; + + +mubuff = [0.7]; + +for m = 1:numel(mubuff) + + Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_buff",mubuff(m)); + + [EQ_sig] = Eq.process(Scope_sig,Symbols); + + %%%%% DEMAP %%%%%% + Rx_bits = PAMmapper(4,0).demap(EQ_sig); + + % BER + [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + disp(['BER: ',sprintf('%.1E',ber)]); + +end + +% fs = 2*112e9; +% N = 2^(nextpow2(length(Eq.error(1,:)))-6); +% [p_lin,w] = pwelch(Eq.error(1,:),hanning(N),N/2,N,fs,"centered","power","mean"); +% p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power" + +% figure(123); % If figure does not exist, create new figure +% hold on +% plot(w.*1e-9,p_dbm,'DisplayName',['bla'],'LineWidth',1); +% xlabel("Frequency in GHz"); +% %ylabel("Power/frequency (dB/Hz)"); +% ylabel("Power (dBm)"); +% xlim([-fs/2 fs/2].*1e-9) +% edgetick = 2^(nextpow2(fs*1e-9)); +% xticks([-edgetick:16:edgetick]); +% xlim([-244, 244]) +% ylim([-120,-0]); +% yticks([-200:10:10]); +% legend + + +%%%% Look at Pam levels %%%%% +if 1 + a = PAMmapper(4,0).separate_pamlevels(EQ_sig); + figure(15);hold on;scatter(1:EQ_sig.length,a,1,'.'); +end + + diff --git a/projects/MPI_August/auswertung/rop_curve.m b/projects/MPI_August/auswertung/rop_curve.m new file mode 100644 index 0000000..c768bd7 --- /dev/null +++ b/projects/MPI_August/auswertung/rop_curve.m @@ -0,0 +1,35 @@ +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sir = wh.parameter.sir.values(1); +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pn_key = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; + +cols = linspecer(4); + +cnt = 0; +for pnk = pn_key + + cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop); + ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,sir,laser_linewidth,pnk,rop); + ber_adapt = wh.getStoValue('ber_adapt',M,datarate,sir,laser_linewidth,pnk,rop); + ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop); + + % Create the initial plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rop,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:)); + plot(rop,ber_dcavg',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:)); + plot(rop,ber_adapt',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:)); + plot(rop,ber_derem',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber derem",'Color',cols(4,:)); + +end + +yline(3.8e-3,'DisplayName','HD-FEC'); +xlabel('Signal to Interference Ratio (dB)'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. ROP'); +set(gca,'yscale','log'); +grid on; +legend diff --git a/projects/MPI_August/bits.mat b/projects/MPI_August/bits.mat new file mode 100644 index 0000000..489efd4 Binary files /dev/null and b/projects/MPI_August/bits.mat differ diff --git a/projects/MPI_August/mpi_uncorrelated.m b/projects/MPI_August/mpi_uncorrelated.m index 716fe91..6b83e01 100644 --- a/projects/MPI_August/mpi_uncorrelated.m +++ b/projects/MPI_August/mpi_uncorrelated.m @@ -2,25 +2,21 @@ %% Parameter to simulate and save params = struct; -params.M = [4]; +params.M = [8]; params.datarate = [224]; -params.sir = [24]; %decibel = attenuation of interference path -params.laser_linewidth = [1e6]; -params.pn_key = [11]; -params.vbias_rel = [0.5]; -params.rop = -5; - -params.clipfactor = [1.5]; -params.rrcalpha = [0.1]; +params.sir = [30]; %decibel = attenuation of interference path +params.laser_linewidth = [1e6]; +params.pn_key = [1]; +params.rop = [-12:0]; name = ['wh_',strrep(num2str(now),'.','')]; wh = DataStorage(params); -wh.addStorage("ber"); -wh.addStorage("rop"); -wh.addStorage("txpapr"); -wh.addStorage("er"); +wh.addStorage("ber_ffe"); +wh.addStorage("ber_dcavg"); +wh.addStorage("ber_adapt"); +wh.addStorage("ber_dcrem"); %% Init Params link_length = 10000; %meter @@ -32,7 +28,9 @@ disp(['Start Simulation of ',num2str(endcnt),' loops...']) tic for M = wh.parameter.M.values for datarate = wh.parameter.datarate.values - for rrcalpha = wh.parameter.rrcalpha.values + for pn_key = wh.parameter.pn_key.values + + %% SETUP HERE: %% kover = 8; @@ -44,157 +42,150 @@ for M = wh.parameter.M.values % MAIN SIGNAL %%%%% Symbol Generation %%%%%% - [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",2).process(); + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key).process(); %%%%% AWG %%%%%% El_sig = M8199.process(Digi_sig); El_sig = El_sig.*0.7222; - El_sig.signal = awgn(El_sig.signal,20,'measured',1); + El_sig.signal = awgn(El_sig.signal,20,'measured',pn_key); %%%%% Lowpass before Modulator %%%%%% - El_sig = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig); - + El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig); % INTERFERENCE SIGNAL %%%%% Symbol Generation %%%%%% - [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",1).process(); + [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key*2).process(); %%%%% AWG %%%%%% El_sig_i = M8199.process(Digi_sig_i); El_sig_i = El_sig_i.*0.7222; - El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',2); + El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',pn_key*2); %%%%% Lowpass before Modulator %%%%%% - El_sig_i = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); - El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig_i); + El_sig_i = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); + El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig_i); for laser_linewidth = wh.parameter.laser_linewidth.values - for pn_key = wh.parameter.pn_key.values - for vbias_rel = wh.parameter.vbias_rel.values + + % MAIN SIGNAL + %%%%% MODULATE E/O CONVERSION %%%%%% + vbias_rel = 0.5; + u_pi = 2.9; + vbias = -vbias_rel*u_pi; + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig); + Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + Opt_sig = Optfilter.process(Opt_sig); + er = Opt_sig.extinctionratio(fsym,M); + + + % INTERFERENCE SIGNAL + %%%%% MODULATE E/O CONVERSION %%%%%% + u_pi = 2.9; + vbias = -vbias_rel*u_pi; + [Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i); + Opt_sig_i = Optfilter.process(Opt_sig_i); - % MAIN SIGNAL - %%%%% MODULATE E/O CONVERSION %%%%%% - u_pi = 2.9; - vbias = -vbias_rel*u_pi; - [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig); - Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); - Opt_sig = Optfilter.process(Opt_sig); + j_ = wh.parameter.sir.length; + i_ = wh.parameter.rop.length; + ber=zeros(j_,i_); + patten=zeros(j_,i_); + + for j = 1:j_ + sir = wh.parameter.sir.values(j); - % INTERFERENCE SIGNAL - %%%%% MODULATE E/O CONVERSION %%%%%% - u_pi = 2.9; - vbias = -vbias_rel*u_pi; - [Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i); - Opt_sig_i = Optfilter.process(Opt_sig_i); + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig_i = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i); + + Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i); + + %%%%% ADD Interference and Main Signal %%%%%% + Opt_sig = Opt_sig_i + Opt_sig; + + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + % % MPI Channel + % Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir); - j_ = wh.parameter.sir.length; - i_ = wh.parameter.rop.length; - ber=zeros(j_,i_); - patten=zeros(j_,i_); + % Receiver ROP curve + parfor i = 1:i_ + rop=wh.parameter.rop.values(i); - for j = 1:j_ - sir = wh.parameter.sir.values(j); + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + patten(j,i) = Rx_sig.power; + + %%%%%% Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); + + %%%%%% Lowpass PhDiode %%%%%% + Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); + % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); - %%%%% Interference Signal Fiber Prop %%%%%% - Opt_sig_i = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i); + %%%%%% Scope %%%%%% + fadc = 256e9; + Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); - Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i); + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); - %%%%% ADD Interference and Main Signal %%%%%% - Opt_sig = Opt_sig_i + Opt_sig; + %%%%% EQUALIZE %%%%%% + Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_ffe(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_ffe(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); - %%%%% Interference Signal Fiber Prop %%%%%% - Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_buff",0.7); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_dcavg(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_dcavg(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"buffer_length",85); [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_adapt(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_adapt(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); - % % MPI Channel - % Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir); - - % Receiver ROP curve - for i = 1:i_ - rop=wh.parameter.rop.values(i); - - % Set ROP - Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); - patten(j,i) = Rx_sig.power; - - %%%%%% Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); - - %%%%%% Lowpass PhDiode %%%%%% - Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); - % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); - - - %%%%%% Scope %%%%%% - fadc = 256e9; - Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); - Scpe_sig.plot("displayname","SIgnal after Scope","fignum",999); - - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); - - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); - - %%%%% EQUALIZE %%%%%% - % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); - - % Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",25,"dfe_order",2,"sps",2,"decide",1); - - % Eq = EQ("Ne",[25,3,3],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",1); - - Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,3,3],"sps",2,"decide",0); - [EQ_sig] = Eq.process(Scpe_sig,Symbols); - - % EQ_sig.normalize("mode","rms").plot('fignum',23,'displayname','before eq') - - - %%%%% DEMAP %%%%%% - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - - %%%% Look at Pam levels %%%%% - if 1 - a = PAMmapper(M,0).separate_pamlevels(EQ_sig); - figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.'); - end - - % BER - [~,errors_bm,ber(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - cnt = cnt+1; - disp(['BER: ',sprintf('%.1E',ber(j,i)),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); - - - end - end - - for j = 1:j_ - sir = wh.parameter.sir.values(j); - for i = 1:i_ - rop=wh.parameter.rop.values(i); - wh.addValueToStorage(ber(j,i),'ber',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); - wh.addValueToStorage(patten(j,i),'rop',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); - wh.addValueToStorage(er,'er',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); - end - - end - + Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_dc",0.07); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + [~,errors_bm,ber_dcrem(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + disp(['BER: ',sprintf('%.1E',ber_dcrem(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); end - toc - disp(['Simulated: ',num2str(cnt/endcnt*100),' %']); - save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh"); end + + for j = 1:j_ + sir = wh.parameter.sir.values(j); + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + + wh.addValueToStorage(ber_ffe(j,i) ,'ber_ffe',M,datarate,sir,laser_linewidth,pn_key,rop); + wh.addValueToStorage(ber_dcavg(j,i),'ber_dcavg',M,datarate,sir,laser_linewidth,pn_key,rop); + wh.addValueToStorage(ber_adapt(j,i),'ber_adapt',M,datarate,sir,laser_linewidth,pn_key,rop); + wh.addValueToStorage(ber_dcrem(j,i),'ber_dcrem',M,datarate,sir,laser_linewidth,pn_key,rop); + + end + + end + + + toc + disp(['Simulated: ',num2str(cnt/endcnt*100),' %']); + save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh"); end end end @@ -202,5 +193,43 @@ end +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sir = wh.parameter.sir.values(1); +laser_linewidth = wh.parameter.laser_linewidth.values(1); +pn_key = wh.parameter.pn_key.values; +rop = wh.parameter.rop.values; + +cols = linspecer(4); + +cnt = 0; +for pnk = pn_key + + cnt = cnt+1; + ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop); + ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,sir,laser_linewidth,pnk,rop); + ber_adapt = wh.getStoValue('ber_adapt',M,datarate,sir,laser_linewidth,pnk,rop); + ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop); + + % Create the initial plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rop,ber_ffe',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:)); + plot(rop,ber_dcavg',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:)); + plot(rop,ber_adapt',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:)); + plot(rop,ber_derem',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber derem",'Color',cols(4,:)); + +end + +yline(3.8e-3,'DisplayName','HD-FEC'); +xlabel('Signal to Interference Ratio (dB)'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. SIR (MPI)'); +set(gca,'yscale','log'); +grid on; +legend + + + diff --git a/projects/MPI_August/scpe_sig.mat b/projects/MPI_August/scpe_sig.mat new file mode 100644 index 0000000..a7eaba8 Binary files /dev/null and b/projects/MPI_August/scpe_sig.mat differ diff --git a/projects/MPI_August/symbols.mat b/projects/MPI_August/symbols.mat new file mode 100644 index 0000000..2cf1201 Binary files /dev/null and b/projects/MPI_August/symbols.mat differ